## (C) (cc by-sa) Wouter van Atteveldt & Jan Kleinnijenhuis, file generated juni 24 2015

# Multilevel 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)
##   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        7

Now, 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: 100

So, 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.3838

So, 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  <.0001
screenreg(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.05

So, 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)
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:base':
##
##     crossprod, tcrossprod
##
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
##     lmList
popdata$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.215
anova(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 ' ' 1

If 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), ]
##    (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.896863
tail(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