Epsy587: Random Intercept Models in R
Carolyn J. Anderson
12/29/2022
- Set Up
- Descriptive Statistics
- Graphing Data
- Create some variables
- Model Fitting
- Model 4: center ses & add back in mean ses
- Model 5: Lots of level 1 variables
- Model 6: Lots of level 1 & 2 variables
- Model 7: Kitchen Sink
- Quick comparison
In this document, analyses will be done using High School and Beyond data set. These analyses are reported in the lecture notes on Random Intercept Models.
Set Up
Libraries
The major R package that we will use is lme4. We will use additional packages later in the semester, but lme4 is sufficient for most of what we'll do in the class.
The output will not include degrees of freedom, test statistics, and p-values for tests of fixed effects, we will also use the package lmerTest.
The package texreg will be used here to obtain a way to compare models. An alternative to texreg is stargazer.
The package lattice is required for some of the graphics.
library(lme4)
library(lmerTest)
library(texreg)
library(stargazer)
library(lattice)I will also set my working directory here
setwd('D:\\Dropbox\\edps587\\lectures\\3 randomintercept')Read in Data
In HSB data for this example are contained in two files
- HSB1data.txt that contains level 1 (tudent level) data.
- HSB2data.txt contains level 2 (school level) data.
These two files will be read in and then merged.
hsb1
hsb1 <- read.table(file="HSB1data.txt", header=TRUE)
# variable names
names(hsb1)## [1] "id" "minority" "female" "ses" "mathach"# dimensions of data frame
dim(hsb1)## [1] 7185 5# Top & bottom of the file
head(hsb1)## id minority female ses mathach
## 1 1224 0 1 -1.528 5.876
## 2 1224 0 1 -0.588 19.708
## 3 1224 0 0 -0.528 20.349
## 4 1224 0 0 -0.668 8.781
## 5 1224 0 0 -0.158 17.898
## 6 1224 0 0 0.022 4.583tail(hsb1)## id minority female ses mathach
## 7180 9586 1 1 1.612 20.967
## 7181 9586 0 1 1.512 20.402
## 7182 9586 0 1 -0.038 14.794
## 7183 9586 0 1 1.332 19.641
## 7184 9586 0 1 -0.008 16.241
## 7185 9586 0 1 0.792 22.733hsb2
hsb2 <- read.table(file="HSB2data.txt", header=TRUE)
names(hsb2)## [1] "id" "size" "sector" "pracad" "disclim" "himinty" "meanses"dim(hsb2)## [1] 160 7# Top & bottom of the file
head(hsb2)## id size sector pracad disclim himinty meanses
## 1 1224 842 0 0.35 1.597 0 -0.428
## 2 1288 1855 0 0.27 0.174 0 0.128
## 3 1296 1719 0 0.32 -0.137 1 -0.420
## 4 1308 716 1 0.96 -0.622 0 0.534
## 5 1317 455 1 0.95 -1.694 1 0.351
## 6 1358 1430 0 0.25 1.535 0 -0.014tail(hsb2)## id size sector pracad disclim himinty meanses
## 155 9347 1067 1 0.58 -0.905 0 0.218
## 156 9359 1184 1 0.69 -0.475 0 0.360
## 157 9397 1314 0 0.44 -0.231 0 0.140
## 158 9508 1119 1 0.52 -1.138 0 -0.132
## 159 9550 1532 0 0.45 0.791 0 0.059
## 160 9586 262 1 1.00 -2.416 0 0.627Merging hsb1 and hsb2
We will make a long file.
I know that the data sets are already sorted by the school id, which will be the matching variable.
hsb <- merge(hsb1,hsb2, by=c('id'))
dim(hsb)## [1] 7185 11names(hsb)## [1] "id" "minority" "female" "ses" "mathach" "size"
## [7] "sector" "pracad" "disclim" "himinty" "meanses"Descriptive Statistics
The number of schools can be found in a number of ways, for example
(nschools <- nrow(hsb2))## [1] 160#or
(nschool2 <- length(unique(hsb$id)))## [1] 160Number of students per school,
table(hsb$id)##
## 1224 1288 1296 1308 1317 1358 1374 1433 1436 1461 1462 1477 1499 1637 1906 1909
## 47 25 48 20 48 30 28 35 44 33 57 62 53 27 53 28
## 1942 1946 2030 2208 2277 2305 2336 2458 2467 2526 2626 2629 2639 2651 2655 2658
## 29 39 47 60 61 67 47 57 52 57 38 57 42 38 52 45
## 2755 2768 2771 2818 2917 2990 2995 3013 3020 3039 3088 3152 3332 3351 3377 3427
## 47 25 55 42 43 48 46 53 59 21 39 52 38 39 45 49
## 3498 3499 3533 3610 3657 3688 3705 3716 3838 3881 3967 3992 3999 4042 4173 4223
## 53 38 48 64 51 43 45 41 54 41 52 53 46 64 44 45
## 4253 4292 4325 4350 4383 4410 4420 4458 4511 4523 4530 4642 4868 4931 5192 5404
## 58 65 53 33 25 41 32 48 58 47 63 61 34 58 28 57
## 5619 5640 5650 5667 5720 5761 5762 5783 5815 5819 5838 5937 6074 6089 6144 6170
## 66 57 45 61 53 52 37 29 25 50 31 29 56 33 43 21
## 6291 6366 6397 6415 6443 6464 6469 6484 6578 6600 6808 6816 6897 6990 7011 7101
## 35 58 60 54 30 29 57 35 56 56 44 55 49 53 33 28
## 7172 7232 7276 7332 7341 7342 7345 7364 7635 7688 7697 7734 7890 7919 8009 8150
## 44 52 53 48 51 58 56 44 51 54 32 22 51 37 47 44
## 8165 8175 8188 8193 8202 8357 8367 8477 8531 8627 8628 8707 8775 8800 8854 8857
## 49 33 30 43 35 27 14 37 41 53 61 48 48 32 32 64
## 8874 8946 8983 9021 9104 9158 9198 9225 9292 9340 9347 9359 9397 9508 9550 9586
## 36 58 51 56 55 53 31 36 19 29 57 53 47 35 29 59And for a summary table of descriptive statistics ingoring clustering,
stargazer(hsb, type="html")##
## <table style="text-align:center"><tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Pctl(25)</td><td>Pctl(75)</td><td>Max</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">id</td><td>7,185</td><td>5,277.898</td><td>2,499.578</td><td>1,224</td><td>3,020</td><td>7,342</td><td>9,586</td></tr>
## <tr><td style="text-align:left">minority</td><td>7,185</td><td>0.275</td><td>0.446</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>
## <tr><td style="text-align:left">female</td><td>7,185</td><td>0.528</td><td>0.499</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>
## <tr><td style="text-align:left">ses</td><td>7,185</td><td>0.0001</td><td>0.779</td><td>-3.758</td><td>-0.538</td><td>0.602</td><td>2.692</td></tr>
## <tr><td style="text-align:left">mathach</td><td>7,185</td><td>12.748</td><td>6.878</td><td>-2.832</td><td>7.275</td><td>18.317</td><td>24.993</td></tr>
## <tr><td style="text-align:left">size</td><td>7,185</td><td>1,056.862</td><td>604.172</td><td>100</td><td>565</td><td>1,436</td><td>2,713</td></tr>
## <tr><td style="text-align:left">sector</td><td>7,185</td><td>0.493</td><td>0.500</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>
## <tr><td style="text-align:left">pracad</td><td>7,185</td><td>0.534</td><td>0.251</td><td>0.000</td><td>0.320</td><td>0.700</td><td>1.000</td></tr>
## <tr><td style="text-align:left">disclim</td><td>7,185</td><td>-0.132</td><td>0.944</td><td>-2.416</td><td>-0.817</td><td>0.460</td><td>2.756</td></tr>
## <tr><td style="text-align:left">himinty</td><td>7,185</td><td>0.280</td><td>0.449</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>
## <tr><td style="text-align:left">meanses</td><td>7,185</td><td>0.006</td><td>0.414</td><td>-1.188</td><td>-0.317</td><td>0.333</td><td>0.831</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr></table>Summary statistics for math by school,
mmean <- aggregate(mathach~id, data=hsb, FUN="mean")
msd <- aggregate(mathach~id, data=hsb, sd)
mmin <- aggregate(mathach~id, data=hsb, min)
mmax <- aggregate(mathach~id, data=hsb, max)
math.stat <- cbind(mmean,msd[,2], mmin[,2], mmax[,2])
names(math.stat) <- c('Mean Math', 'SD math', 'Min', 'Max')
head(math.stat)## Mean Math SD math Min Max NA
## 1 1224 9.715447 7.592785 -2.832 23.584
## 2 1288 13.510800 7.021843 1.575 23.578
## 3 1296 7.635958 5.351070 -1.353 23.172
## 4 1308 16.255500 6.114241 2.512 24.993
## 5 1317 13.177687 5.462586 3.220 23.736
## 6 1358 11.206233 5.876480 -1.347 21.142Graphing Data
Below are a number of graphs that we looked at in class.
Plot of mean ses by id --> variability between schools
plot(hsb$id, hsb$meanses, type='p', pch=21,
main='HSB: Group Mean SES by School ID',
xlab='School ID',
ylab='Group Mean SES')
Plot of student ses by id --> variability within schools
plot(hsb$id, hsb$ses, type='p', pch=21,
main='HSB: Student SES by School ID',
xlab='School ID',
ylab='Group Mean SES')
Just one school
school.1224 <- hsb[ which(hsb$id==1224), ]
plot(school.1224$ses, school.1224$mathach, type='p', pch=21,
main="HSB: Math by SES for School 1224",
xlab='Student SES',
ylab='Math Achievement')
reg.1224 <- lm(school.1224$mathach~school.1224$ses, data=school.1224)
abline(reg.1224)
Panel plot of a random sample of 20
Will need this for graph
hsb$cses <- hsb$ses - hsb$meanses
# Random sample of 20 schools
groups <- unique(hsb$id)[sample(1:160,20)]
subset <- hsb[hsb$id%in%groups,]
xyplot(mathach ~ ses | id, data=subset, col.line='black',
type=c('p'),
main='Varability in Math ~ SES relationship')
Panel plot of a random sample of 20 with regression lines
xyplot(mathach ~ ses | id, data=subset, col.line='black',
type=c('p','r'),
main='Varability in Math ~ SES relationship')
Panel using school centered ses
xyplot(mathach ~ cses | id, data=subset, col.line='black',
type=c('p','r'),
main='Math Achievement by School Centered SES',
xlab='School Centered SES',
ylab='Math Achievement Test Scores')
All regression lines in one plot
hsb$schid <- as.factor(hsb$id)
mathbysch <- lm(formula = mathach ~ 1*schid + ses*schid, data = hsb)
summary(mathbysch)##
## Call:
## lm(formula = mathach ~ 1 * schid + ses * schid, data = hsb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.2662 -4.3241 0.1164 4.4434 18.3617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.80513 1.07892 10.015 < 2e-16 ***
## ses 2.50858 1.42433 1.761 0.078243 .
## schid1288 2.30981 1.63811 1.410 0.158571
## schid1296 -2.71135 1.50564 -1.801 0.071779 .
## schid1308 5.38383 2.31095 2.330 0.019850 *
## schid1317 1.93263 1.49342 1.294 0.195675
## schid1358 0.50077 1.54571 0.324 0.745966
## schid1374 -1.02794 1.57352 -0.653 0.513602
## schid1433 7.59375 1.93963 3.915 9.13e-05 ***
## schid1436 6.40551 1.68127 3.810 0.000140 ***
## schid1461 1.79224 1.83710 0.976 0.329307
## schid1462 -0.86438 1.63846 -0.528 0.597825
## schid1477 3.22695 1.33789 2.412 0.015892 *
## schid1499 -1.45126 1.46669 -0.989 0.322463
## schid1637 -1.58240 1.93355 -0.818 0.413161
## schid1906 4.08035 1.53228 2.663 0.007765 **
## schid1909 2.91000 1.62495 1.791 0.073366 .
## schid1942 7.24481 2.15312 3.365 0.000770 ***
## schid1946 2.08878 1.45109 1.439 0.150067
## schid2030 0.80939 1.46000 0.554 0.579340
## schid2208 3.48380 1.44486 2.411 0.015928 *
## schid2277 -3.04284 1.60341 -1.898 0.057772 .
## schid2305 -0.15854 1.49082 -0.106 0.915314
## schid2336 4.87138 1.52493 3.194 0.001407 **
## schid2458 2.50705 1.37360 1.825 0.068019 .
## schid2467 0.37680 1.47283 0.256 0.798088
## schid2526 6.19572 1.41288 4.385 1.18e-05 ***
## schid2626 2.85731 1.46414 1.952 0.051035 .
## schid2629 4.13325 1.35396 3.053 0.002276 **
## schid2639 -4.79735 2.05301 -2.337 0.019482 *
## schid2651 -0.02971 1.46288 -0.020 0.983795
## schid2655 5.20761 1.66503 3.128 0.001770 **
## schid2658 1.43796 1.53994 0.934 0.350452
## schid2755 5.35972 1.60260 3.344 0.000829 ***
## schid2768 0.27005 1.62429 0.166 0.867961
## schid2771 2.48783 1.45899 1.705 0.088207 .
## schid2818 2.76317 1.43728 1.923 0.054584 .
## schid2917 -1.91950 1.64127 -1.170 0.242235
## schid2990 6.77560 1.63697 4.139 3.53e-05 ***
## schid2995 -0.82690 1.43411 -0.577 0.564232
## schid3013 1.97552 1.36488 1.447 0.147834
## schid3020 3.20957 1.36566 2.350 0.018791 *
## schid3039 5.13100 1.83264 2.800 0.005128 **
## schid3088 -0.81267 1.56344 -0.520 0.603223
## schid3152 2.31798 1.36792 1.695 0.090209 .
## schid3332 2.23115 1.79822 1.241 0.214740
## schid3351 -0.14055 1.50037 -0.094 0.925370
## schid3377 -2.05558 1.60519 -1.281 0.200385
## schid3427 8.98516 1.39791 6.428 1.38e-10 ***
## schid3498 5.67178 1.66749 3.401 0.000674 ***
## schid3499 2.02493 1.62384 1.247 0.212441
## schid3533 -0.43801 1.40641 -0.311 0.755476
## schid3610 4.19475 1.32624 3.163 0.001569 **
## schid3657 1.14131 1.54740 0.738 0.460804
## schid3688 3.22872 1.58226 2.041 0.041331 *
## schid3705 -0.74530 1.44138 -0.517 0.605118
## schid3716 2.03032 1.51871 1.337 0.181309
## schid3838 5.17079 1.36976 3.775 0.000161 ***
## schid3881 0.83901 1.45328 0.577 0.563741
## schid3967 1.97532 1.40096 1.410 0.158593
## schid3992 3.64354 1.45401 2.506 0.012238 *
## schid3999 0.46366 1.40372 0.330 0.741177
## schid4042 2.82945 1.40445 2.015 0.043982 *
## schid4173 1.81183 1.41450 1.281 0.200274
## schid4223 4.05123 1.41427 2.865 0.004189 **
## schid4253 -1.54158 1.40816 -1.095 0.273668
## schid4292 1.98114 1.43138 1.384 0.166379
## schid4325 2.57028 1.36363 1.885 0.059488 .
## schid4350 0.38311 1.52841 0.251 0.802084
## schid4383 0.20816 1.62999 0.128 0.898386
## schid4410 2.40179 1.44299 1.664 0.096067 .
## schid4420 3.73250 1.56158 2.390 0.016865 *
## schid4458 -3.80592 1.99749 -1.905 0.056777 .
## schid4511 2.60846 1.34873 1.934 0.053152 .
## schid4523 -2.32543 1.39637 -1.665 0.095889 .
## schid4530 -0.76610 1.51463 -0.506 0.613011
## schid4642 3.41613 1.33581 2.557 0.010569 *
## schid4868 1.04048 1.59120 0.654 0.513203
## schid4931 2.65308 1.40754 1.885 0.059486 .
## schid5192 -0.29347 1.57613 -0.186 0.852298
## schid5404 3.60813 1.88843 1.911 0.056091 .
## schid5619 2.40119 1.41430 1.698 0.089591 .
## schid5640 3.03094 1.36691 2.217 0.026631 *
## schid5650 3.45313 1.40740 2.454 0.014170 *
## schid5667 1.49680 1.42701 1.049 0.294258
## schid5720 3.39685 1.36331 2.492 0.012739 *
## schid5761 1.33681 1.42067 0.941 0.346753
## schid5762 -7.69105 2.76223 -2.784 0.005378 **
## schid5783 1.80637 1.58023 1.143 0.253033
## schid5815 -1.48153 2.19088 -0.676 0.498920
## schid5819 0.97556 1.40243 0.696 0.486692
## schid5838 2.59445 1.55763 1.666 0.095831 .
## schid5937 5.29408 2.03331 2.604 0.009243 **
## schid6074 3.39713 1.39636 2.433 0.015006 *
## schid6089 4.61900 1.51455 3.050 0.002299 **
## schid6144 -1.04797 1.53047 -0.685 0.493533
## schid6170 4.82816 1.77779 2.716 0.006628 **
## schid6291 1.37906 1.73028 0.797 0.425470
## schid6366 4.38120 1.39244 3.146 0.001660 **
## schid6397 2.58829 1.35153 1.915 0.055523 .
## schid6415 1.71611 1.37690 1.246 0.212677
## schid6443 -1.59200 1.62905 -0.977 0.328476
## schid6464 -3.03487 1.98818 -1.526 0.126942
## schid6469 6.32969 1.65971 3.814 0.000138 ***
## schid6484 2.21940 1.51317 1.467 0.142497
## schid6578 2.42375 1.50858 1.607 0.108178
## schid6600 1.15061 1.35044 0.852 0.394230
## schid6808 -1.30010 1.41807 -0.917 0.359275
## schid6816 3.01764 1.48769 2.028 0.042558 *
## schid6897 3.04094 1.44295 2.107 0.035115 *
## schid6990 -4.36426 1.44500 -3.020 0.002535 **
## schid7011 3.16745 1.50988 2.098 0.035958 *
## schid7101 1.03573 1.57341 0.658 0.510386
## schid7172 -2.45010 1.46812 -1.669 0.095190 .
## schid7232 2.18822 1.37403 1.593 0.111305
## schid7276 1.57909 1.36522 1.157 0.247455
## schid7332 3.09940 1.44341 2.147 0.031806 *
## schid7341 -0.73142 1.38249 -0.529 0.596784
## schid7342 0.81469 1.48396 0.549 0.583025
## schid7345 0.39338 1.34940 0.292 0.770663
## schid7364 3.39019 1.42314 2.382 0.017237 *
## schid7635 3.66547 1.40684 2.605 0.009195 **
## schid7688 7.59556 1.38536 5.483 4.34e-08 ***
## schid7697 4.10672 1.58794 2.586 0.009725 **
## schid7734 2.58722 1.82176 1.420 0.155601
## schid7890 -2.80691 1.56661 -1.792 0.073223 .
## schid7919 2.21900 1.70243 1.303 0.192470
## schid8009 2.47684 1.61701 1.532 0.125633
## schid8150 4.10403 1.48732 2.759 0.005807 **
## schid8165 5.08931 1.44793 3.515 0.000443 ***
## schid8175 1.16970 1.53199 0.764 0.445180
## schid8188 1.32457 1.56238 0.848 0.396586
## schid8193 5.83953 1.44927 4.029 5.65e-05 ***
## schid8202 0.58965 1.49306 0.395 0.692910
## schid8357 3.86372 1.59358 2.425 0.015353 *
## schid8367 -6.25875 1.94684 -3.215 0.001311 **
## schid8477 2.46471 1.48900 1.655 0.097915 .
## schid8531 1.36939 1.54514 0.886 0.375511
## schid8627 -0.11740 1.36836 -0.086 0.931630
## schid8628 5.89540 1.34224 4.392 1.14e-05 ***
## schid8707 1.55269 1.39934 1.110 0.267214
## schid8775 -0.99797 1.45950 -0.684 0.494141
## schid8800 -1.64775 1.82633 -0.902 0.366971
## schid8854 -5.09813 1.83355 -2.780 0.005443 **
## schid8857 4.07693 1.51268 2.695 0.007052 **
## schid8874 2.64583 1.55668 1.700 0.089241 .
## schid8946 0.16407 1.39210 0.118 0.906183
## schid8983 0.59417 1.42824 0.416 0.677410
## schid9021 2.30978 1.60740 1.437 0.150772
## schid9104 4.91627 1.92137 2.559 0.010527 *
## schid9158 -0.70717 1.43587 -0.493 0.622378
## schid9198 7.00277 1.74567 4.012 6.10e-05 ***
## schid9225 3.13095 1.51963 2.060 0.039404 *
## schid9292 -0.07461 2.40088 -0.031 0.975211
## schid9340 1.68832 1.75011 0.965 0.334734
## schid9347 2.16513 1.36774 1.583 0.113468
## schid9359 4.76061 1.45830 3.264 0.001102 **
## schid9397 -0.77780 1.40874 -0.552 0.580881
## schid9508 3.31515 1.50878 2.197 0.028037 *
## schid9550 0.07760 1.56086 0.050 0.960350
## schid9586 3.01995 1.57368 1.919 0.055021 .
## schid1288:ses 0.74687 2.33332 0.320 0.748912
## schid1296:ses -1.43262 1.97356 -0.726 0.467919
## schid1308:ses -2.38256 3.22857 -0.738 0.460564
## schid1317:ses -1.23467 2.13415 -0.579 0.562925
## schid1358:ses 2.55943 2.22678 1.149 0.250437
## schid1374:ses 1.34574 2.19366 0.613 0.539586
## schid1433:ses -0.65429 2.25481 -0.290 0.771691
## schid1436:ses -0.90802 2.15450 -0.421 0.673438
## schid1461:ses 3.75792 2.10272 1.787 0.073955 .
## schid1462:ses -3.33739 1.99605 -1.672 0.094570 .
## schid1477:ses -1.27798 1.82982 -0.698 0.484941
## schid1499:ses 1.12615 1.83964 0.612 0.540452
## schid1637:ses 0.60822 2.11412 0.288 0.773587
## schid1906:ses -0.36308 1.97595 -0.184 0.854217
## schid1909:ses 0.34621 2.17007 0.160 0.873250
## schid1942:ses -2.41920 2.60206 -0.930 0.352547
## schid1946:ses 1.07725 2.00262 0.538 0.590650
## schid2030:ses -1.09660 1.93800 -0.566 0.571520
## schid2208:ses 0.12806 1.94124 0.066 0.947406
## schid2277:ses -4.52361 1.84800 -2.448 0.014396 *
## schid2305:ses -3.29069 1.82288 -1.805 0.071085 .
## schid2336:ses -0.60361 1.99445 -0.303 0.762170
## schid2458:ses 0.44811 1.88184 0.238 0.811791
## schid2467:ses 0.62855 2.18602 0.288 0.773713
## schid2526:ses -2.34908 1.94608 -1.207 0.227442
## schid2626:ses 1.59110 2.27863 0.698 0.485033
## schid2629:ses -2.28623 1.82841 -1.250 0.211197
## schid2639:ses -3.13869 2.09015 -1.502 0.133232
## schid2651:ses 2.39048 2.10969 1.133 0.257215
## schid2655:ses 2.71846 1.96496 1.383 0.166566
## schid2658:ses 0.12132 2.01603 0.060 0.952016
## schid2755:ses -1.94808 2.01086 -0.969 0.332688
## schid2768:ses 1.00370 1.97985 0.507 0.612201
## schid2771:ses 1.75961 2.14607 0.820 0.412291
## schid2818:ses 0.29383 2.08559 0.141 0.887966
## schid2917:ses -1.37273 1.75762 -0.781 0.434822
## schid2990:ses -1.18404 1.94416 -0.609 0.542529
## schid2995:ses -1.07627 1.75075 -0.615 0.538741
## schid3013:ses 1.33122 2.25710 0.590 0.555350
## schid3020:ses -0.85490 1.87408 -0.456 0.648280
## schid3039:ses 0.44709 2.39099 0.187 0.851676
## schid3088:ses -0.71724 1.88286 -0.381 0.703267
## schid3152:ses 0.25967 1.74123 0.149 0.881456
## schid3332:ses -0.47763 2.23143 -0.214 0.830517
## schid3351:ses -0.05354 1.84308 -0.029 0.976825
## schid3377:ses -3.25544 1.94163 -1.677 0.093656 .
## schid3427:ses -2.99675 1.94053 -1.544 0.122564
## schid3498:ses -2.63967 2.03762 -1.295 0.195204
## schid3499:ses -1.51620 2.12860 -0.712 0.476304
## schid3533:ses -2.82035 2.17571 -1.296 0.194919
## schid3610:ses 0.44727 1.86805 0.239 0.810780
## schid3657:ses 1.22733 1.79993 0.682 0.495341
## schid3688:ses -0.97186 2.23346 -0.435 0.663477
## schid3705:ses -1.35008 1.94908 -0.693 0.488536
## schid3716:ses 3.35521 1.85019 1.813 0.069808 .
## schid3838:ses -1.91068 1.88447 -1.014 0.310662
## schid3881:ses -0.11788 2.28965 -0.051 0.958943
## schid3967:ses 0.80249 1.96293 0.409 0.682683
## schid3992:ses -1.97071 1.98880 -0.991 0.321768
## schid3999:ses 1.05839 1.73366 0.610 0.541553
## schid4042:ses -0.81496 1.86568 -0.437 0.662259
## schid4173:ses 0.85709 2.03937 0.420 0.674300
## schid4223:ses -0.02200 2.07409 -0.011 0.991539
## schid4253:ses -2.90813 1.83270 -1.587 0.112604
## schid4292:ses -2.66919 1.83901 -1.451 0.146707
## schid4325:ses 0.24747 1.75933 0.141 0.888142
## schid4350:ses 1.86334 2.13808 0.871 0.383512
## schid4383:ses 3.67161 2.55073 1.439 0.150072
## schid4410:ses 0.25161 2.11021 0.119 0.905091
## schid4420:ses 0.45008 2.13376 0.211 0.832945
## schid4458:ses -1.37674 1.97483 -0.697 0.485734
## schid4511:ses -2.46607 1.98367 -1.243 0.213842
## schid4523:ses -0.12779 1.89205 -0.068 0.946152
## schid4530:ses -0.86116 1.88798 -0.456 0.648314
## schid4642:ses 0.76380 1.84544 0.414 0.678970
## schid4868:ses -1.22211 2.05816 -0.594 0.552674
## schid4931:ses -1.59674 1.84706 -0.864 0.387357
## schid5192:ses -0.90509 2.04041 -0.444 0.657360
## schid5404:ses -1.29435 2.14745 -0.603 0.546704
## schid5619:ses 2.74895 1.90061 1.446 0.148123
## schid5640:ses 1.31916 1.98941 0.663 0.507294
## schid5650:ses -1.82796 1.84619 -0.990 0.322146
## schid5667:ses 1.01438 1.88893 0.537 0.591275
## schid5720:ses -0.04228 1.90514 -0.022 0.982297
## schid5761:ses 0.59943 1.85689 0.323 0.746847
## schid5762:ses -3.52268 2.42247 -1.454 0.145944
## schid5783:ses 0.60547 2.06417 0.293 0.769282
## schid5815:ses 0.50942 2.59139 0.197 0.844160
## schid5819:ses -0.53606 2.02509 -0.265 0.791242
## schid5838:ses -0.65553 2.27965 -0.288 0.773695
## schid5937:ses -1.46897 2.45967 -0.597 0.550380
## schid6074:ses -0.97949 1.93037 -0.507 0.611881
## schid6089:ses -0.81613 2.08251 -0.392 0.695148
## schid6144:ses 0.26169 1.92946 0.136 0.892119
## schid6170:ses 2.30320 2.17944 1.057 0.290646
## schid6291:ses 1.47229 2.21286 0.665 0.505861
## schid6366:ses -0.99106 1.87232 -0.529 0.596601
## schid6397:ses 0.25042 1.76272 0.142 0.887031
## schid6415:ses 1.02149 1.87230 0.546 0.585371
## schid6443:ses -3.25193 2.03987 -1.594 0.110942
## schid6464:ses -1.50509 2.31470 -0.650 0.515565
## schid6469:ses -0.75329 1.92349 -0.392 0.695346
## schid6484:ses -1.90290 2.06380 -0.922 0.356541
## schid6578:ses -0.11804 1.93329 -0.061 0.951317
## schid6600:ses 2.19571 1.84002 1.193 0.232791
## schid6808:ses -0.23248 1.83242 -0.127 0.899048
## schid6816:ses -1.15586 1.84187 -0.628 0.530320
## schid6897:ses 1.07191 1.84629 0.581 0.561548
## schid6990:ses -1.56089 1.72986 -0.902 0.366917
## schid7011:ses 2.56607 2.23770 1.147 0.251527
## schid7101:ses -1.21313 2.01720 -0.601 0.547599
## schid7172:ses -1.51410 1.97357 -0.767 0.442996
## schid7232:ses 2.49302 2.05216 1.215 0.224474
## schid7276:ses 1.26478 1.77762 0.712 0.476796
## schid7332:ses -0.04538 1.94385 -0.023 0.981376
## schid7341:ses -0.80488 1.74336 -0.462 0.644323
## schid7342:ses -1.49612 2.01193 -0.744 0.457129
## schid7345:ses 1.70334 1.73433 0.982 0.326069
## schid7364:ses -2.24909 2.31765 -0.970 0.331873
## schid7635:ses -0.06011 1.90781 -0.032 0.974867
## schid7688:ses -2.39224 2.05023 -1.167 0.243325
## schid7697:ses 0.62764 2.27534 0.276 0.782676
## schid7734:ses 3.52665 2.05743 1.714 0.086555 .
## schid7890:ses -3.16455 2.02869 -1.560 0.118830
## schid7919:ses 1.48079 2.36005 0.627 0.530391
## schid8009:ses -0.95171 2.13225 -0.446 0.655365
## schid8150:ses -2.69429 2.07623 -1.298 0.194440
## schid8165:ses -0.70634 1.98655 -0.356 0.722180
## schid8175:ses -0.89621 2.10015 -0.427 0.669584
## schid8188:ses 1.88901 2.18394 0.865 0.387093
## schid8193:ses -0.17337 2.16073 -0.080 0.936051
## schid8202:ses 1.19732 2.05121 0.584 0.559432
## schid8357:ses 0.16720 1.83514 0.091 0.927408
## schid8367:ses -2.25821 2.64174 -0.855 0.392682
## schid8477:ses 1.30357 1.89867 0.687 0.492375
## schid8531:ses 0.80965 1.99919 0.405 0.685500
## schid8627:ses -0.63902 1.85433 -0.345 0.730398
## schid8628:ses -1.27719 1.96170 -0.651 0.515029
## schid8707:ses 0.88295 1.79905 0.491 0.623592
## schid8775:ses -1.50713 1.94202 -0.776 0.437740
## schid8800:ses 0.05954 2.01597 0.030 0.976438
## schid8854:ses -0.56974 1.96540 -0.290 0.771914
## schid8857:ses -1.70636 2.02152 -0.844 0.398646
## schid8874:ses 1.58772 2.02195 0.785 0.432338
## schid8946:ses -0.81810 1.77996 -0.460 0.645806
## schid8983:ses -1.12264 1.95764 -0.573 0.566347
## schid9021:ses 0.01558 1.99352 0.008 0.993766
## schid9104:ses -1.01460 2.32244 -0.437 0.662222
## schid9158:ses 1.35263 1.81523 0.745 0.456202
## schid9198:ses 0.10197 2.21707 0.046 0.963318
## schid9225:ses 0.37731 1.99467 0.189 0.849973
## schid9292:ses -1.74991 3.09529 -0.565 0.571858
## schid9340:ses 0.80092 2.45686 0.326 0.744441
## schid9347:ses 0.17741 1.84995 0.096 0.923602
## schid9359:ses -3.34206 2.04455 -1.635 0.102175
## schid9397:ses -0.06214 2.05105 -0.030 0.975832
## schid9508:ses 1.44521 2.31249 0.625 0.532019
## schid9550:ses 1.38336 2.03924 0.678 0.497561
## schid9586:ses -0.83650 1.95373 -0.428 0.668552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.06 on 6865 degrees of freedom
## Multiple R-squared: 0.2583, Adjusted R-squared: 0.2238
## F-statistic: 7.495 on 319 and 6865 DF, p-value: < 2.2e-16fschmath <- fitted(mathbysch)
school.list <- unique(hsb$schid)
length(school.list)## [1] 160coef <- coefficients(mathbysch)
intercepts <- c(0, coef[3:161])
slopes <- c(0, coef[162:320])
# Plot all regressions in one plot
plot(data = hsb, mathach~ses,type = 'n',
ylim = c(-10, 30),
xlim = c(-3, 3),
cex.main = 1.5,
xlab = 'SES',
ylab = "Math Scores",
main = "Separate Regression for Each School "
)
for(i in 1:length(school.list)){
abline(a = (coef[1] + intercepts[i]), b = (coef[2]+slopes[i]), col = 'blue')
par<- par(new=F)
}
abline(v=c(0), lwd=2, color='black')## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
Create some variables
The data file already has mean SES, but need school mean centered ses
We might want to have a variables that indicates the proportion (mean) of students that are minority within each school. There will be 160 such proportions and we will need to merge them into our main data set.
meanminority <- aggregate(minority ~ id , data=hsb, FUN="mean")
names(meanminority) <- c('id', 'meanminority')
hsb <- merge(hsb,meanminority, by =c('id'))
head(hsb)## id minority female ses mathach size sector pracad disclim himinty
## 1 1224 0 1 -1.528 5.876 842 0 0.35 1.597 0
## 2 1224 0 1 -0.588 19.708 842 0 0.35 1.597 0
## 3 1224 0 0 -0.528 20.349 842 0 0.35 1.597 0
## 4 1224 0 0 -0.668 8.781 842 0 0.35 1.597 0
## 5 1224 0 0 -0.158 17.898 842 0 0.35 1.597 0
## 6 1224 0 0 0.022 4.583 842 0 0.35 1.597 0
## meanses cses schid meanminority
## 1 -0.428 -1.10 1224 0.08510638
## 2 -0.428 -0.16 1224 0.08510638
## 3 -0.428 -0.10 1224 0.08510638
## 4 -0.428 -0.24 1224 0.08510638
## 5 -0.428 0.27 1224 0.08510638
## 6 -0.428 0.45 1224 0.08510638The same goes for gender,
meanfemale <- aggregate(female ~ id , data=hsb, FUN="mean")
names(meanfemale) <- c('id', 'meanfemale')
hsb <- merge(hsb,meanfemale, by =c('id'))
head(hsb)## id minority female ses mathach size sector pracad disclim himinty
## 1 1224 0 1 -1.528 5.876 842 0 0.35 1.597 0
## 2 1224 0 1 -0.588 19.708 842 0 0.35 1.597 0
## 3 1224 0 0 -0.528 20.349 842 0 0.35 1.597 0
## 4 1224 0 0 -0.668 8.781 842 0 0.35 1.597 0
## 5 1224 0 0 -0.158 17.898 842 0 0.35 1.597 0
## 6 1224 0 0 0.022 4.583 842 0 0.35 1.597 0
## meanses cses schid meanminority meanfemale
## 1 -0.428 -1.10 1224 0.08510638 0.5957447
## 2 -0.428 -0.16 1224 0.08510638 0.5957447
## 3 -0.428 -0.10 1224 0.08510638 0.5957447
## 4 -0.428 -0.24 1224 0.08510638 0.5957447
## 5 -0.428 0.27 1224 0.08510638 0.5957447
## 6 -0.428 0.45 1224 0.08510638 0.5957447If you want to save data so that you don't got through steps to create is again, you can write the results to a file. There are other file formats that you can use (e.g., csv)
write.table(hsb, 'hsb.txt', row.names=F, na='.')Model Fitting
We will start simple and work to more complex
Empty/Null HLM
This model has no predictor variables. The formula that is entered is very simple to what you would enter for lm; however, what is entered here is a linear mixed model. The random part of the model is indicated by the "(1 | id)". The value "1" indicates an interecpt so (1 | id) means random intercept; that is \(U_{0j}\).
By default a fixed intercept, \(\gamma_{00}\), is included; however, I made it explicit below.
Note that I included "REML=FALSE". I want to perform maximum likelihood estimation. Late we will talk about different estimation methods.
model0 <- lmer(mathach ~ 1 + (1 | id), data=hsb, REML=FALSE)
summary(model0)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 47121.8 47142.4 -23557.9 47115.8 7182
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06262 -0.75365 0.02676 0.76070 2.74184
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.553 2.925
## Residual 39.148 6.257
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6371 0.2436 157.6209 51.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The values in the table labeled "Random effects" are the variances and standard deviations of the random effects. The row labeled "id (intercept)" contains the \(\hat{\tau}_{0}^2\) and the row labeled "residual" contains \(\hat{\sigma}^2\). The columns labeled "Standard deviation" contain \(\sqrt{\hat{\tau}_{0}^2}\)) and \(\sqrt{\hat{\sigma}^2}\)).
The values in the table labeled "Fixed effects" is the \(\hat{\gamma}_{00}\) an, it's standard error, The columns, "df", "t value" and "Pr(>|t|)" were generated by lmerTest. Note there is information about this at the top of the output.
ICC
I know that I've seen a packge that computes this, but it's easy to do so. There are (at least) three ways to compute this that range from brute force to writing your own function.
Cut-n-paste
Cut the values of \(\hat{\tau}_{0}^2\) and \(\hat{\sigma}^2\) from the summary of the model and use these
# Cut and paste
icc <- 8.553/(8.553+39.148)
icc## [1] 0.1793044Extract rather than cut and paste
The next method extracts information and puts it into formula for an ICC
vars <- as.data.frame(VarCorr(model0))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
icc## [1] 0.1793109Write a function
The last method, estentially take the last method and turns it into a little function
icc.lmer <- function(model){
vars <- as.data.frame(VarCorr(model))[4]
total <- sum(vars)
tau00 <- vars[1,1]
icc <- tau00/total
return(icc)
}
# example of how to use this
icc.lmer(model0)## [1] 0.1793109Note that the first method will contain some round off error; however, it's good enough.
Random intercept with Fixed Slope for SES
Just add "ses" to the fixed part of the model,
model1 <- lmer(mathach ~ 1 + ses + (1 | id), data=hsb, REML=FALSE)
summary(model1)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + ses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46649.0 46676.5 -23320.5 46641.0 7181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1263 -0.7277 0.0220 0.7578 2.9186
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 4.729 2.175
## Residual 37.030 6.085
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6576 0.1873 149.1760 67.57 <2e-16 ***
## ses 2.3915 0.1057 6837.3052 22.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ses 0.003icc.lmer(model1)## [1] 0.1132352Since we now have two models, it would be nice to have an easy way to compare them
screenreg(list(model0,model1), single.row=TRUE,
custom.model.names=c("Null","+SES"),
custom.header=list("Random Intercept Models"=1:2)
)##
## ===============================================================
## Random Intercept Models
## ------------------------------------------
## Null +SES
## ---------------------------------------------------------------
## (Intercept) 12.64 (0.24) *** 12.66 (0.19) ***
## ses 2.39 (0.11) ***
## ---------------------------------------------------------------
## AIC 47121.81 46649.00
## BIC 47142.45 46676.52
## Log Likelihood -23557.91 -23320.50
## Num. obs. 7185 7185
## Num. groups: id 160 160
## Var: id (Intercept) 8.55 4.73
## Var: Residual 39.15 37.03
## ===============================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Random intercept with level 1 and level 2 predictors
This is where being able to write the linear mixed model is good thing. The level 2 variable, meanses, is just added to the fixed part of the model.
model2 <- lmer(mathach ~ 1 + ses + meanses + (1 | id), data=hsb, REML=FALSE)
summary(model2)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + ses + meanses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46573.8 46608.2 -23281.9 46563.8 7180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1675 -0.7257 0.0176 0.7553 2.9445
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.647 1.627
## Residual 37.014 6.084
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6616 0.1484 155.6417 85.315 <2e-16 ***
## ses 2.1912 0.1087 7022.4974 20.165 <2e-16 ***
## meanses 3.6745 0.3754 184.9826 9.788 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ses
## ses 0.004
## meanses -0.005 -0.289Use School centered SES
as a level 1 predictor,
hsb$cses <- hsb$ses - hsb$meanses
model3 <- lmer(mathach ~ 1 + cses + (1 | id), data=hsb, REML=FALSE)
summary(model3)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46728.4 46755.9 -23360.2 46720.4 7181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09692 -0.73196 0.01945 0.75739 2.91423
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.612 2.935
## Residual 37.005 6.083
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6494 0.2437 157.7181 51.90 <2e-16 ***
## cses 2.1912 0.1086 7023.0025 20.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cses 0.003icc.lmer(model3)## [1] 0.1887876Some Graphics:
We will plot both the population average models as well as school specific ("conditional") regression line(s).
We need to get fixed coefficients from the output,
intercept <- fixef(model3)[1]
slope <- fixef(model3)[2]For the conditional/school specific we will need estimates of the random effects, \(\hat{U}_{0j}\).
intercept.j <- ranef(model3)
student.int <- intercept.j$id[1]We will start by setting up a graphics window (i.e., type='n'). In the first graph, we only plot the population average regression line. In the 2nd figure, we graph both the popuation average line and the school specific regressions.
Popuation Average Regression Line
plot(mathach~cses,type = 'n', data = hsb,
ylim = c(-10, 30),
xlim = c(-3, 3),
cex.main = 1.5,
xlab = 'School Centered SES',
ylab = "Fitted Math Achievement Scores",
main = "Population Average Regression "
)
# Population average line
for(i in min(hsb$cses):max(hsb$cses)){
abline(a=intercept, b=slope, color='black',lwd=4)
}## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
## Warning in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): "color" is
## not a graphical parameter
Population and School Specific
plot(mathach~cses,type = 'n', data = hsb,
ylim = c(-10, 30),
xlim = c(-3, 3),
cex.main = 1.5,
xlab = 'School Centered SES',
ylab = "Fitted Math Achievement Scores",
main = "Population & School Specific Regressions "
)
# draw in school specific line
for (i in 1:length(school.list)){
abline(a=(intercept+student.int[i,]), b=slope, type='l', col=i)
}
# Population average line
for(i in min(hsb$cses):max(hsb$cses)){
abline(a=intercept, b=slope, col='black',lwd=4)
}
The black line is the population average.
Model 4: center ses & add back in mean ses
model4 <- lmer(mathach ~ 1 + cses + meanses + (1 | id), data=hsb, REML=FALSE)
summary(model4)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46573.8 46608.2 -23281.9 46563.8 7180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1675 -0.7257 0.0176 0.7553 2.9445
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.647 1.627
## Residual 37.014 6.084
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6616 0.1484 155.6417 85.31 <2e-16 ***
## cses 2.1912 0.1087 7022.4974 20.16 <2e-16 ***
## meanses 5.8656 0.3594 155.3036 16.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses
## cses 0.004
## meanses -0.004 0.000icc.lmer(model4)## [1] 0.06673899Model 5: Lots of level 1 variables
model5 <- lmer(mathach ~ 1 + cses + female + minority + meanses + (1 | id), data=hsb, REML=FALSE)
summary(model5)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46347.3 46395.4 -23166.6 46333.3 7178
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2173 -0.7222 0.0364 0.7629 2.9090
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.396 1.548
## Residual 35.886 5.990
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.0483 0.1749 316.2213 80.317 < 2e-16 ***
## cses 1.9265 0.1084 7070.5416 17.766 < 2e-16 ***
## female -1.2185 0.1607 6014.8042 -7.582 3.93e-14 ***
## minority -2.7282 0.2032 4182.5246 -13.427 < 2e-16 ***
## meanses 4.8085 0.3524 164.8798 13.647 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty
## cses -0.070
## female -0.485 0.051
## minority -0.327 0.155 0.016
## meanses -0.088 0.033 0.042 0.202icc.lmer(model5)## [1] 0.06259292Model 6: Lots of level 1 & 2 variables
model6 <- lmer(mathach ~ 1 + cses + female + minority + meanses + himinty + pracad
+ disclim + sector + size +(1 | id), data=hsb, REML=FALSE)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(model6)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + himinty +
## pracad + disclim + sector + size + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46289.5 46372.1 -23132.8 46265.5 7173
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2734 -0.7212 0.0386 0.7570 2.9654
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.292 1.137
## Residual 35.880 5.990
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.131e+01 5.016e-01 1.768e+02 22.547 < 2e-16 ***
## cses 1.904e+00 1.085e-01 7.053e+03 17.540 < 2e-16 ***
## female -1.261e+00 1.580e-01 4.997e+03 -7.979 1.82e-15 ***
## minority -2.989e+00 2.111e-01 6.363e+03 -14.158 < 2e-16 ***
## meanses 3.082e+00 4.219e-01 1.527e+02 7.305 1.44e-11 ***
## himinty 2.437e-01 3.182e-01 1.898e+02 0.766 0.444721
## pracad 3.003e+00 7.992e-01 1.586e+02 3.757 0.000241 ***
## disclim -3.868e-01 1.795e-01 1.678e+02 -2.155 0.032583 *
## sector 8.391e-01 3.863e-01 1.516e+02 2.172 0.031391 *
## size 7.430e-04 2.139e-04 1.604e+02 3.474 0.000659 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss himnty pracad disclm sector
## cses -0.011
## female -0.227 0.051
## minority -0.009 0.161 0.019
## meanses 0.472 0.015 0.005 0.092
## himinty 0.139 -0.055 -0.033 -0.331 0.413
## pracad -0.726 -0.002 0.056 -0.027 -0.587 -0.171
## disclim -0.322 0.000 0.091 -0.026 -0.032 -0.052 0.221
## sector -0.161 -0.006 0.027 -0.046 0.008 -0.132 -0.342 0.487
## size -0.608 -0.004 0.024 -0.030 -0.158 -0.183 0.097 -0.031 0.274
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingicc.lmer(model6)## [1] 0.03476522Model 7: Kitchen Sink
We add even more variables. Note this is not how we should go about modeling the data; however, this is for showing you how to fit complex models.
We will ignore the warning messages for now...
model7 <- lmer(mathach ~ 1 + cses + female + minority + meanses + meanminority + meanfemale + himinty + pracad + disclim + sector + size +(1 | id), data=hsb, REML=FALSE)## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(model7)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + meanminority +
## meanfemale + himinty + pracad + disclim + sector + size + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46289.4 46385.7 -23130.7 46261.4 7171
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2596 -0.7209 0.0369 0.7583 3.0062
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.249 1.118
## Residual 35.876 5.990
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.187e+01 5.874e-01 1.590e+02 20.204 < 2e-16 ***
## cses 1.912e+00 1.087e-01 7.024e+03 17.599 < 2e-16 ***
## female -1.163e+00 1.679e-01 7.024e+03 -6.928 4.66e-12 ***
## minority -2.924e+00 2.194e-01 7.024e+03 -13.327 < 2e-16 ***
## meanses 2.927e+00 4.391e-01 1.502e+02 6.665 4.71e-10 ***
## meanminority -8.695e-01 7.987e-01 1.818e+02 -1.089 0.277719
## meanfemale -8.599e-01 4.926e-01 1.860e+02 -1.745 0.082563 .
## himinty 6.837e-01 4.851e-01 1.610e+02 1.409 0.160622
## pracad 2.867e+00 8.044e-01 1.587e+02 3.564 0.000482 ***
## disclim -4.470e-01 1.841e-01 1.676e+02 -2.427 0.016269 *
## sector 8.569e-01 3.884e-01 1.508e+02 2.206 0.028888 *
## size 7.414e-04 2.134e-04 1.613e+02 3.474 0.000658 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss mnmnrt menfml himnty pracad
## cses 0.001
## female 0.000 0.054
## minority 0.000 0.167 0.019
## meanses 0.368 0.000 0.000 0.000
## meanminorty -0.023 -0.046 -0.005 -0.275 0.310
## meanfemale -0.534 -0.018 -0.341 -0.006 0.013 0.018
## himinty 0.120 0.000 0.000 0.000 0.019 -0.759 -0.061
## pracad -0.684 0.000 0.000 0.000 -0.577 -0.095 0.153 -0.045
## disclim -0.392 0.000 0.000 0.000 -0.055 -0.090 0.245 0.024 0.257
## sector -0.171 0.000 0.000 0.000 -0.042 -0.160 0.073 0.034 -0.304
## size -0.544 0.000 0.000 0.000 -0.182 -0.106 0.067 -0.040 0.116
## disclm sector
## cses
## female
## minority
## meanses
## meanminorty
## meanfemale
## himinty
## pracad
## disclim
## sector 0.496
## size -0.003 0.289
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingicc.lmer(model7)## [1] 0.03364665Quick comparison
screenreg(list(model0,model5,model6,model7),
single.row=TRUE,
custom.model.names=c("Empty/Null", "Add Level 1", "Level 1 & 2","Kitchen Sink"),
custom.head=list("Simple"=1, "Complex"=2:4)
)##
## ===========================================================================================================
## Simple Complex
## -------------------- ----------------------------------------------------------------
## Empty/Null Add Level 1 Level 1 & 2 Kitchen Sink
## -----------------------------------------------------------------------------------------------------------
## (Intercept) 12.64 (0.24) *** 14.05 (0.17) *** 11.31 (0.50) *** 11.87 (0.59) ***
## cses 1.93 (0.11) *** 1.90 (0.11) *** 1.91 (0.11) ***
## female -1.22 (0.16) *** -1.26 (0.16) *** -1.16 (0.17) ***
## minority -2.73 (0.20) *** -2.99 (0.21) *** -2.92 (0.22) ***
## meanses 4.81 (0.35) *** 3.08 (0.42) *** 2.93 (0.44) ***
## himinty 0.24 (0.32) 0.68 (0.49)
## pracad 3.00 (0.80) *** 2.87 (0.80) ***
## disclim -0.39 (0.18) * -0.45 (0.18) *
## sector 0.84 (0.39) * 0.86 (0.39) *
## size 0.00 (0.00) *** 0.00 (0.00) ***
## meanminority -0.87 (0.80)
## meanfemale -0.86 (0.49)
## -----------------------------------------------------------------------------------------------------------
## AIC 47121.81 46347.27 46289.51 46289.38
## BIC 47142.45 46395.43 46372.07 46385.70
## Log Likelihood -23557.91 -23166.63 -23132.76 -23130.69
## Num. obs. 7185 7185 7185 7185
## Num. groups: id 160 160 160 160
## Var: id (Intercept) 8.55 2.40 1.29 1.25
## Var: Residual 39.15 35.89 35.88 35.88
## ===========================================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05