Motivation

When data are collected in clusters, or multiple measurements are made on the same subjects, the assumption of independent observations no longer holds. It is likely the clusters of data, or repeated measures, are correlated or related in some way. In addition, when our dependent variable is constrained to only two discrete values (eg, Success and Failure), we cannot rely on Normally distributed errors. A common method for addressing these issues is to use Generalized Linear Mixed Models (GLMM).

Simulate data

Understanding a GLMM analysis can be facilitated by simulating data for a GLMM. Let’s pretend we want data with the following features:

We want to simulate data such that y varies depending on treat, wl, and period. Let’s simulate data for 100 subjects that have 34 x 3 measurements:

n <- 100
id <- seq(n)
treat <- c("G", "M", "S")
wl <- c("L", "T", "H")
period <- 1:5
trial <- 1:3
d <- expand.grid(period = period, treat = treat, wl = wl, 
                 trial = trial, id = id)
# treatment "M" has 3 periods of measure instead of 5
d <- d[!(d$treat == "M" & d$period %in% c(4,5)), ]
# treatment "S" has 2 workloads instead of 3
d <- d[!(d$treat == "S" & d$wl=="T"), ]
# arrange for inspection
d <- d[order(d$id, d$treat, d$wl, d$period, d$trial), 
       c("id", "trial", "treat", "wl", "period")]

Now we simulate the response using a logistic regression model:

\[\pi_{i} = \frac{\text{exp}(\beta_{0} + \beta_1x_{1} + \cdots + \beta_px_{p} + b_i)}{1 + \text{exp}(\beta_0 + \beta_1x_{1} + \cdots + \beta_px_{p} + b_i)}\]

where \(b_i\) is the random effect of subject i. \(b_i\) is assumed to be drawn from a normal distribution with mean 0 and some constant standard deviation. The estimated random effect is the estimate of that constant standard deviation.

The formula above is available as the plogis function in R.

Below we generate interactions for period:treat and period:wl.

set.seed(111)
# generate random effect for each subject
b <- rnorm(n, mean = 0, sd = 0.3)

# generate linear predictor
m <- 0.1 + 
  0.3*(d$period) + 
  0.2*(d$treat=="M") + 0.3*(d$treat=="S") +
  0.2*(d$wl=="T")+ 0.3*(d$wl=="H") +
  0.2*d$period*(d$treat=="M") + -0.2*d$period*(d$treat=="S") +
  0.1*d$period*(d$wl=="T") + -0.3*d$period*(d$wl=="H") +
  b[d$id]

# generate probability using logistic function
p <- plogis(m)

# generate 0/1 using random binomial distribution
d$y <- rbinom(nrow(d), size = 1, prob = p)

Fit model

Now we’re ready to analyze the data and see if we can recover the TRUE values we used to simulate the data. We have to use lme4 for this; there is no generalized mixed-effect methods in the nlme package.

Let’s fit the “correct” model. (Note: I set nAGQ = 0, which results in a faster but less exact form of parameter estimation. This is mainly to speed up the creation of this document.)

library(lme4)
# correct model; formula matches formula used above
m <- glmer(y ~ period + treat + wl + period:treat + period:wl + (1|id),
           data = d, family = binomial, nAGQ = 0)
summary(m, corr = FALSE)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y ~ period + treat + wl + period:treat + period:wl + (1 | id)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##  11842.2  11921.7  -5910.1  11820.2    10189 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8988 -1.0541  0.4929  0.7022  1.6294 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.07543  0.2746  
## Number of obs: 10200, groups:  id, 100
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.02250    0.10596  -0.212  0.83182    
## period         0.35294    0.03424  10.307  < 2e-16 ***
## treatM         0.21035    0.14019   1.500  0.13349    
## treatS         0.38300    0.12291   3.116  0.00183 ** 
## wlT            0.25083    0.14338   1.749  0.08022 .  
## wlH            0.49900    0.10933   4.564 5.01e-06 ***
## period:treatM  0.19173    0.06343   3.023  0.00250 ** 
## period:treatS -0.23598    0.03845  -6.137 8.44e-10 ***
## period:wlT     0.06626    0.05391   1.229  0.21908    
## period:wlH    -0.39611    0.03617 -10.950  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice the estimated coefficients are pretty close to the “true” values. For example, the coefficient forperiod:treatS is estimated as -0.2359774, which is pretty close to -0.2, the value we used above (ie, -0.2*d$period*(d$treat=="S")). Likewise the estimated random effect Std.Dev. is estimated as 0.2746448, which is pretty close to 0.3, the true value used in the simulation (ie, b <- rnorm(n, mean = 0, sd = 0.3)).

With a fitted model, we can create effect plots.

library(effects)
plot(Effect(c("period", "treat"), mod = m), 
     rug = FALSE, multiline = TRUE, ci.style = "bands", 
     type="response")

