Generalized linear mixed models in R: nitty-gritty

Licensed under the Creative Commons attribution-noncommercial license. Please share & remix noncommercially, mentioning its origin.

library(lme4)
library(broom)
library(broom.mixed)
library(dotwhisker)

Overview

Most aspects of GLMMs are carried over from LMMs (random effects) and GLMs (families and links).

integration methods

You do have to decide on an approximation method.

penalized quasi-likelihood
- fastest/least accurate
- uses quasi-likelihood (no AIC etc.?)
Laplace approximation
- fast, flexible compromise
adaptive Gauss-Hermite quadrature
- slowest
- most accurate
- most limited

In lme4, use the nAGQ= argument; nAGQ=1 (default) corresponds to Laplace approximation

integration rules

it only really matters in relatively extreme cases (small numbers of binary obs per group)
use the most accurate method that’s available and feasible

Laplace-approximation diagnostics

library(lattice)
aspect <- 0.6
xlab <- "z"; ylab <- "density"; type <- c("g","l"); scaled <- FALSE
mm <- readRDS("../../data/toenail_lapldiag.rds")
print(xyplot(y ~ zvals|id, data=mm,
             type=type, aspect=aspect,
             xlab=xlab,ylab=ylab,
             as.table=TRUE,
             panel=function(x,y,...){
    if (!scaled) {
        panel.lines(x, dnorm(x), lty=2)
    } else {
        panel.abline(h=1, lty=2)
    }
    panel.xyplot(x,y,...)
}))

comparing integration methods

g1 <- glmer(incidence/size ~ period + (1|herd),
            family=binomial,
            data=cbpp,
            weights=size)
g2 <- update(g1,nAGQ=5)
g3 <- update(g1,nAGQ=10)
g4 <- MASS:::glmmPQL(incidence/size ~ period,
                     random = ~1|herd,
                     data=cbpp,
                     family=binomial,
                     weights=size,
                     verbose=FALSE)
dwplot(list(Laplace=g1,AGQ5=g2,AGQ10=g3,glmmPQL=g4))

overdispersion

reminder: “too much” variance
only applies to families with estimated variance:
e.g. Poisson, binomial (\(N>1\))
sum^2 of Pearson residuals / residual degrees of freedom
e.g. aods3::gof()
observation-level random effects

aods3::gof(g1)

## D  = 73.4743, df = 51, P(>D) = 0.02133054
## X2 = 63.0674, df = 51, P(>X2) = 0.1196646

dealing with overdispersion

quasi-likelihood
observation-level random effects (Elston et al. 2001)
overdispersed distributions
(negative binomial, beta-binomial, etc.)
glmmTMB/brms (Brooks et al. 2017)

library(glmmTMB)
g5 <- glmmTMB(incidence/size ~ period + (1|herd),
              family=binomial,
              data=cbpp,
              weights=size)
g6 <- update(g5, family=betabinomial)

dotwhisker::dwplot(list(lme4=g1,glmmTMB=g5,glmmTMB_BB=g6),
                   effects="fixed")+
    geom_vline(xintercept=0,linetype=2)

other diagnostics

assumption of Normal residuals may not be very good
simulated residuals (DHARMa)

DHARMa::simulateResiduals(g1,plot=TRUE)

zero-inflation

glmmTMB etc.
Owls example

zipm3 <- glmmTMB(count~spp * mined + (1|site),
                ziformula=~spp * mined,
                Salamanders, family="poisson")
## can even include random effects in zero-inflation model!
zipm4 <- update(zipm3, ziformula = ~spp*mined+(1|site))
zipm5 <- update(zipm3, ziformula = ~(1|site))

library(bbmle)

## Loading required package: stats4

AICtab(zipm3,zipm4,zipm5)

##       dAIC df
## zipm4  0.0 30
## zipm3 25.9 29
## zipm5 28.1 17

Brooks, Mollie E., Kasper Kristensen, Koen J. van Benthem, Arni Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin Mächler, and Benjamin M. Bolker. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling.” R Journal 9 (2): 378–400. https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf.

Elston, D. A., R. Moss, T. Boulinier, C. Arrowsmith, and X. Lambin. 2001. “Analysis of Aggregation, a Worked Example: Numbers of Ticks on Red Grouse Chicks.” Parasitology 122 (5): 563–69.