Note: This page is designed to show the how multilevel model can be done using R and to be able to compare the results with those in the book.
On this page we will use the lmer function which is found in the lme4 package. There are several other possible choices but we will go with lmer.
The data were downloaded in Stata format from here and imported into R using the foreign library from a directory called rdata on the local computer. This page is updated using R 2.11.1 in January, 2011.
library(foreign) imm10
Table 3.2, page 46. OLS regression lines over 10 schools.
tmp<-by(imm10, schnum, function(x) lm(math ~ homework, data=x)) tmat<-t(sapply(tmp,coef)) tmat (Intercept) homework 1 50.68354 -3.553797 2 49.01229 -2.920123 3 38.75000 7.909091 4 34.39382 5.592664 5 53.93863 -4.718412 6 49.25896 -2.486056 7 59.21022 1.094640 8 36.05535 6.496310 9 38.52000 5.860000 10 37.71392 6.335052
Two equations at the top of page 47.
# get value of public for each school ptmp<-by(imm10, schnum, function(x) lm(public~1, data=x)) pmat<-as.matrix(sapply(ptmp,coef)) # combine public with intercept and slope a<-cbind(tmat,pmat) a (Intercept) homework 1 50.68354 -3.553797 1 2 49.01229 -2.920123 1 3 38.75000 7.909091 1 4 34.39382 5.592664 1 5 53.93863 -4.718412 1 6 49.25896 -2.486056 1 7 59.21022 1.094640 0 8 36.05535 6.496310 1 9 38.52000 5.860000 1 10 37.71392 6.335052 1 summary(lm(a[,1]~a[,3])) Call: lm(formula = a[, 1] ~ a[, 3]) Residuals: Min 1Q Median 3Q Max -8.754 -5.232 -2.199 6.050 10.791 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 59.210 7.435 7.964 4.51e-05 a[, 3] -16.063 7.837 -2.050 0.0745 Residual standard error: 7.435 on 8 degrees of freedom Multiple R-squared: 0.3443, Adjusted R-squared: 0.2624 F-statistic: 4.201 on 1 and 8 DF, p-value: 0.07455 summary(lm(a[,2]~a[,3])) Call: lm(formula = a[, 2] ~ a[, 3]) Residuals: Min 1Q Median 3Q Max -6.776 -4.869 1.768 4.159 5.852 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.0946 5.2680 0.208 0.841 a[, 3] 0.9626 5.5530 0.173 0.867 Residual standard error: 5.268 on 8 degrees of freedom Multiple R-squared: 0.003742, Adjusted R-squared: -0.1208 F-statistic: 0.03005 on 1 and 8 DF, p-value: 0.8667
Equation near bottom of page 47 and Table 3.3.
library(lme4) lmer(math ~ homework + (homework|schnum), REML=FALSE) Linear mixed model fit by maximum likelihood Formula: math ~ homework + (homework | schnum) AIC BIC logLik deviance REMLdev 1781 1803 -884.7 1769 1764 Random effects: Groups Name Variance Std.Dev. Corr schnum (Intercept) 61.806 7.8617 homework 19.979 4.4697 -0.804 Residual 43.067 6.5625 Number of obs: 260, groups: schnum, 10 Fixed effects: Estimate Std. Error t value (Intercept) 44.773 2.603 17.199 homework 2.049 1.472 1.391 Correlation of Fixed Effects: (Intr) homework -0.803
Equation near bottom of page 49 and Table 3.4.
lmer(math ~ homework + public + (homework|schnum), REML=FALSE) Linear mixed model fit by maximum likelihood Formula: math ~ homework + public + (homework | schnum) AIC BIC logLik deviance REMLdev 1765 1790 -875.4 1751 1744 Random effects: Groups Name Variance Std.Dev. Corr schnum (Intercept) 40.677 6.3778 homework 21.683 4.6565 -0.982 Residual 42.955 6.5540 Number of obs: 260, groups: schnum, 10 Fixed effects: Estimate Std. Error t value (Intercept) 58.056 2.695 21.546 homework 1.941 1.525 1.273 public -14.651 1.831 -8.000 Correlation of Fixed Effects: (Intr) homwrk homework -0.772 public -0.602 0.010
Equation at the bottom of page 50 and Table 3.5. The negative value for the interaction coefficient in the book is probably a typo error, it should be positive.
lmer(math ~ homework + public + homework:public + (homework|schnum), REML=FALSE) Linear mixed model fit by maximum likelihood Formula: math ~ homework + public + homework:public + (homework | schnum) AIC BIC logLik deviance REMLdev 1767 1795 -875.4 1751 1739 Random effects: Groups Name Variance Std.Dev. Corr schnum (Intercept) 40.503 6.3642 homework 21.577 4.6451 -0.982 Residual 42.954 6.5540 Number of obs: 260, groups: schnum, 10 Fixed effects: Estimate Std. Error t value (Intercept) 59.2102 6.5976 8.975 homework 1.0946 4.6688 0.234 public -15.9419 6.9775 -2.285 homework:public 0.9472 4.9385 0.192 Correlation of Fixed Effects: (Intr) homwrk public homework -0.966 public -0.946 0.913 homwrk:pblc 0.913 -0.945 -0.965
Table 3.6, page 52.
# quietly rerun model from page 47 m1<-lmer(math ~ homework + (homework|schnum), REML=FALSE) # random intercepts and slopes for each school rcoef<-coef(m1) rcoef
$schnum (Intercept) homework 1 50.27101 -3.143116 2 48.88642 -2.753799 3 39.19512 7.566056 4 35.15670 5.393500 5 53.08124 -3.737865 6 48.58612 -1.765064 7 58.05519 1.336116 8 37.15232 6.059531 9 39.16960 5.430173 10 38.17276 6.101294
Figure 3.8, page 53.
n<-nrow(rcoef[[1]]) x<-cbind(0, 7) y<-cbind(seq(1:n), seq(1:n)) y[,1]<-rcoef[[1]][,1] y[,2]<-rcoef[[1]][,1]+rcoef[[1]][,2]*7 plot(x, y[1,], type="l", ylab="predicte", ylim=c(20,100)) for (i in 2:n){points(x, y[i,], type="l", lty=i)}