First of all, be warned, the terminology surrounding multilevel models is vastly inconsistent. For example, multilevel models themselves may be referred to as hierarchical linear models, random effects models, multilevel models, random intercept models, random slope models, or pooling models. Depending on the discipline, software used, and the academic literature many of these terms may be referring to the same general modeling strategy. In this tutorial I will attempt to provide a user guide to multilevel modeling by demonstrating how to fit multilevel models in R and by attempting to connect the model fitting procedure to commonly used terminology used regarding these models.
We will cover the following topics:
The structure and methods of merMod objects
Extracting random effects of merMod objects
Plotting and interpreting merMod objects
If you haven’t already, make sure you head over to the Getting Started With Multilevel Models tutorial in order to ensure you have set up your environment correctly and installed all the necessary packages. The tl;dr is that you will need:
Multilevel models are appropriate for a particular kind of data structure where units are nested within groups (generally 5+ groups) and where we want to model the group structure of the data. We will use data from Jon Starkweather at the University of North Texas on student social attributes within schools and classes. Visit the excellent tutorial available here for more.
library(lme4) # load library
library(arm) # convenience functions for regression in R
Here we see we have two possible grouping variables – class and school. Let’s explore them a bit further:
table(lmm.data$class)
##
## a b c d
## 300 300 300 300
table(lmm.data$school)
##
## I II III IV V VI
## 200 200 200 200 200 200
table(lmm.data$class, lmm.data$school)
##
## I II III IV V VI
## a 50 50 50 50 50 50
## b 50 50 50 50 50 50
## c 50 50 50 50 50 50
## d 50 50 50 50 50 50
This is a perfectly balanced dataset. In all likelihood you aren’t working with a perfectly balanced dataset, but we’ll explore the implications for that in the future. For now, let’s plot the data a bit. Using the excellent xyplot function in the lattice package, we can explore the relationship between schools and classes across our variables.
require(lattice)
xyplot(extro ~ open + social + agree | class, data = lmm.data,
By school we see that students are tightly grouped, but that school I and school VI show substantially more dispersion than the other schools. The same pattern among our predictors holds between schools as it did between classes. Let’s put it all together:
xyplot(extro ~ open + social + agree | school + class, data = lmm.data,
In the last tutorial we fit a series of random intercept models to our nested data. We will examine the lmerMod object produced when we fit this model in much more depth in order to understand how to work with mixed effect models in R. We start by fitting a the basic example below grouped by class:
MLexamp1 <- lmer(extro ~ open + agree + social + (1|school), data=lmm.data)
class(MLexamp1)
## [1] "lmerMod"
## attr(,"package")
## [1] "lme4"
First, we see that MLexamp1 is now an R object of the class lmerMod. This lmerMod object is an S4 class, and to explore its structure, we use slotNames:
Within the lmerMod object we see a number of objects that we may wish to explore. To look at any of these, we can simply type MLexamp1@ and then the slot name itself. For example:
MLexamp1@call # returns the model call
## lmer(formula = extro ~ open + agree + social + (1 | school),
## data = lmm.data)
MLexamp1@beta # returns the betas
## [1] 59.116514 0.009751 0.027788 -0.002151
class(MLexamp1@frame) # returns the class for the frame slot
## [1] "data.frame"
head(MLexamp1@frame) # returns the model frame
## extro open agree social school
## 1 63.69 43.43 38.03 75.06 IV
## 2 69.48 46.87 31.49 98.13 VI
## 3 79.74 32.27 40.21 116.34 VI
## 4 62.97 44.41 30.51 90.47 IV
## 5 64.25 36.86 37.44 98.52 IV
## 6 50.97 46.26 38.83 75.22 I
The merMod object has a number of methods available – too many to enumerate here. But, we will go over some of the more common in the list below:
A common need is to extract the fixed effects from a merMod object. fixef extracts a named numeric vector of the fixed effects, which is handy.
fixef(MLexamp1)
## (Intercept) open agree social
## 59.116514 0.009751 0.027788 -0.002151
If you want to get a sense of the p-values or statistical significance of these parameters, first consult the lme4 help by running ?mcmcsamp for a rundown of various ways of doing this. One convenient way built into the package is:
confint(MLexamp1, level = 0.99)
## 0.5 % 99.5 %
## .sig01 4.91840 23.88758
## .sigma 2.53287 2.81456
## (Intercept) 46.27751 71.95610
## open -0.02465 0.04415
## agree -0.01164 0.06721
## social -0.01493 0.01063
From here we can see first that our fixed effect parameters overlap 0 indicating no evidence of an effect. We can also see that .sig01, which is our estimate of the variability in the random effect, is very large and very widely defined. This indicates we may have a lack of precision between our groups - either because the group effect is small between groups, we have too few groups to get a more precise estimate, we have too few units within each group, or a combination of all of the above.
Another common need is to extract the residual standard error, which is necessary for calculating effect sizes. To get a named vector of the residual standard error:
sigma(MLexamp1)
## [1] 2.671
For example, it is common practice in education research to standardize fixed effects into “effect sizes” by dividing the fixed effect paramters by the residual standard error, which can be accomplished in lme4 easily:
fixef(MLexamp1) / sigma(MLexamp1)
## (Intercept) open agree social
## 22.1336705 0.0036508 0.0104042 -0.0008055
From this, we can see that our predictors of openness, agreeableness and social are virtually useless in predicting extroversion – as our plots showed. Let’s turn our attention to the random effects next.
Explore Group Variation and Random Effects
In all likelihood you fit a mixed-effect model because you are directly interested in the group-level variation in your model. It is not immediately clear how to explore this group level variation from the results of summary.merMod. What we get from this output is the variance and the standard deviation of the group effect, but we do not get effects for individual groups. This is where the ranef function comes in handy.
ranef(MLexamp1)
## $school
## (Intercept)
## I -14.091
## II -6.183
## III -1.971
## IV 1.966
## V 6.331
## VI 13.948
Running the ranef function gives us the intercepts for each school, but not much additional information – for example the precision of these estimates. To do that, we need some additional commands:
re1 <- ranef(MLexamp1, condVar=TRUE) # save the ranef.mer object
class(re1)
## [1] "ranef.mer"
attr(re1[[1]], which = "postVar")
## , , 1
##
## [,1]
## [1,] 0.03565
##
## , , 2
##
## [,1]
## [1,] 0.03565
##
## , , 3
##
## [,1]
## [1,] 0.03565
##
## , , 4
##
## [,1]
## [1,] 0.03565
##
## , , 5
##
## [,1]
## [1,] 0.03565
##
## , , 6
##
## [,1]
## [1,] 0.03565
The ranef.mer object is a list which contains a data.frame for each group level. The dataframe contains the random effects for each group (here we only have an intercept for each school). When we ask lme4 for the conditional variance of the random effects it is stored in an attribute of those dataframes as a list of variance-covariance matrices.
This structure is indeed complicated, but it is powerful as it allows for nested, grouped, and cross-level random effects. Also, the creators of lme4 have provided users with some simple shortcuts to get what they are really interested in out of a ranef.mer object.
This graphic shows a dotplot of the random effect terms, also known as a caterpillar plot. Here you can clearly see the effects of each school on extroversion as well as their standard errors to help identify how distinct the random effects are from one another. Interpreting random effects is notably tricky, but for assistance I would recommend looking at a few of these resources:
Gelman and Hill 2006 - Data Analysis Using Regression and Multilevel/Hierarchical Techniques
John Fox - An R Compantion to Applied Regression web appendix
Using Simulation and Plots to Explore Random Effects
A common econometric approach is to create what are known as empirical Bayes estimates of the group-level terms. Unfortunately there is not much agreement about what constitutes a proper standard error for the random effect terms or even how to consistently define empirical Bayes estimates.1 However, in R there are a few additional ways to get estimates of the random effects that provide the user with information about the relative sizes of the effects for each unit and the precision in that estimate. To do this, we use the sim function in the arm package.2
# A function to extract simulated estimates of random effect paramaters from
# lme4 objects using the sim function in arm
# whichel = the character for the name of the grouping effect to extract estimates for
# nsims = the number of simulations to pass to sim
# x = model object
REsim <- function(x, whichel=NULL, nsims){
require(plyr)
mysim <- sim(x, n.sims = nsims)
if(missing(whichel)){
dat <- plyr::adply(mysim@ranef[[1]], c(2, 3), plyr::each(c(mean, median, sd)))
warning("Only returning 1st random effect because whichel not specified")
} else{
dat <- plyr::adply(mysim@ranef[[whichel]], c(2, 3), plyr::each(c(mean, median, sd)))
}
return(dat)
}
REsim(MLexamp1, whichel = "school", nsims = 1000)
## X1 X2 mean median sd
## 1 I (Intercept) -14.053 -14.093 3.990
## 2 II (Intercept) -6.141 -6.122 3.988
## 3 III (Intercept) -1.934 -1.987 3.981
## 4 IV (Intercept) 2.016 2.051 3.993
## 5 V (Intercept) 6.378 6.326 3.981
## 6 VI (Intercept) 13.992 14.011 3.987
The REsim function returns for each school the level name X1, the estimate name, X2, the mean of the estimated values, the median, and the standard deviation of the estimates.
Another convenience function can help us plot these results to see how they compare to the results of dotplot:
# Dat = results of REsim
# scale = factor to multiply sd by
# var = character of "mean" or "median"
# sd = character of "sd"
plotREsim <- function(dat, scale, var, sd){
require(eeptools)
dat[, sd] <- dat[, sd] * scale
dat[, "ymax"] <- dat[, var] + dat[, sd]
dat[, "ymin"] <- dat[, var] - dat[, sd]
dat[order(dat[, var]), "id"] <- c(1:nrow(dat))
ggplot(dat, aes_string(x = "id", y = var, ymax = "ymax",
ymin = "ymin")) +
geom_pointrange() + theme_dpi() +
labs(x = "Group", y = "Effect Range", title = "Effect Ranges") +
This presents a more conservative view of the variation between random effect components. Depending on how your data was collected and your research question, alternative ways of estimating these effect sizes are possible. However, proceed with caution.3
Another approach recommended by the authors of lme4 involves the RLRsim package. Using this package we can test whether or not inclusion of the random effects improves model fit and we can evaluate the p-value of additional random effect terms using a likelihood ratio test based on simulation.4
library(RLRsim)
m0 <- lm(extro ~ agree + open + social, data =lmm.data) # fit the null model
exactLRT(m = MLexamp1, m0 = m0)
##
## simulated finite sample distribution of LRT. (p-value based on
## 10000 simulated values)
##
## data:
## LRT = 2958, p-value < 2.2e-16
Here exactLRT issues a warning because we originally fit the model with REML instead of full maximum likelihood. Fortunately, the refitML function in lme4 allows us to easily refit our model using full maximum likelihood to conduct an exact test easily.
mA <- refitML(MLexamp1)
exactLRT(m= mA, m0 = m0)
##
## simulated finite sample distribution of LRT. (p-value based on
## 10000 simulated values)
##
## data:
## LRT = 2958, p-value < 2.2e-16
Here we can see that the inclusion of our grouping variable is significant, even though the effect of each individual group may be substantively small and/or imprecisely measured. This is important in understanding the correct specification of the model. Our next tutorial will cover specification tests like this in more detail.
How do interpret the substantive impact of our random effects? This is often critical in observation work trying to use a multilevel structure to understand the impact that the grouping can have on the individual observation. To do this we select 12 random cases and then we simulate their predicted value of extro if they were placed in each of the 6 schools. Note, that this is a very simple simulation just using the mean of the fixed effect and the conditional mode of the random effect and not replicating or sampling to get a sense of the variability. This will be left as an exercise to the reader and/or a future tutorial!
# Simulate
# Let's create 12 cases of students
#
#sample some rows
simX <- sample(lmm.data$id, 12)
simX <- lmm.data[lmm.data$id %in% simX, c(3:5)] # get their data
Here we can clearly see that within each school the cases are relatively the same indicating that the group effect is larger than the individual effects.
These plots are useful in demonstrating the relative importance of group and individual effects in a substantive fashion. Even more can be done to make the the graphs more informative, such as placing references to the total variability of the outcome and also looking at the distance moving groups moves each observation from its true value.
Conclusion
lme4 provides a very powerful object-oriented toolset for dealing with mixed effect models in R. Understanding model fit and confidence intervals of lme4 objects requires some diligent research and the use of a variety of functions and extensions of lme4 itself. In our next tutorial we will explore how to identify a proper specification of a random-effect model and Bayesian extensions of the lme4 framework for difficult to specify models. We will also explore the generalized linear model framework and the glmer function for generalized linear modeling with multi-levels.
Appendix
print(sessionInfo(),locale=FALSE)
## R version 3.1.0 (2014-04-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
[See message from lme4 co-author Doug Bates on this subject]. (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q4/002984.html)↩
Andrew Gelman and Yu-Sung Su (2014). arm: Data Analysis Using Regression and Multilevel/Hierarchical Models. R package version 1.7-03. http://CRAN.R-project.org/package=arm↩
WikiDot FAQ from the R Mixed Models Mailing List↩
There are also an extensive series of references available in the References section of the help by running ?exactLRT and ?exactRLRT.↩
发帖
雷达卡
京公网安备 11010802022788号
