## (C) (cc by-sa) Wouter van Atteveldt & Jan Kleinnijenhuis, file generated juni 24 2015Multilevel Modeling with R
*Caveat* I am not an expert on multilevel modelling. Please take this document as a source of inspiration rather than as a definitive set of answers.This document gives an overview of two commonly used packages for multi-level modelling in R (also called ‘mixed models’ or ‘random effects models’).
Since the yearly time series data we used so far are not suitable for multilevel analysis, let’s take the textbook data of Joop Hox on popularity of pupils in schools: (see also http://www.ats.ucla.edu/stat/examples/ma_hox/)
library(foreign)
popdata<-read.dta("http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta")
head(popdata)## pupil school popular sex texp const teachpop
## 1 1 1 8 girl 24 1 7
## 2 2 1 7 boy 24 1 7
## 3 3 1 7 girl 24 1 6
## 4 4 1 9 girl 24 1 6
## 5 5 1 8 girl 24 1 7
## 6 6 1 7 boy 24 1 7Now, we can model a time series model with only the random intercept at the school level:
library(nlme)
m = lme(popular ~ sex + texp, random=~1|school, popdata)
summary(m)## Linear mixed-effects model fit by REML
## Data: popdata
## AIC BIC logLik
## 4454.358 4482.355 -2222.179
##
## Random effects:
## Formula: ~1 | school
## (Intercept) Residual
## StdDev: 0.6971072 0.6781765
##
## Fixed effects: popular ~ sex + texp
## Value Std.Error DF t-value p-value
## (Intercept) 3.560682 0.17147861 1899 20.764582 0
## sexgirl 0.844669 0.03095057 1899 27.290909 0
## texp 0.093445 0.01085419 98 8.609134 0
## Correlation:
## (Intr) sexgrl
## sexgirl -0.088
## texp -0.905 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.35852403 -0.67969990 0.02435564 0.59332190 3.78514578
##
## Number of Observations: 2000
## Number of Groups: 100So, popularity of a course is determined by both gender and teacher experience. Let’s try a varying slopes model, with teacher experience also differing per school, and see whether that is a significant improvement:
m2 = lme(popular ~ sex + texp, random=~texp|school, popdata)
anova(m, m2)## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 5 4454.358 4482.355 -2222.179
## m2 2 7 4456.443 4495.639 -2221.222 1 vs 2 1.915354 0.3838So, although the log likelihood of m2 is slightly better, it also uses more degrees of freedom and the BIC is higher, indicating a worse model. The anova output means that this change is not significant.
Next, let’s have a look at the slope of the gender effect. First, a useful tool can be a visual inspection of the slope for a random sample of schools, just to get an idea of variation. First, take a random sample of 12 schools from the list of unique school ids:
schools = sample(unique(popdata$school), size=12, replace=F)
sample = popdata[popdata$school %in% schools, ]Now, we can use the xyplot function from the lattice package:
library(lattice)
xyplot(popular~sex|as.factor(school),type=c("p","g","r"), col.line="black", data=sample)
So, (at least in my sample) there is considerable variation: in some schools gender has almost no effect, but in other schools the slope is relatively steep and generally positive (meaning girls have higher popularity). Let’s test whether a model with a random slope on gender is a significant improvement:
library(texreg)## Version: 1.35
## Date: 2015-04-25
## Author: Philip Leifeld (University of Konstanz)
##
## Please cite the JSS article in your publications -- see citation("texreg").m2 = lme(popular ~ sex + texp, random=~sex|school, popdata)
anova(m, m2)## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 5 4454.358 4482.355 -2222.179
## m2 2 7 4289.895 4329.091 -2137.948 1 vs 2 168.4634 <.0001screenreg(list(m, m2))##
## ==========================================
## Model 1 Model 2
## ------------------------------------------
## (Intercept) 3.56 *** 3.34 ***
## (0.17) (0.16)
## sexgirl 0.84 *** 0.84 ***
## (0.03) (0.06)
## texp 0.09 *** 0.11 ***
## (0.01) (0.01)
## ------------------------------------------
## AIC 4454.36 4289.89
## BIC 4482.36 4329.09
## Log Likelihood -2222.18 -2137.95
## Num. obs. 2000 2000
## Num. groups 100 100
## ==========================================
## *** p < 0.001, ** p < 0.01, * p < 0.05So, m2 is indeed a significant improvement.
The lme4 package
lme4 is a package that gives a bit more flexibility in specifying time series. Specifically, it allows us to specify a binomial family, i.e. logistic regression. The following dichotomized the popularity and checks whether the effect of gender on popularity is dependent on the school:
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:base':
##
## crossprod, tcrossprod
##
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
## lmListpopdata$dich = cut(popdata$popular, 2, labels=c("lo","hi"))
m = glmer(dich ~ sex + (1|school), popdata, family="binomial")
m2 = glmer(dich ~ sex + (1 + sex|school), popdata, family="binomial")
summary(m2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: dich ~ sex + (1 + sex | school)
## Data: popdata
##
## AIC BIC logLik deviance df.resid
## 1536.2 1564.2 -763.1 1526.2 1995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7807 -0.3262 -0.1125 0.2658 5.3141
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 8.495 2.915
## sexgirl 3.400 1.844 0.16
## Number of obs: 2000, groups: school, 100
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0970 0.3512 -5.970 2.37e-09 ***
## sexgirl 2.8582 0.2969 9.626 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sexgirl -0.215anova(m, m2)## Data: popdata
## Models:
## m: dich ~ sex + (1 | school)
## m2: dich ~ sex + (1 + sex | school)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m 3 1572.0 1588.8 -783.00 1566.0
## m2 5 1536.2 1564.2 -763.12 1526.2 39.757 2 2.328e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1If we would like to see for which schools the effect of gender were the strongest, we can use the ranef function to get the intercepts and slopes per group, and order them by slope:
effects = ranef(m2)$school
effects = effects[order(effects$sexgirl), ]
head(effects)## (Intercept) sexgirl
## 38 1.9159645 -3.513633
## 53 0.5027642 -3.224843
## 10 -1.2437362 -2.345330
## 80 1.9168907 -2.131811
## 12 0.1724678 -2.115681
## 61 -0.6707056 -1.896863tail(effects)## (Intercept) sexgirl
## 32 -0.7054251 1.752564
## 8 0.4938311 1.757972
## 7 0.1999551 1.798121
## 16 0.5681741 1.872205
## 54 -0.4788237 2.131774
## 52 -0.4730479 2.361743