Subset of data from the Television School and Family Smoking Prevention and Cessation Project (TVSFP).: http://www.hsph.harvard.edu/fitzmaur/ala2e/tvsfp.txt
Linear Mixed-Effects Models Using R: http://www.springer.com/statistics/statistical+theory+and+methods/book/978-1-4614-3899-1
Students are nested within classrooms (cid), and classrooms are nested within schools (sid).
## Load
library(foreign)
tvsfp <- read.dta("http://www.hsph.harvard.edu/fitzmaur/ala2e/tvsfp.dta")
## Summary
summary(tvsfp)
sid cid curriculum tvprevent prescore postscore
Min. :193 Min. :193101 Min. :0.000 Min. :0.000 Min. :0.00 Min. :0.00
1st Qu.:405 1st Qu.:405101 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:1.00 1st Qu.:2.00
Median :415 Median :415105 Median :0.000 Median :0.000 Median :2.00 Median :3.00
Mean :422 Mean :422042 Mean :0.477 Mean :0.499 Mean :2.07 Mean :2.66
3rd Qu.:509 3rd Qu.:509106 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:3.00 3rd Qu.:4.00
Max. :515 Max. :515113 Max. :1.000 Max. :1.000 Max. :6.00 Max. :7.00
## Load nlme
library(nlme)
## Fit
## outcome2 ~ outcome1 + tx.1 + tx.2 + interaction
res.lme <- lme(fixed = postscore ~ prescore + curriculum + tvprevent + curriculum:tvprevent,
data = tvsfp,
random = ~ 1 | sid/cid, # random intercepts level3 = sid, level2 = cid
method = "REML"
)
summary(res.lme)
Linear mixed-effects model fit by REML
Data: tvsfp
AIC BIC logLik
5389 5432 -2687
Random effects:
Formula: ~1 | sid
(Intercept)
StdDev: 0.1966
Formula: ~1 | cid %in% sid
(Intercept) Residual
StdDev: 0.2543 1.266
Fixed effects: postscore ~ prescore + curriculum + tvprevent + curriculum:tvprevent
Value Std.Error DF t-value p-value
(Intercept) 1.7020 0.12543 1464 13.569 0.0000
prescore 0.3054 0.02589 1464 11.794 0.0000
curriculum 0.6413 0.16095 24 3.985 0.0005
tvprevent 0.1821 0.15724 24 1.158 0.2583
curriculum:tvprevent -0.3309 0.22459 24 -1.474 0.1536
Correlation:
(Intr) prescr crrclm tvprvn
prescore -0.442
curriculum -0.634 0.015
tvprevent -0.645 0.008 0.501
curriculum:tvprevent 0.448 0.005 -0.716 -0.700
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.49875 -0.69757 -0.01721 0.68241 3.14602
Number of Observations: 1600
Number of Groups:
sid cid %in% sid
28 135
intervals(res.lme)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 1.4559 1.7020 1.9480
prescore 0.2546 0.3054 0.3562
curriculum 0.3091 0.6413 0.9735
tvprevent -0.1424 0.1821 0.5066
curriculum:tvprevent -0.7945 -0.3309 0.1326
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: sid
lower est. upper
sd((Intercept)) 0.1034 0.1966 0.3738
Level: cid
lower est. upper
sd((Intercept)) 0.165 0.2543 0.392
Within-group standard error:
lower est. upper
1.221 1.266 1.312
## Load lme4
library(lme4)
## group IDs need to be in factors for lme4
tvsfp <- within(tvsfp, {
sid <- factor(sid)
cid <- factor(cid)
})
## Fit
## outcome2 ~ outcome1 + tx.1 + tx.2 + interaction + (1|level3/level2)
res.lmer <- lmer(formula = postscore ~ prescore + curriculum + tvprevent + curriculum:tvprevent + (1|sid/cid),
data = tvsfp,
family = gaussian(link = "identity")
)
res.lmer
Linear mixed model fit by REML
Formula: postscore ~ prescore + curriculum + tvprevent + curriculum:tvprevent + (1 | sid/cid)
Data: tvsfp
AIC BIC logLik deviance REMLdev
5389 5432 -2687 5358 5373
Random effects:
Groups Name Variance Std.Dev.
cid:sid (Intercept) 0.0647 0.254
sid (Intercept) 0.0386 0.197
Residual 1.6023 1.266
Number of obs: 1600, groups: cid:sid, 135; sid, 28
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.7020 0.1254 13.57
prescore 0.3054 0.0259 11.79
curriculum 0.6413 0.1609 3.99
tvprevent 0.1821 0.1572 1.16
curriculum:tvprevent -0.3309 0.2245 -1.47
Correlation of Fixed Effects:
(Intr) prescr crrclm tvprvn
prescore -0.442
curriculum -0.634 0.015
tvprevent -0.645 0.008 0.501
crrclm:tvpr 0.448 0.005 -0.716 -0.700
The total variance of the outcome was partitioned into three components:
The within-classroom correlation of student outcomes is (0.039 + 0.065)/(0.065 + 0.039 + 1.602) = 0.061
The within-school correlation of student outcomes (different classrooms) is (0.039)/(0.065 + 0.039 + 1.602) = 0.023