|
|
|
||||
| [stat/R/examples/mlm_imm/header.htm] |
|
||||
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 http://www.ats.ucla.edu/stat/examples/imm/imm10 and imported into R using the foreign library from a directory called rdata on the local computer.
library(foreign) imm10 <- read.dta(file="~/rdata/imm10.dta") attach(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
amat<-cbind(tmat,pmat)
amat
(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), method="ML")
Linear mixed-effects model fit by maximum likelihood
Formula: math ~ homework + (homework | schnum)
AIC BIC logLik MLdeviance REMLdeviance
1779 1797 -884.7 1769 1764
Random effects:
Groups Name Variance Std.Dev. Corr
schnum (Intercept) 62.195 7.8864
homework 19.982 4.4702 -0.803
Residual 43.051 6.5613
number of obs: 260, groups: schnum, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 44.772 2.611 17.150
homework 2.049 1.472 1.392
Correlation of Fixed Effects:
(Intr)
homework -0.802
Equation near bottom of page 49 and Table 3.4.
lmer(math ~ homework + public + (homework|schnum), method="ML")
Linear mixed-effects model fit by maximum likelihood
Formula: math ~ homework + public + (homework | schnum)
AIC BIC logLik MLdeviance REMLdeviance
1763 1784 -875.4 1751 1744
Random effects:
Groups Name Variance Std.Dev. Corr
schnum (Intercept) 40.669 6.3772
homework 21.689 4.6571 -0.982
Residual 42.956 6.5541
number of obs: 260, groups: schnum, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 58.056 2.694 21.552
homework 1.941 1.525 1.272
public -14.651 1.830 -8.005
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), method="ML")
Linear mixed-effects model fit by maximum likelihood
Formula: math ~ homework + public + homework:public + (homework | schnum)
AIC BIC logLik MLdeviance REMLdeviance
1766 1791 -876.1 1752 1742
Random effects:
Groups Name Variance Std.Dev. Corr
schnum (Intercept) 40.046 6.3282
homework 22.576 4.7514 -1.000
Residual 44.021 6.6348
number of obs: 260, groups: schnum, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 59.2102 6.5686 9.014
homework 1.0946 4.7750 0.229
public -15.8672 6.9458 -2.284
homework:public 0.8899 5.0494 0.176
Correlation of Fixed Effects:
(Intr) homwrk public
homework -0.982
public -0.946 0.929
homwrk:pblc 0.929 -0.946 -0.981
Warning message:
In .local(x, ..., value) :
Estimated variance-covariance for factor 'schnum' is singular
Table 3.6, page 52.
# quietly rerun model from page 47 m1<-lmer(math ~ homework + (homework|schnum), method="ML") # random intercepts and slopes for each school rcoef<-coef(m1) rcoef An object of class "coef.lmer" [[1]] (Intercept) homework 1 50.27035 -3.143235 2 48.88295 -2.753085 3 39.19630 7.566027 4 35.14984 5.395286 5 53.08193 -3.737538 6 48.58528 -1.765018 7 58.06732 1.332620 8 37.14563 6.061955 9 39.16759 5.430831 10 38.17021 6.102127
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)}
UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services