Working with generalized linear mixed models (GLMM) and linear mixed models (LMM) has become increasingly easy with the advances in the lme4 package recently. As we have found ourselves using these models more and more within our work, we, the authors, have developed a set of tools for simplifying and speeding up common tasks for interacting with merMod objects from lme4. This package provides those tools.
As the complexity of the model fit grows, it becomes harder and harder to interpret the substantive effect of parameters in the model.
Let’s start with a medium-sized example model using the InstEval data provided by the lme4 package. These data represent university lecture evaluations at ETH Zurich made by students. In this data, s is an individual student, d is an individual lecturer, studage is the semester the student is enrolled, lectage is how many semesters back the lecture with the rating took place, dept is the department of the lecture, and y is an integer 1:5 representing the ratings of the lecture from “poor” to “very good”:
library(lme4)
head(InstEval)
#> s d studage lectage service dept y
#> 1 1 1002 2 2 0 2 5
#> 2 1 1050 2 1 1 6 2
#> 3 1 1582 2 2 0 2 5
#> 4 1 2050 2 2 1 3 3
#> 5 2 115 2 1 0 5 2
#> 6 2 756 2 1 0 5 4
str(InstEval)
#> 'data.frame': 73421 obs. of 7 variables:
#> $ s : Factor w/ 2972 levels "1","2","3","4",..: 1 1 1 1 2 2 3 3 3 3 ...
#> $ d : Factor w/ 1128 levels "1","6","7","8",..: 525 560 832 1068 62 406 3 6 19 75 ...
#> $ studage: Ord.factor w/ 4 levels "2"<"4"<"6"<"8": 1 1 1 1 1 1 1 1 1 1 ...
#> $ lectage: Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 2 1 2 2 1 1 1 1 1 1 ...
#> $ service: Factor w/ 2 levels "0","1": 1 2 1 2 1 1 2 1 1 1 ...
#> $ dept : Factor w/ 14 levels "15","5","10",..: 14 5 14 12 2 2 13 3 3 3 ...
#> $ y : int 5 2 5 3 2 4 4 5 5 4 ...Starting with a simple model:
After fitting the model we can make use of the first function provided by merTools, fastdisp which modifies the function arm:::display to more quickly display a summary of the model without calculating the model sigma:
library(merTools)
fastdisp(m1)
#> lmer(formula = y ~ service + lectage + studage + (1 | d) + (1 |
#> s), data = InstEval)
#> coef.est coef.se
#> (Intercept) 3.22 0.02
#> service1 -0.07 0.01
#> lectage.L -0.19 0.02
#> lectage.Q 0.02 0.01
#> lectage.C -0.02 0.01
#> lectage^4 -0.02 0.01
#> lectage^5 -0.04 0.02
#> studage.L 0.10 0.02
#> studage.Q 0.01 0.02
#> studage.C 0.02 0.02
#>
#> Error terms:
#> Groups Name Std.Dev.
#> s (Intercept) 0.33
#> d (Intercept) 0.52
#> Residual 1.18
#> ---
#> number of obs: 73421, groups: s, 2972; d, 1128
#> AIC = 237655We see some interesting effects. First, our decision to include student and lecturer effects seems justified as there is substantial variance within these groups. Second, there do appear to be some effects by age and for lectures given as a service by an outside lecturer. Let’s look at these in more detail. One way to do this would be to plot the coefficients together in a line to see which deviate from 0 and in what direction. To get a confidence interval for our fixed effect coefficients we have a number of options that represent a tradeoff between coverage and computation time – see confint.merMod for details.
An alternative is to simulate values of the fixed effects from the posterior using the function arm::sim. Our next tool, FEsim, is a convenience wrapper to do this and provide an informative data frame of the results.
feEx <- FEsim(m1, 1000)
cbind(feEx[,1] , round(feEx[, 2:4], 3))
#> feEx[, 1] mean median sd
#> 1 (Intercept) 3.225 3.225 0.020
#> 2 service1 -0.070 -0.070 0.013
#> 3 lectage.L -0.186 -0.186 0.017
#> 4 lectage.Q 0.024 0.024 0.012
#> 5 lectage.C -0.025 -0.025 0.013
#> 6 lectage^4 -0.020 -0.019 0.014
#> 7 lectage^5 -0.039 -0.039 0.015
#> 8 studage.L 0.096 0.096 0.018
#> 9 studage.Q 0.005 0.005 0.017
#> 10 studage.C 0.017 0.017 0.016We can present these results graphically, using ggplot2: