• Understand the structural components of generalized linear mixed models

• Understand the options for fitting GLMMs and their pros and cons

• Understand some of the diagnostics available for evaluating GLMM fits

GLMMs combine the flexibility of non-normal distributions (GLMs) with the ability to address correlations among observations and nested data structures (LMMs)

Good news

• these extensions follow similar methods to GLMs and LMMs

Bad news

• these models are on the frontier of statistical research

• existing documentation is rather technical

• multiple approaches for fitting models; some with different results

Just like GLMs, GLMMs have three components:

1. Distribution of the data \(f(y; \theta)\)

2. Link function \(g(\eta)\)

3. Linear predictor \(\eta\)

We can write the linear predictor for GLMs as

\[ \eta = \beta_0 + \beta_1 x_1 + ... + \beta_k x_k \\ \Downarrow \\ \eta = \mathbf{X} \boldsymbol{\beta} \]

where the \(\beta_i\) are fixed effects of the covariates \(x_i\)

For GLMMs, our linear predictor also includes random effects

\[ \eta = \beta_0 + \beta_1 x_1 + \dots + \beta_k x_k + \alpha_0 + \alpha_1 z_1 + \dots + \alpha_l z_l\\ \Downarrow \\ \eta = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} \]

where the \(\beta_i\) are fixed effects of the covariates \(x_i\)

Survival of fish \(s_{i,j}\) as a function of length \(x_{i,j}\) in some location \(j\)

\[ \begin{aligned} \text{data distribution:} & ~~~ y_{i,j} \sim \text{Binomial}(N_{i,j}, s_{i,j}) \\ \text{link function:} & ~~~ \text{logit}(s_{i,j}) = \text{log}\left(\frac{s_{i,j}}{1-s_{i,j}}\right) = \mu_{i,j} \\ \text{linear model:} & ~~~ \mu_{i,j} = (\beta_0 + \alpha_j) + \beta_1x_{i,j} \\ & ~~~ \alpha_j \sim \text{N}(0, \sigma_{\delta}^2) \end{aligned} \]

Best practices suggest we try to keep things simple

Why? Because GLMMs involve solving an integral with no analytical solution

Recall that we think of likelihoods in terms of the observed data

But the random effects in our model are unobserved random variables, so we need to integrate them out of the likelihood

The likelihood for a GLMM involves integrating over all possible random effects

\[ \mathcal{L}(y; \boldsymbol{\beta}, \phi, \boldsymbol{\nu}) = \prod_{i} \int \underbrace{f_d(y; \boldsymbol{\beta}, \phi, \boldsymbol{\alpha})}_{\text{distn for data}} ~ \underbrace{f_r(\boldsymbol{\alpha}; \boldsymbol{\nu})}_{\text{distn for RE}} d \boldsymbol{\alpha} \]

If \(f(y; \boldsymbol{\beta}, \phi, \boldsymbol{\alpha})\) is not Gaussian, we cannot remove it from the likelihood, which makes it very difficult to compute

To avoid the integral, we will consider 3 methods that approximate the likelihood

They all have pros and cons so it’s not possible to pick the “best”

Penalized quasi-likelihood (PQL) uses a Taylor series expansion to approximate the linear predictor as an LMM

\[ \begin{align} g(\boldsymbol{\mu}) &= \boldsymbol{\eta} \\ &= \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} \\ &\Downarrow \\ g(\mathbf{y}) & \approx g(\boldsymbol{\mu}) + g'(\boldsymbol{\mu}) (\mathbf{y} - \boldsymbol{\mu}) \\ & \approx \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + g'(\boldsymbol{\mu}) \boldsymbol{\epsilon} \end{align} \]

The conditional variance of the data in a GLMM is then

