Multilevel Models
Julian Faraway
[1] "Created: Wed Apr 1 16:50:39 2015"See the introduction for an overview.
library(lme4)
library(ggplot2)
options(digits=5,show.signif.stars=FALSE)Load in and plot the data:
data(jsp, package="faraway")
jspr <- jsp[jsp$year==2,]
ggplot(jspr, aes(x=raven, y=math))+xlab("Raven Score")+ylab("Math Score")+geom_point(position = position_jitter())
ggplot(jspr, aes(x=social, y=math))+xlab("Social Class")+ylab("Math Score")+geom_boxplot()
Fit the model
jspr$craven <- jspr$raven-mean(jspr$raven)
mmod <- lmer(math ~ craven*social+(1|school)+(1|school:class),jspr)
print(summary(mmod), corr=FALSE)Linear mixed model fit by REML ['lmerMod']
Formula: math ~ craven * social + (1 | school) + (1 | school:class)
Data: jspr
REML criterion at convergence: 5921.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.987 -0.525 0.147 0.624 2.634
Random effects:
Groups Name Variance Std.Dev.
school:class (Intercept) 1.18 1.08
school (Intercept) 3.15 1.77
Residual 27.14 5.21
Number of obs: 953, groups: school:class, 90; school, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 31.9112 1.1955 26.69
craven 0.6058 0.1885 3.21
social2 0.0236 1.2722 0.02
social3 -0.6307 1.3089 -0.48
social4 -1.9670 1.1971 -1.64
social5 -1.3585 1.3002 -1.04
social6 -2.2687 1.3737 -1.65
social7 -2.5518 1.4055 -1.82
social8 -3.3950 1.8014 -1.88
social9 -0.8313 1.2535 -0.66
craven:social2 -0.1321 0.2058 -0.64
craven:social3 -0.2243 0.2189 -1.02
craven:social4 0.0358 0.1949 0.18
craven:social5 -0.1503 0.2089 -0.72
craven:social6 -0.0386 0.2326 -0.17
craven:social7 0.3983 0.2318 1.72
craven:social8 0.2560 0.2615 0.98
craven:social9 -0.0810 0.2055 -0.39anova(mmod)Analysis of Variance Table
Df Sum Sq Mean Sq F value
craven 1 10160 10160 374.33
social 8 610 76 2.81
craven:social 8 565 71 2.60As expected there is no inferential output.
Compute confidence intervals
confint(mmod, method="boot") 2.5 % 97.5 %
sd_(Intercept)|school:class 0.000000 1.83289
sd_(Intercept)|school 0.698365 2.36535
sigma 4.964268 5.46959
(Intercept) 29.343825 34.26939
craven 0.252737 0.95004
social2 -2.371178 2.60305
social3 -2.977756 1.71295
social4 -4.179090 0.32526
social5 -3.888741 1.13666
social6 -5.159254 0.27126
social7 -5.376949 0.16934
social8 -6.867827 0.44423
social9 -3.332531 1.64995
craven:social2 -0.503540 0.25935
craven:social3 -0.647204 0.21972
craven:social4 -0.324462 0.42358
craven:social5 -0.551636 0.25193
craven:social6 -0.532220 0.38311
craven:social7 -0.043332 0.84925
craven:social8 -0.216729 0.77682
craven:social9 -0.420089 0.36742We notice that there is some evidence that there is little variation between classes within schools.
pbkrtest package
We can test the fixed effects in this model by a sequence of comparisons. First we can test the three way interactions:
library(pbkrtest)
mmod <- lmer(math ~ raven*social*gender+(1|school)+(1|school:class),data=jspr,REML=FALSE)
mmod2 <- lmer(math ~ (raven+social+gender)^2+(1|school)+(1|school:class),data=jspr,REML=FALSE)
KRmodcomp(mmod, mmod2)F-test with Kenward-Roger approximation; computing time: 0.34 sec.
large : math ~ raven + social + gender + (1 | school) + (1 | school:class) +
raven:social + raven:gender + social:gender + raven:social:gender
small : math ~ (raven + social + gender)^2 + (1 | school) + (1 | school:class)
stat ndf ddf F.scaling p.value
Ftest 0.85 8.00 895.01 1 0.56We see that the outcome is very similar to that seen in the text although due to the imbalance we do not expect an exact match. Several additional comparisons would be necessary to test all the terms.
Several of the subsequent tests can be executed using this method. For example, we can test the compositional effect:
schraven <- lm(raven ~ school, jspr)$fit
jspr$craven <- jspr$raven-mean(jspr$raven)
mmod <- lmer(math ~ craven+social+(1|school)+(1|school:class),jspr, REML=FALSE)
mmodc <- lmer(math ~ craven+social+schraven+(1|school)+ (1|school:class), jspr,REML=FALSE)
KRmodcomp(mmodc,mmod)F-test with Kenward-Roger approximation; computing time: 0.14 sec.
large : math ~ craven + social + schraven + (1 | school) + (1 | school:class)
small : math ~ craven + social + (1 | school) + (1 | school:class)
stat ndf ddf F.scaling p.value
Ftest 0.18 1.00 55.35 1 0.68The outcome is almost the same as before.
Testing the random effects
The RLRsim package cannot be used here because there is more than one random effect term.
We can implement the parametric bootstrap directly to test whether there is any variance in the classes
First we compute the observed LRT statistic
mmod <- lmer(math ~ craven+social+(1|school)+(1|school:class),jspr, REML=FALSE)
rmod <- lmer(math ~ craven+social+(1|school),jspr, REML=FALSE)
(olrt <- 2*(logLik(mmod)-logLik(rmod)))'log Lik.' 1.8853 (df=13)Now we resample 1000 times:
lrstat <- numeric(1000)
for(i in 1:1000){
rmath <- unlist(simulate(rmod))
bmod <- refit(mmod, rmath)
smod <- refit(rmod, rmath)
lrstat[i] <- 2*(logLik(bmod)-logLik(smod))
}
mean(lrstat > olrt)[1] 0.08There is some evidence of negligible variation of classes within schools. Perhaps teachers don’t matter?
lmerTest package
library(lmerTest)
mmod <- lmer(math ~ craven*social+(1|school)+(1|school:class),jspr, REML=FALSE)
print(summary(mmod), corr=FALSE)Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
approximations to degrees of freedom [merModLmerTest]
Formula: math ~ craven * social + (1 | school) + (1 | school:class)
Data: jspr
AIC BIC logLik deviance df.resid
5949.4 6051.5 -2953.7 5907.4 932
Scaled residuals:
Min 1Q Median 3Q Max
-4.024 -0.530 0.149 0.631 2.655
Random effects:
Groups Name Variance Std.Dev.
school:class (Intercept) 1.18 1.09
school (Intercept) 3.01 1.74
Residual 26.64 5.16
Number of obs: 953, groups: school:class, 90; school, 48
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 31.9088 1.1838 892.0000 26.95 <2e-16
craven 0.6057 0.1868 922.0000 3.24 0.0012
social2 0.0252 1.2603 939.0000 0.02 0.9840
social3 -0.6275 1.2967 938.0000 -0.48 0.6285
social4 -1.9639 1.1859 942.0000 -1.66 0.0981
social5 -1.3555 1.2881 941.0000 -1.05 0.2929
social6 -2.2681 1.3609 944.0000 -1.67 0.0959
social7 -2.5514 1.3924 945.0000 -1.83 0.0672
social8 -3.3927 1.7846 934.0000 -1.90 0.0576
social9 -0.8282 1.2417 945.0000 -0.67 0.5050
craven:social2 -0.1319 0.2039 921.0000 -0.65 0.5180
craven:social3 -0.2244 0.2168 927.0000 -1.03 0.3010
craven:social4 0.0359 0.1931 922.0000 0.19 0.8525
craven:social5 -0.1500 0.2070 920.0000 -0.72 0.4689
craven:social6 -0.0385 0.2304 927.0000 -0.17 0.8673
craven:social7 0.3986 0.2296 926.0000 1.74 0.0829
craven:social8 0.2563 0.2591 922.0000 0.99 0.3229
craven:social9 -0.0808 0.2036 922.0000 -0.40 0.6917anova(mmod)Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
craven 5587 5587 1 936 209.7 <2e-16
social 569 71 8 937 2.7 0.0066
craven:social 565 71 8 925 2.7 0.0071Can test the random effects:
rand(mmod)Analysis of Random effects Table:
Chi.sq Chi.DF p.value
school 7.14 1 0.008
school:class 2.39 1 0.122but we prefer the parametric bootstrap result because of problems with asymptotic distribution of the LRT in such cases.
MCMCglmm
We use the expanded prior for the reasons established in our one random effect example.
As in the nested example, we need to make sure the classes within the schools have a unique label.
library(MCMCglmm)
jspr$classch <- paste(jspr$school,jspr$class,sep=".")
eprior <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02,alpha.V=1000),G2=list(V=1,nu=0.02,alpha.V=1000)))
bmod <- MCMCglmm(math ~ craven+social,random=~school+classch,data=jspr,verbose=FALSE,prior=eprior,pr=TRUE)
xyplot(log(bmod$VCV))
The diagnostics look OK (although others should be checked). We examine the summary:
summary(bmod)
Iterations = 3001:12991
Thinning interval = 10
Sample size = 1000
DIC: 5924
G-structure: ~school
post.mean l-95% CI u-95% CI eff.samp
school 3.38 0.974 6.7 872
~classch
post.mean l-95% CI u-95% CI eff.samp
classch 1.29 7.33e-06 3.17 660
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 27.7 25.3 30.4 1000
Location effects: math ~ craven + social
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 32.02080 29.79315 33.98047 821 <0.001
craven 0.58466 0.53016 0.65310 1000 <0.001
social2 -0.38737 -2.49112 1.80880 902 0.730
social3 -0.77981 -2.94580 1.52172 779 0.508
social4 -2.11819 -4.22201 -0.00712 892 0.056
social5 -1.40664 -3.92407 0.70140 1000 0.228
social6 -2.36371 -4.84082 0.08640 1000 0.066
social7 -3.06553 -5.66022 -0.58605 1000 0.026
social8 -3.55622 -6.81500 -0.14157 1000 0.038
social9 -0.91316 -3.10904 1.35909 1252 0.404We see that the raven effect on math scores is clear while the effect of social class is more muted. We are more confident that there is a difference between schools since the credible interval is further away from zero but we are less clear about the classes within schools difference as the credible interval goes down to negligible values. Looking at posteriors distributions for the random components in the model can help make these distinctions clearer:
We can view the posterior distributions for the SDs corresponding to the three components of variation:
library(reshape2)
ref <- data.frame(sqrt(bmod$VCV))
rdf <- melt(ref, value.name="math")
ggplot(rdf, aes(x=math,color=variable))+geom_density()
I have taken square roots to convert from variances to SDs as this is more natural. The plot shows the error SD is concentrated around 5 points - this represents the student variability on the tests. The school SD has a mode just below 2 points indicating the variability between schools is considerably less than the variability between students. Finally, the class SD looks a bit smaller still and might be negligible although we are not sure of this.
Here is the posterior for the (centered) raven effect on math scores:
ref <- data.frame(bmod$Sol)
ggplot(ref, aes(x=craven))+geom_density()
It is quite concentrated around just below 0.6 so we a quite confident regarding the size of this effect. We can also get the social class posteriors:
rdf <- melt(ref[,3:10], value.name="math", variable.name="social")
ggplot(rdf, aes(x=math, color=social))+geom_density()
We can see that the differences from social class 1 for classes 4,5,6,7 and 8 are almost clearly negative but there is more overlap for 2,3 and 9. We can see that the worst difference between classes is about 4 points - slightly less than the student SD noticed earlier.
We can also examine the school posteriors:
rdf <- melt(ref[,11:58], value.name="math", variable.name="school")
ggplot(rdf, aes(x=math, group=school))+geom_density()
We have not attempted to distinguish the schools because there are too many to label sepearately. We can see that while there is some clear variation between schools, there is a substantial amount of overlap meaning that any “league table” will only moderately distinguish schools. But since people want to see it, here is our league table formed from the posterior means:
sort(colMeans(ref[,11:58]))school.29 school.28 school.21 school.1 school.45 school.40 school.9
-2.381192 -2.361399 -2.352584 -2.163779 -2.074173 -2.027070 -1.943258
school.4 school.22 school.39 school.35 school.33 school.5 school.14
-1.204783 -1.042523 -1.011145 -0.975269 -0.966785 -0.694506 -0.682613
school.48 school.12 school.2 school.20 school.42 school.18 school.16
-0.610552 -0.497106 -0.429663 -0.298724 -0.227178 -0.217221 -0.206756
school.23 school.44 school.13 school.49 school.32 school.6 school.46
-0.179811 -0.103258 0.037598 0.047559 0.060598 0.140163 0.222069
school.8 school.15 school.27 school.17 school.50 school.37 school.11
0.326869 0.372703 0.385509 0.421031 0.487143 0.514095 0.528140
school.19 school.3 school.26 school.30 school.25 school.7 school.47
0.652933 0.800845 1.086979 1.087040 1.143102 1.209974 1.567239
school.41 school.24 school.31 school.38 school.36 school.34
1.620961 2.107995 2.272928 2.281633 2.415627 2.486207 Note that there has been some renumbering of the schools due to some missing numbers in the sequence. Some care is necessary in identifying which school is which.
Check on versions of R packages used:
sessionInfo()R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] reshape2_1.2.2 MCMCglmm_2.21 ape_3.1-4 coda_0.16-1
[5] lattice_0.20-29 lmerTest_2.0-20 pbkrtest_0.4-1 MASS_7.3-33
[9] ggplot2_0.9.3.1 lme4_1.1-7 Rcpp_0.11.3 Matrix_1.1-4
[13] knitr_1.6 rmarkdown_0.2.68
loaded via a namespace (and not attached):
[1] bitops_1.0-6 boot_1.3-11 caTools_1.16
[4] cluster_1.15.2 colorspace_1.2-4 corpcor_1.6.7
[7] dichromat_2.0-0 digest_0.6.4 evaluate_0.5.3
[10] formatR_0.10 Formula_1.1-1 gdata_2.13.3
[13] gplots_2.13.0 grid_3.1.1 gtable_0.1.2
[16] gtools_3.3.1 Hmisc_3.14-3 htmltools_0.2.4
[19] KernSmooth_2.23-12 labeling_0.2 latticeExtra_0.6-26
[22] minqa_1.2.3 munsell_0.4.2 nlme_3.1-117
[25] nloptr_1.0.4 numDeriv_2012.9-1 parallel_3.1.1
[28] plyr_1.8.1 proto_0.3-10 RColorBrewer_1.0-5
[31] scales_0.2.3 splines_3.1.1 stringr_0.6.2
[34] survival_2.37-7 tensorA_0.36 tools_3.1.1
[37] yaml_2.1.11