The probability of success for the S group never really deviates from about 0.65. The probabilities of success for the M and G groups increase over time, with M having a higher probability of success.

Notice the M group is getting predictions for periods 4 and 5 for which we have no data. We can save the effect predictions into a data frame and subset out the predictions for period 4 and 5 for the M group and re-draw the plot. I find ggplot easier for this task.

eff_df <- as.data.frame(Effect(c("period", "treat"), mod = m, ))
eff_df <- eff_df[!(eff_df$treat == "M" & eff_df$period %in% c(4,5)),]
library(ggplot2)
ggplot(eff_df) +
  aes(x = period, y = fit, color = treat, fill = treat) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 1/4)

We could also remove group M altogether using the same strategy above.

Let’s fit a wrong model. We fit no interactions and specify a random effect for period in addition to the intercept.

m2 <- glmer(y ~ period + treat + wl + (period|id),
           data = d, family = binomial, nAGQ = 0)
summary(m2, corr = FALSE)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y ~ period + treat + wl + (period | id)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##  12092.6  12157.7  -6037.3  12074.6    10191 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7991 -1.0903  0.5142  0.6897  1.5117 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  id     (Intercept) 0.126147 0.35517       
##         period      0.001071 0.03273  -1.00
## Number of obs: 10200, groups:  id, 100
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.67681    0.07619   8.883  < 2e-16 ***
## period       0.10035    0.01723   5.824 5.76e-09 ***
## treatM       0.48836    0.05981   8.166 3.20e-16 ***
## treatS      -0.32552    0.05311  -6.129 8.82e-10 ***
## wlT          0.34735    0.06690   5.192 2.08e-07 ***
## wlH         -0.58059    0.04882 -11.892  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can compare the models via AIC since one is not nested within the other. The first model is appreciably better with a lower AIC.

AIC(m, m2)
##    df      AIC
## m  11 11842.17
## m2  9 12092.59

If we want we can simulate a random effect for period. Notice we draw 100 random normal values from N(0, 0.4) for each subject, multiply by period, and to the model.

set.seed(222)
# generate random effect for each subject
b0 <- rnorm(n, mean = 0, sd = 0.3)
b1 <- rnorm(n, mean = 0, sd = 0.4)

# generate linear predictor
m <- 0.1 + 
  0.3*d$period + 
  0.2*(d$treat=="M") + 
  0.3*(d$treat=="S") +
  0.2*(d$wl=="T")+ 
  0.3*(d$wl=="H") +
  0.2*d$period*(d$treat=="M") + 
  -0.2*d$period*(d$treat=="S") +
  0.1*d$period*(d$wl=="T") + 
  -0.3*d$period*(d$wl=="H") +
  b0[d$id] + b1[d$id]*d$period # random effects

# generate probability using logistic function
p <- plogis(m)

# generate 0/1 using random binomial distribution
d$y <- rbinom(nrow(d), size = 1, prob = p)

Now we fit the correct model. Notice the double bars in the random effect formula syntax (||). That says do not estimate the correlation between the random effects for period and the intercept. That is correct because they are not correlated. We did not specify that in the simulation.

# correct model
m3 <- glmer(y ~ period + treat + wl + period:treat + period:wl + 
              (period||id),
            data = d, family = binomial, nAGQ = 0)

Let’s fit the “wrong” model that allows the random effects to be correlated. (ie, just one bar in the random effect syntax)

# wrong model
m4 <- glmer(y ~ period + treat + wl + period:treat + period:wl + 
              (period|id),
            data = d, family = binomial, nAGQ = 0)
VarCorr(m4)
##  Groups Name        Std.Dev. Corr  
##  id     (Intercept) 0.42159        
##         period      0.39078  -0.030

The estimated correlation is close to 0. Comparing AIC values is not super decisive. The larger model with correlation has a slightly higher AIC.

AIC(m3, m4)
##    df      AIC
## m3 12 10906.72
## m4 13 10908.69

In real life we never know the correct model. The above was meant to demonstrate the assumptions we make when we fit a model.

One-way ANOVA style model

We can combine the treat and wl groups into one factor using the interaction function. I call this variable WL.

# droplevels drops "S.T" as a level since there is no Sudden-Transition grouping.
d$WL <- droplevels(interaction(d$treat, d$wl))

Let’s simulate new data. Notice we don’t interact period with WL. This is mainly to help me simulate data. An interaction would have resulted in a much longer formula!

set.seed(333)
# generate random effects for each subject
b0 <- rnorm(n, mean = 0, sd = 0.3)
b1 <- rnorm(n, mean = 0, sd = 0.2)

