LMM_Examples
This example is from West, Welch, and Galecki (2014)’s “Linear Mixed Models: A Practical Guide Using Statistical Software (Second Edition).
setwd("C:/Liem/GEOG515/Spring15/Labs/")
class <- read.table("classroom.txt", header=TRUE, sep="\t")
str(class)## 'data.frame': 1190 obs. of 12 variables:
## $ sex : int 1 0 1 0 0 1 0 0 1 0 ...
## $ minority: int 1 1 1 1 1 1 1 1 1 1 ...
## $ mathkind: int 448 460 511 449 425 450 452 443 422 480 ...
## $ mathgain: int 32 109 56 83 53 65 51 66 88 -7 ...
## $ ses : num 0.46 -0.27 -0.03 -0.38 -0.03 0.76 -0.03 0.2 0.64 0.13 ...
## $ yearstea: num 1 1 1 2 2 2 2 2 2 2 ...
## $ mathknow: num NA NA NA -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 -0.11 ...
## $ housepov: num 0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 0.082 ...
## $ mathprep: num 2 2 2 3.25 3.25 3.25 3.25 3.25 3.25 3.25 ...
## $ classid : int 160 160 160 217 217 217 217 217 217 217 ...
## $ schoolid: int 1 1 1 1 1 1 1 1 1 1 ...
## $ childid : int 1 2 3 4 5 6 7 8 9 10 ...library(nlme)
#classData: 1,190 first-grade students sampled from 312 classrooms in 107 schools. Dependent variable: MATHGAIN, change in student math achievement scores from the spring of kindergarten to the spring of first grade. Three-level data set: Level 1: students Level 2: classrooms Level 3: schools Question: the contribution of selected student-level, classroom-level and school-level covariates to the variation in students’ math achievement gain.
Descriptive Summary
# Level 1 Descriptive Statistics
attach(class)
level1 <- class[1:5]
dim(level1)## [1] 1190 5summary(level1)## sex minority mathkind mathgain
## Min. :0.0000 Min. :0.0000 Min. :290.0 Min. :-110.00
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:439.2 1st Qu.: 35.00
## Median :1.0000 Median :1.0000 Median :466.0 Median : 56.00
## Mean :0.5059 Mean :0.6773 Mean :466.7 Mean : 57.57
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:495.0 3rd Qu.: 77.00
## Max. :1.0000 Max. :1.0000 Max. :629.0 Max. : 253.00
## ses
## Min. :-1.61000
## 1st Qu.:-0.49000
## Median :-0.03000
## Mean :-0.01298
## 3rd Qu.: 0.39750
## Max. : 3.21000# Level 2 Descriptive Statistics
level2 <- aggregate(class,list(classid = class$classid),mean)
summary(level2$yearstea)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 4.00 10.00 12.28 20.00 40.00dim(level2)## [1] 312 13# Level 3 Descriptive Statistics
level3 <- aggregate(class,list(schoolid = class$schoolid),mean)
summary(level3$housepov)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0120 0.0855 0.1480 0.1941 0.2645 0.5640dim(level3)## [1] 107 13Descriptive Summary #2
class.nomiss <- subset(class, !is.na(mathknow))
# Level 1 Descriptive Statistics
level1 <- data.frame(class.nomiss$sex,class.nomiss$minority,class.nomiss$mathkind,class.nomiss$mathgain,class.nomiss$ses)
summary(level1)## class.nomiss.sex class.nomiss.minority class.nomiss.mathkind
## Min. :0.0000 Min. :0.0000 Min. :290.0
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:440.0
## Median :1.0000 Median :1.0000 Median :467.0
## Mean :0.5023 Mean :0.6661 Mean :467.2
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:495.0
## Max. :1.0000 Max. :1.0000 Max. :629.0
## class.nomiss.mathgain class.nomiss.ses
## Min. :-84.00 Min. :-1.610000
## 1st Qu.: 35.00 1st Qu.:-0.490000
## Median : 56.00 Median :-0.030000
## Mean : 57.84 Mean :-0.009306
## 3rd Qu.: 77.00 3rd Qu.: 0.400000
## Max. :253.00 Max. : 3.210000dim(level1)## [1] 1081 5# Level 2 Descriptive Statistics
level2.agg <- aggregate(class.nomiss,list(classid = class.nomiss$classid),mean)
level2 <- data.frame(level2.agg$yearstea,level2.agg$mathknow,level2.agg$mathprep)
summary(level2)## level2.agg.yearstea level2.agg.mathknow level2.agg.mathprep
## Min. : 0.00 Min. :-2.50000 Min. :1.000
## 1st Qu.: 4.00 1st Qu.:-0.76000 1st Qu.:2.000
## Median :10.00 Median :-0.19000 Median :2.330
## Mean :12.28 Mean :-0.08025 Mean :2.585
## 3rd Qu.:20.00 3rd Qu.: 0.62000 3rd Qu.:3.000
## Max. :40.00 Max. : 2.61000 Max. :6.000dim(level2)## [1] 285 3# Level 3 Descriptive Statistics
level3.agg <- aggregate(class.nomiss,list(schoolid = class.nomiss$schoolid),mean)
summary(level3.agg$housepov)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0120 0.0850 0.1480 0.1954 0.2660 0.5640dim(level3.agg)## [1] 105 13Boxplot
class.first6 <- class[class$schoolid <= 6,]
par(mfrow=c(3,2))
for (i in 1:6)
{boxplot(class.first6$mathgain[class.first6$schoolid==i] ~
class.first6$classid[class.first6$schoolid==i])}
Step 1: Fit the random effects model (Model 4.1). Fit a three-level model with a fixed intercept, and random effects associated with the intercept for classrooms (Level 2) and schools (Level 3)
# Model 1.
model1.fit <- lme(mathgain ~ 1, random = ~1 | schoolid/classid, class, method = "REML")
summary(model1.fit)## Linear mixed-effects model fit by REML
## Data: class
## AIC BIC logLik
## 11776.76 11797.09 -5884.382
##
## Random effects:
## Formula: ~1 | schoolid
## (Intercept)
## StdDev: 8.802955
##
## Formula: ~1 | classid %in% schoolid
## (Intercept) Residual
## StdDev: 9.961301 32.06609
##
## Fixed effects: mathgain ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 57.42703 1.443288 878 39.78903 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.64405126 -0.59835943 -0.03364297 0.53341988 5.63351835
##
## Number of Observations: 1190
## Number of Groups:
## schoolid classid %in% schoolid
## 107 312random.effects(model1.fit)## Level: schoolid
## (Intercept)
## 1 0.94154109
## 2 2.54801842
## 3 13.67745892
## 4 -6.46372188
## 5 0.60989716
## 6 8.73455373
## 7 0.28100790
## 8 -3.82861616
## 9 -6.14219882
## 10 8.28958950
## 11 -4.25424964
## 12 -2.87495377
## 13 -3.19716802
## 14 -5.66477744
## 15 7.42631397
## 16 8.61692569
## 17 6.27077279
## 18 5.77431854
## 19 -0.54546698
## 20 1.76742229
## 21 -0.51596550
## 22 -3.20256313
## 23 3.39536399
## 24 -5.44437195
## 25 -7.95706094
## 26 -0.32978307
## 27 -0.05904167
## 28 -3.27506003
## 29 -3.88074723
## 30 -4.43955611
## 31 -4.52723091
## 32 -1.30685665
## 33 1.02255801
## 34 -1.76219657
## 35 -9.56164373
## 36 3.89665085
## 37 -2.81971486
## 38 -3.40953392
## 39 -11.97371445
## 40 -7.17939444
## 41 3.76628837
## 42 0.49235426
## 43 -0.38871182
## 44 5.58679610
## 45 1.85903460
## 46 -1.23451548
## 47 5.25613514
## 48 -1.87602007
## 49 7.29736993
## 50 -5.22710212
## 51 -0.16560814
## 52 -5.21010573
## 53 4.49336766
## 54 -6.64626729
## 55 -1.84254086
## 56 -2.73987034
## 57 4.11334513
## 58 -1.50519347
## 59 -2.76046255
## 60 6.31255622
## 61 -5.55703847
## 62 -1.12136456
## 63 -0.36205614
## 64 3.16128551
## 65 -0.75350027
## 66 0.67818549
## 67 8.36809077
## 68 10.45885216
## 69 1.82148268
## 70 11.05849097
## 71 -2.70316611
## 72 -2.45634190
## 73 -2.12297284
## 74 0.15103997
## 75 7.02422044
## 76 -7.63101220
## 77 -1.05746143
## 78 -9.50512896
## 79 1.83437487
## 80 4.67552552
## 81 3.08625622
## 82 -0.48622634
## 83 -0.48933015
## 84 9.99237239
## 85 -1.42241419
## 86 -6.90004269
## 87 -4.26527037
## 88 -5.03141935
## 89 -1.37145791
## 90 -0.78729502
## 91 -4.06475433
## 92 1.53452208
## 93 -7.83154617
## 94 11.62601917
## 95 2.16314934
## 96 -3.84140092
## 97 1.87140739
## 98 -0.94527041
## 99 -2.63475660
## 100 0.83441688
## 101 2.38097612
## 102 -1.19677800
## 103 12.48049756
## 104 7.85093027
## 105 1.34616791
## 106 -5.37603618
## 107 -2.70187672
##
## Level: classid %in% schoolid
## (Intercept)
## 1/160 1.638501853
## 1/217 -0.432870842
## 2/197 -1.694663767
## 2/211 2.513497117
## 2/307 2.443870810
## 3/11 3.869973397
## 3/137 5.564542344
## 3/145 8.102749126
## 3/228 -0.023458144
## 4/48 3.770040373
## 4/179 -12.046751442
## 5/299 0.780965307
## 6/254 5.717086533
## 6/286 0.429675529
## 6/291 5.037719410
## 7/89 -0.607982947
## 7/107 -1.589717469
## 7/162 4.415322750
## 7/280 -1.857795403
## 8/147 7.746709115
## 8/265 0.320723689
## 8/279 -12.969925565
## 9/17 -2.383891616
## 9/58 -0.905164328
## 9/283 -4.575949075
## 10/178 3.443248126
## 10/208 -0.199581353
## 10/278 5.133704453
## 10/303 2.237339387
## 11/86 1.766399577
## 11/113 -0.712322637
## 11/131 -0.594185716
## 11/172 -11.327581017
## 11/195 1.622912056
## 11/201 -5.263757676
## 11/233 0.072803622
## 11/285 4.301751377
## 11/293 4.686469557
## 12/20 -1.246500201
## 12/49 0.103136694
## 12/81 3.027025581
## 12/251 0.474749087
## 12/302 -6.039751750
## 13/64 -0.270094300
## 13/166 -3.823837408
## 14/41 -5.260331191
## 14/141 -4.153732364
## 14/167 2.160390354
## 15/5 10.310039084
## 15/85 6.509854440
## 15/92 2.565187934
## 15/119 -2.275344035
## 15/260 -7.600439805
## 16/206 7.341322308
## 16/256 3.692537997
## 17/7 -1.117529168
## 17/10 6.843694309
## 17/87 -0.853500152
## 17/177 4.483225886
## 17/216 -1.326248584
## 18/31 7.393939137
## 19/98 -1.068050033
## 19/110 -1.429166200
## 19/170 1.798752917
## 20/30 -1.446681191
## 20/94 1.395475213
## 20/101 2.314367029
## 21/6 -5.699329342
## 21/130 2.113559130
## 21/158 2.925083162
## 22/159 -4.100840070
## 23/43 -1.756211006
## 23/47 4.944172980
## 23/163 1.159757077
## 24/55 -8.912239296
## 24/77 -5.983131156
## 24/93 4.754048848
## 24/150 3.169874749
## 25/97 -2.745418832
## 25/297 -7.443493313
## 26/78 4.477742943
## 26/102 -0.695364569
## 26/139 -4.204661265
## 27/104 0.598451171
## 27/190 2.181422210
## 27/230 -7.097008561
## 27/262 4.241533090
## 28/114 -6.378317430
## 28/118 1.482789891
## 28/202 -10.436450808
## 28/203 11.138306945
## 29/115 0.265616159
## 29/223 -5.234862081
## 30/175 -5.684793354
## 31/51 -2.713465536
## 31/156 -0.635529173
## 31/234 9.428509155
## 31/258 -11.876574238
## 32/191 3.861849450
## 32/204 -3.902186371
## 32/252 -1.633075869
## 33/26 7.685943165
## 33/32 -7.784105438
## 33/88 -5.977742511
## 33/100 -1.711751572
## 33/229 9.097028469
## 34/2 2.245520094
## 34/188 -4.501989685
## 35/56 -5.814650532
## 35/122 -6.428908837
## 36/63 -0.942160650
## 36/180 5.931770944
## 37/103 -7.769714459
## 37/125 5.818722871
## 37/210 -0.286021900
## 37/249 -1.373594397
## 38/16 0.726974802
## 38/309 -5.092838254
## 39/79 -8.870142159
## 39/198 -5.654783462
## 39/274 0.224133156
## 39/288 -1.031391814
## 40/9 -13.301731936
## 40/257 4.108611537
## 41/18 0.335025095
## 41/50 -1.307112774
## 41/142 5.794770524
## 42/54 -1.091568080
## 42/70 -0.696982120
## 42/207 -2.048373786
## 42/296 4.467377159
## 43/220 -1.543118569
## 43/227 1.045378189
## 44/176 3.407839319
## 44/219 -1.117870626
## 44/246 5.844760058
## 44/268 -0.980909633
## 45/300 2.380469415
## 46/19 3.770709950
## 46/69 -4.084414820
## 46/71 3.608928903
## 46/133 -1.315448573
## 46/240 -3.560556341
## 47/61 -1.890114156
## 47/82 9.671946549
## 47/155 -1.051420298
## 48/153 -2.402219091
## 49/105 10.799208237
## 49/212 -1.455021767
## 50/25 -4.233258445
## 50/127 0.404136336
## 50/289 -2.864113586
## 51/59 -0.212059045
## 52/255 -0.935164607
## 52/270 -5.736307435
## 53/14 8.238074799
## 53/149 -2.484376454
## 54/259 -6.766335533
## 54/301 -1.744122545
## 55/24 -3.131778931
## 55/35 1.299370086
## 55/117 -0.526940551
## 56/3 -8.760520447
## 56/168 5.252152421
## 57/36 2.416712875
## 57/40 7.212699413
## 57/171 0.656519432
## 57/250 -2.595028151
## 57/272 -2.423819094
## 58/12 -1.927380495
## 59/34 -4.283128354
## 59/108 0.748392276
## 60/121 8.957280988
## 60/242 -0.874135558
## 61/1 -1.324660053
## 61/124 0.011442295
## 61/192 -4.185273788
## 61/267 -1.617223315
## 62/22 -2.638289653
## 62/46 1.202397047
## 63/310 -0.463608140
## 64/38 3.506534722
## 64/294 0.541449903
## 65/95 -5.771300467
## 65/128 4.059624091
## 65/284 0.746829169
## 66/91 1.711355531
## 66/218 -1.392160238
## 66/269 0.549212331
## 67/4 9.425842968
## 67/73 3.488496917
## 67/245 -2.199109428
## 68/57 -4.633179219
## 68/244 12.859455901
## 68/263 -2.144649412
## 68/271 0.699153761
## 68/298 6.611641954
## 69/90 4.408774698
## 69/290 -2.076390021
## 70/83 6.600554007
## 70/146 4.830233145
## 70/152 4.711225306
## 70/261 -0.042730249
## 70/311 -1.939029598
## 71/80 -1.609982860
## 71/120 4.967170818
## 71/209 3.280295073
## 71/264 -6.972947927
## 71/304 -3.125903841
## 72/29 1.503320095
## 72/148 -3.960535048
## 72/243 -0.688098626
## 73/213 -1.789512653
## 73/276 -0.928926214
## 74/232 0.068260818
## 74/295 0.125143888
## 75/39 -0.234784280
## 75/42 8.028897822
## 75/154 0.007887291
## 75/181 1.192421299
## 76/8 -7.532456412
## 76/28 2.800341259
## 76/136 -1.636441312
## 76/277 -6.138101761
## 76/287 2.735247211
## 77/27 -1.602296209
## 77/53 2.158009211
## 77/189 -0.275073889
## 77/305 -1.634704604
## 78/13 -4.198711517
## 78/111 -1.179856005
## 78/183 -6.792625424
## 79/23 -0.815092914
## 79/112 -0.123180340
## 79/253 3.287166213
## 80/222 7.282200896
## 80/308 -1.295251808
## 81/116 4.446829881
## 81/174 -0.494919279
## 82/164 -0.241108467
## 82/169 3.566423171
## 82/184 -3.886375191
## 82/225 -0.061545964
## 83/165 3.628677095
## 83/235 -4.255257929
## 84/76 2.111236986
## 84/143 7.180101650
## 84/161 3.503763123
## 85/66 -6.479963566
## 85/84 0.791681106
## 85/157 1.808999375
## 85/186 -6.846658218
## 85/196 8.904558591
## 86/65 -2.686669853
## 86/132 -5.656246079
## 86/275 -0.492498213
## 87/21 -1.159070888
## 87/72 -1.682728398
## 87/106 -4.717819479
## 87/187 2.097996012
## 88/67 -5.079214237
## 88/312 -1.363452233
## 89/292 -1.756133864
## 90/52 -7.464534113
## 90/238 6.456413166
## 91/15 -1.704067242
## 91/45 -9.426357826
## 91/151 3.905026860
## 91/231 2.020533613
## 92/173 -0.802684555
## 92/199 -2.196857723
## 92/248 4.964477670
## 93/123 -2.111045964
## 93/126 -8.713452296
## 93/182 2.995797198
## 93/194 -2.199491076
## 94/134 8.410009581
## 94/224 -4.183283959
## 94/282 3.712382276
## 94/306 6.947857137
## 95/140 -2.036851863
## 95/144 4.806736213
## 96/109 -1.277120746
## 96/185 -3.659937290
## 96/237 0.018194556
## 97/75 2.396312613
## 98/221 -1.210406365
## 99/205 0.303142848
## 99/214 -2.710017595
## 99/226 0.599840430
## 99/236 4.159023472
## 99/266 -5.725760421
## 100/74 -4.247473826
## 100/99 -0.023009562
## 100/135 5.481382961
## 100/281 -0.142439705
## 101/68 -6.282521700
## 101/138 5.986047800
## 101/215 -4.984150961
## 101/273 8.329433551
## 102/44 -6.990513576
## 102/241 5.458055042
## 103/37 0.096148434
## 103/193 1.794552347
## 103/247 14.090412615
## 104/60 0.087269520
## 104/129 9.965743698
## 105/33 0.968308482
## 105/200 0.755441862
## 106/62 -6.883943809
## 107/96 -3.746854794
## 107/239 0.287137098intervals(model1.fit)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 54.59433 57.42703 60.25972
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: schoolid
## lower est. upper
## sd((Intercept)) 5.827003 8.802955 13.29878
## Level: classid
## lower est. upper
## sd((Intercept)) 6.591329 9.961301 15.05425
##
## Within-group standard error:
## lower est. upper
## 30.60181 32.06609 33.60042Decide whether to keep the random intercepts for classrooms (Model 1 vs. Model 1A).
# Model 1A.
model1A.fit <- lme(mathgain ~ 1, random = ~1 | schoolid, class, method = "REML")
anova(model1.fit, model1A.fit)## Model df AIC BIC logLik Test L.Ratio p-value
## model1.fit 1 4 11776.76 11797.09 -5884.382
## model1A.fit 2 3 11782.67 11797.91 -5888.335 1 vs 2 7.904762 0.0049Step 2: Build the Level 1 Model by adding Level 1 Covariates (Model 1 vs. Model 2). Add fixed effects associated with covariates measured on the students to the Level 1 Model to obtain Model 2, and evaluate the reduction in the residual variance.
# Model 2.
model2.fit <- lme(mathgain ~ mathkind + sex + minority + ses, random = ~1 | schoolid/classid,
class, na.action = "na.omit", method = "REML")
summary(model2.fit)## Linear mixed-effects model fit by REML
## Data: class
## AIC BIC logLik
## 11401.8 11442.42 -5692.902
##
## Random effects:
## Formula: ~1 | schoolid
## (Intercept)
## StdDev: 8.671991
##
## Formula: ~1 | classid %in% schoolid
## (Intercept) Residual
## StdDev: 9.12604 27.10286
##
## Fixed effects: mathgain ~ mathkind + sex + minority + ses
## Value Std.Error DF t-value p-value
## (Intercept) 282.79034 10.853234 874 26.055860 0.0000
## mathkind -0.46980 0.022266 874 -21.099524 0.0000
## sex -1.25119 1.657730 874 -0.754762 0.4506
## minority -8.26213 2.340113 874 -3.530655 0.0004
## ses 5.34638 1.241094 874 4.307794 0.0000
## Correlation:
## (Intr) mthknd sex minrty
## mathkind -0.978
## sex -0.044 -0.031
## minority -0.307 0.163 -0.018
## ses 0.140 -0.168 0.019 0.163
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.82574352 -0.61102781 -0.03366662 0.55380954 4.26781224
##
## Number of Observations: 1190
## Number of Groups:
## schoolid classid %in% schoolid
## 107 312# remove sex because it's not significant
model2.fit <- lme(mathgain ~ mathkind + minority + ses, random = ~1 | schoolid/classid,
class, na.action = "na.omit", method = "REML")
# Model 1: ML estimation.
model1.ml.fit <- lme(mathgain ~ 1, random = ~1 | schoolid/classid, class,
na.action = "na.omit",method = "ML")
# Model 2: ML estimation.
model2.ml.fit <- lme(mathgain ~ mathkind + minority + ses,
random = ~1 | schoolid/classid, class,
na.action = "na.omit", method = "ML")
anova(model2.ml.fit)## numDF denDF F-value p-value
## (Intercept) 1 875 1924.5228 <.0001
## mathkind 1 875 409.9155 <.0001
## minority 1 875 18.6378 <.0001
## ses 1 875 18.6713 <.0001anova(model1.ml.fit, model2.ml.fit)## Model df AIC BIC logLik Test L.Ratio p-value
## model1.ml.fit 1 4 11779.33 11799.66 -5885.666
## model2.ml.fit 2 7 11405.53 11441.10 -5695.766 1 vs 2 379.7994 <.0001Step 3: Build the Level 2 Model by adding Level 2 covariates (Model 3). Add fixed effects associated with the covariates measured on the Level 2 clusters (classrooms) to the Level 2 model to create Model 3 and decide whether to retain the effects of the Level 2 covariates in the model.
# Model 3.
model3.fit <- update(model2.fit, fixed = ~ mathkind + sex +
minority + ses + yearstea + mathprep + mathknow)
summary(model3.fit)## Linear mixed-effects model fit by REML
## Data: class
## AIC BIC logLik
## 10335 10389.76 -5156.499
##
## Random effects:
## Formula: ~1 | schoolid
## (Intercept)
## StdDev: 8.671285
##
## Formula: ~1 | classid %in% schoolid
## (Intercept) Residual
## StdDev: 9.310153 26.71766
##
## Fixed effects: mathgain ~ mathkind + sex + minority + ses + yearstea + mathprep + mathknow
## Value Std.Error DF t-value p-value
## (Intercept) 282.02452 11.701687 792 24.101185 0.0000
## mathkind -0.47501 0.022747 792 -20.882471 0.0000
## sex -1.33950 1.718580 792 -0.779423 0.4360
## minority -7.86886 2.418081 792 -3.254177 0.0012
## ses 5.41925 1.275995 792 4.247078 0.0000
## yearstea 0.03974 0.117070 177 0.339435 0.7347
## mathprep 1.09485 1.148493 177 0.953296 0.3417
## mathknow 1.91448 1.147015 177 1.669094 0.0969
## Correlation:
## (Intr) mthknd sex minrty ses yearst mthprp
## mathkind -0.943
## sex -0.069 -0.005
## minority -0.302 0.166 -0.015
## ses 0.121 -0.163 0.020 0.159
## yearstea -0.092 0.005 0.015 0.040 -0.033
## mathprep -0.283 0.052 -0.002 0.017 0.043 -0.168
## mathknow -0.049 0.029 0.007 0.146 -0.015 0.015 0.004
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.3989531 -0.6289628 -0.0413954 0.5564112 4.3029072
##
## Number of Observations: 1081
## Number of Groups:
## schoolid classid %in% schoolid
## 105 285Step 4: Build the Level 3 Model by adding the Level 3 covariate (Model 4). Add a fixed effect associated with the covariate measured on the Level 3 clusters (schools) to the Level 3 model to create Model 4, and evaluate the reduction in the variance component associated with the Level 3 clusters.
# Model 4.
model4.fit <- update(model2.fit, fixed = ~ mathkind + sex + minority + ses + housepov)
anova(model4.fit)## numDF denDF F-value p-value
## (Intercept) 1 874 1875.7391 <.0001
## mathkind 1 874 410.3731 <.0001
## sex 1 874 0.8613 0.3536
## minority 1 874 18.3205 <.0001
## ses 1 874 18.5953 <.0001
## housepov 1 105 1.3251 0.2523summary(model4.fit)## Linear mixed-effects model fit by REML
## Data: class
## AIC BIC logLik
## 11396.06 11441.75 -5689.029
##
## Random effects:
## Formula: ~1 | schoolid
## (Intercept)
## StdDev: 8.818243
##
## Formula: ~1 | classid %in% schoolid
## (Intercept) Residual
## StdDev: 9.030822 27.10018
##
## Fixed effects: mathgain ~ mathkind + sex + minority + ses + housepov
## Value Std.Error DF t-value p-value
## (Intercept) 285.05800 11.020766 874 25.865534 0.0000
## mathkind -0.47086 0.022281 874 -21.132931 0.0000
## sex -1.23460 1.657434 874 -0.744884 0.4565
## minority -7.75588 2.384993 874 -3.251950 0.0012
## ses 5.23971 1.244971 874 4.208703 0.0000
## housepov -11.43923 9.937384 105 -1.151131 0.2523
## Correlation:
## (Intr) mthknd sex minrty ses
## mathkind -0.969
## sex -0.042 -0.032
## minority -0.265 0.153 -0.015
## ses 0.123 -0.165 0.019 0.144
## housepov -0.173 0.035 -0.009 -0.184 0.078
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.79254947 -0.60523530 -0.03270737 0.55242345 4.27266338
##
## Number of Observations: 1190
## Number of Groups:
## schoolid classid %in% schoolid
## 107 312Diagnostics for the final model
library(lattice)
trellis.device(color=F)
qqnorm(model2.fit, ~ranef(.,level=2))
qqnorm(model2.fit, ~ranef(.,level=1))
qqnorm(model2.fit, ~resid(.), plot.it = TRUE)
plot(resid(model2.fit, type="p") ~ fitted(model2.fit))
abline(h = 0, lty = 2)
lines(lowess(resid(model2.fit, type="p") ~ fitted(model2.fit)))