\[ \begin{align} g(\mathbf{y}) &\approx \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + g'(\boldsymbol{\mu}) \boldsymbol{\epsilon} \\ &\Downarrow \\ g(\mathbf{y}) - \mathbf{X} \boldsymbol{\beta} &\approx \mathbf{Z} \boldsymbol{\alpha} + \boldsymbol{\epsilon} g'(\boldsymbol{\mu}) \\ &\Downarrow \\ \text{Var} \left( g(\mathbf{y}) - \mathbf{X} \boldsymbol{\beta} \right) &\approx \text{Var} \left( \mathbf{Z} \boldsymbol{\alpha} \right) + \text{Var} \left( \boldsymbol{\epsilon} g'(\boldsymbol{\mu}) \right) \end{align} \]

The conditional variance of the data in a GLMM is then

\[ \begin{align} g(\mathbf{y}) &\approx \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + g'(\boldsymbol{\mu}) \boldsymbol{\epsilon} \\ &\Downarrow \\ g(\mathbf{y}) - \mathbf{X} \boldsymbol{\beta} &\approx \mathbf{Z} \boldsymbol{\alpha} + \boldsymbol{\epsilon} g'(\boldsymbol{\mu}) \\ &\Downarrow \\ \text{Var} \left( g(\mathbf{y}) - \mathbf{X} \boldsymbol{\beta} \right) &\approx \text{Var} \left( \mathbf{Z} \boldsymbol{\alpha} \right) + \text{Var} \left( \boldsymbol{\epsilon} g'(\boldsymbol{\mu}) \right) \end{align} \]

which is similar to that for an LMM

\[ \text{Var} \left( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \right) = \text{Var} \left( \mathbf{Z} \boldsymbol{\alpha} \right) + \text{Var} \left( \boldsymbol{\epsilon} \right) \]

Pros

• fast, flexible, and widely implemented

Cons

• only asymptotically correct

• biased for Binomial and Poisson with small samples

• inference confounded by approximate likelihood

Laplace approximation is a long standing (1774) method for computing integrals of the form

\[ \int f(x) e^{\lambda g(x)} dx \]

This integrand is quite similar to the likelihood of a GLMM based on exponential distributions

Thus, we only need to find the maximum of \(g(x)\) and its second derivative, and evaluate them at only one point

Pros

• approximation of true likelihood rather than quasi-likelihood

• more accurate than PQL

Cons

• slower and less flexible than PQL

• may be impossible to compute for complex models

Gauss-Hermite quadrature is an expansion of Laplace approximation where the integrand is evaluated at more than one point

Quadrature is a method for numerically approximating an integral as a weighted sum

\[ \int f(u) e^{-u^2}du \approx \sum_i w_i f(u_i) \]

This method works by optimizing the placement and number of the \(u_i\) and the choice of the \(w_i\)

Pros

• More accurate than Laplace

Cons

• Slow and computationally intense

• Limited to a few random effects (one in practice)

Let’s consider a long-term study of invasive brown tree snakes in Guam

Introduced to the island shortly after WWII

Voracious predators on native birds and other vertebrates

Our data consist of counts of the number of eggs per female at 23 locations over 14 years

We are interested in the fixed effect of body size and the random effects of location and year

We’ll begin with only the effects of body size and location

We fit PQL models with MASS::glmPQL()

## load MASS
library(MASS)
## fit model
snakes_pql <- glmmPQL(eggs ~ size, random = ~1 | loc, data = df_eggs,
                      family = poisson)

summary(snakes_pql)

We can fit Laplace models with lme4::glmer()

## load lme4
library(lme4)
## fit model
snakes_lap <- glmer(eggs ~ size + (1 | loc), data = df_eggs,
                    family = poisson)

summary(snakes_lap)

We can fit GHQ models with lme4::glmer(..., nAGQ = pts)

## fit model
snakes_ghq <- glmer(eggs ~ size + (1 | loc), data = df_eggs,
                    family = poisson, nAGQ = 20)

summary(snakes_ghq)

Note: this method only works with one random effect

Here is a summary of the results from the 3 methods

##                 PQL    SE     Laplace   SE       GHQ    SE 
## (Intercept)    1.125 0.117      1.099 0.117     1.099 0.118
## size           0.508 0.079      0.506 0.066     0.506 0.067
## location SD    0.508    NA      0.525    NA     0.525    NA

What if we also want to include the random effect of year?

glmmPQL only allows for nested random effects

glmer(..., nAGQ = pts) only allows for one random effect

We can use the Laplace approximation via glmer

## fit model
snakes_lap_2 <- glmer(eggs ~ size + (1 | loc) + (1 | year),
                      data = df_eggs, family = poisson)

summary(snakes_lap_2)

Diagnostics for GLMMs are similar to those for GLMs, but we are limited in our choices

Recall our goodness of fit test based on the Pearson’s \(\chi^2\)

\[ X^2 = \sum_{i = 1}^n \frac{(O_i - E_i)^2}{E_i} \sim \chi^2_{(n - 1)} \]

where \(O_i\) is the observed count and \(E_i\) is the expected count

For a binomial distribution

\[ X^2 = \sum_{i = 1}^n \frac{(y_i - n_i \hat{p}_i)^2}{n_i \hat{p}_i} \]

For a Poisson distribution

\[ X^2 = \sum_{i = 1}^n \frac{(y_i - \lambda_i)^2}{\lambda_i} \]

\(H_0\): Our model is correctly specified

## Pearson's X^2 statistic
X2 <- sum((eggs - fitted(snakes_lap_2))^2 / fitted(snakes_lap_2))
## likelihood ratio test
pchisq(X2, df = nn - length(coef(snakes_lap_2)),
       lower.tail = FALSE)

## [1] 0.9986993

The \(p\)-value is large so we cannot reject \(H_0\)

For other models, we can calculate the leverages to evaluate potentially extreme values in predictor space

For GLMMs, however, the leverages depend on the estimated variance-covariance matrices of the random effects

For other models, we can calculate Cook’s distances to identify potentially influential data points

For GLMMs, however, the Cook’s distances involve derivatives of the likelihood with respect to the random effects (this is an active area of research)

As with GLMs, we can check GLMMs for evidence of overdispersion, which we estimate as

\[ \hat{c} = \frac{X^2}{n - k} \]

Let’s do so for our snake model applied to another data set

## Pearson's X^2 statistic
X2 <- sum((eggs - fitted(snakes_lap))^2 / fitted(snakes_lap))
## number of parameters
k <- length(coef(snakes_lap)) + length(ranef(snakes_lap))
## overdispersion parameter
(c_hat <- X2 / (nn - k))

## [1] 2.767328

pchisq(deviance(snakes_lap), k, lower.tail = FALSE)

## [1] 5.191758e-216

We can fit neg binomial models using Laplace approximation with lme4::glmer.nb()

## fit model
snakes_lap_nb <- glmer.nb(eggs ~ size + (1 | loc) + (1 | year),
                          data = df_eggs)

summary(snakes_lap_nb)

Adapted from Bolker et al (2009)