# generate linear predictor
m <- 0.1 + 
  0.3*d$period + 
  0.1*(d$WL == "M.L") +
  0.2*(d$WL == "S.L") +
  0.3*(d$WL == "G.T") +
  0.2*(d$WL == "M.T") +
  -0.1*(d$WL == "G.H") +
  -0.2*(d$WL == "M.H") +
  -0.3*(d$WL == "S.H") +
  b0[d$id] + b1[d$id]*d$period # random effects

# generate probability using logistic function
p <- plogis(m)

# generate 0/1 using random binomial distribution
d$y <- rbinom(nrow(d), size = 1, prob = p)

Now fit the “correct” model.

m5 <- glmer(y ~ period + WL + (period||id), data = d, 
            family = binomial, nAGQ = 0)
summary(m5, corr = FALSE)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y ~ period + WL + (period || id)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##  11506.9  11586.5  -5742.5  11484.9    10189 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6853 -1.0018  0.4704  0.6668  1.8751 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.09172  0.3029  
##  id.1   period      0.03969  0.1992  
## Number of obs: 10200, groups:  id, 100
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.23448    0.08447   2.776 0.005505 ** 
## period       0.28918    0.02740  10.556  < 2e-16 ***
## WLM.L       -0.01745    0.09638  -0.181 0.856304    
## WLS.L        0.16381    0.08736   1.875 0.060763 .  
## WLG.T        0.30622    0.08887   3.446 0.000569 ***
## WLM.T        0.12027    0.09777   1.230 0.218666    
## WLG.H       -0.13671    0.08487  -1.611 0.107205    
## WLM.H       -0.38075    0.09386  -4.057 4.98e-05 ***
## WLS.H       -0.35416    0.08360  -4.236 2.27e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing for differences between levels of WL can be carried out using the multcomp package. Setting up a custom contrast takes extra work. Below is an example.

K is the contrast matrix. Each row represents a comparison. “G.L” is the reference level, and R fits treatment contrasts by default, so the coefficients for all other levels are difference from “G.L”. So any comparison with “G.L” is simply that coefficient (1st and 4th rows of K). The number of zeroes and ones (in this case, 9) corresponds to the number of coefficients in our model.

library(multcomp)

# create the contrast for multiple comparisons
K <- rbind("S.L - G.L" = c(0, 0, 0, 1, 0, 0, 0, 0, 0),
           "S.H - G.H" = c(0, 0, 0, 0, 0, 0, -1, 0, 1),
           "M.H - G.H" = c(0, 0, 0, 0, 0, 0, -1, 1, 0),
           "M.L - G.L" = c(0, 0, 1, 0, 0, 0, 0, 0, 0))
mc <- glht(m5, linfct = K)
summary(mc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glmer(formula = y ~ period + WL + (period || id), data = d, family = binomial, 
##     nAGQ = 0)
## 
## Linear Hypotheses:
##                Estimate Std. Error z value Pr(>|z|)  
## S.L - G.L == 0  0.16381    0.08736   1.875   0.2098  
## S.H - G.H == 0 -0.21745    0.08253  -2.635   0.0322 *
## M.H - G.H == 0 -0.24404    0.09296  -2.625   0.0331 *
## M.L - G.L == 0 -0.01745    0.09638  -0.181   0.9995  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

The tests have adjusted p-values to help protect against false positives. The default is the single-step method. There are many others to choose from. The most conservative is the bonferroni correction. Here’s how to specify it, although the results are not much different and thus not shown.

summary(mc, test = adjusted(type = "bonferroni"))

To see confidence intervals on the estimated differences use the confint function.

confint(mc)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: glmer(formula = y ~ period + WL + (period || id), data = d, family = binomial, 
##     nAGQ = 0)
## 
## Quantile = 2.4773
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                Estimate lwr      upr     
## S.L - G.L == 0  0.16381 -0.05259  0.38022
## S.H - G.H == 0 -0.21745 -0.42190 -0.01300
## M.H - G.H == 0 -0.24404 -0.47434 -0.01374
## M.L - G.L == 0 -0.01745 -0.25622  0.22132

The multcomp package has some basic plotting methods for these tests as well. The difference between levels is plotted along with 95% error bars. Overlap with the dotted 0 line indicates uncertainty about whether or not the magnitude is positive or negative.

par(mar = c(5.1, 5.1, 4.1, 2.1))
plot(mc)

Basic diagnostics

Recall the assumption: random effects are drawn from a Normal distribution with constant variance. Two plots to help assess this assumption are dotplot and qqmath. We call these plots on the random effects which we extract with ranef.

The dotplot is good for checking the constant variance assumption. This plot is sometimes called a “caterpillar plot”. We like what we see here. The error bars around the random effects are a constant size. Also we don’t see any random effects too far above or below 0.

re <- ranef(m5)
lattice::dotplot(re)
## $id

To check the normality assumption, call qqmath on each random effect. Again these look good.

lattice::qqmath(re$id$`(Intercept)`)

lattice::qqmath(re$id$period)