• Understand some of the diagnostics available for evaluating GLMM fits

• Become familiar with some of the inferential tools for GLMMs

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

Our data consist of counts of the number of eggs per female at 10 locations over 13 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

library(lme4)

## set control parameters for model fitting
cntrl <- glmerControl(optimizer = "bobyqa", tol = 1e-4, optCtrl = list(maxfun = 100000))

## fixed effect of size plus random effect of location
fit_1 <- glmer(eggs ~ size + (1|loc),
               family = poisson,
               data = df_eggs,
               control = cntrl)

## fixed effect of size plus random effects of location & year
fit_2 <- glmer(eggs ~ size + (1|loc) + (1|year),
               family = poisson,
               data = df_eggs, 
               control = cntrl)

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(fit_2))^2 / fitted(fit_2))

## likelihood ratio test
pchisq(X2, df = nn - length(coef(fit_2)),
       lower.tail = FALSE)

## [1] 0.6447966

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

Check for heteroscedasticity

Use DHARMa::simulateResiduals to scale the residuals

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

We can fit a GLM and evaluate leverage of the fixed 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)

Again, we can fit a GLM and look at the Cook’s distances for the fixed effects

Large snakes with lots of eggs

## large Cook's D
cooks_big <- sort(cooks_d, decreasing = TRUE)[1:2]       

## details about them
df_eggs[as.numeric(names(cooks_big)),]

##     eggs size loc year
## 138   30 0.92   1   12
## 186   28 0.89   1    8

For confidence intervals, we will again use bootstrapping

• We will cover this in lab

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))^2 / fitted(snakes))

## number of parameters
k <- length(coef(snakes)) + length(ranef(snakes))

## overdispersion parameter
(c_hat <- X2 / (nn - k))

## [1] 2.180794

## test for overdispersion
pchisq(deviance(snakes), k, lower.tail = FALSE)

## [1] 0

## neg binomial model for overdispersion
fit_nb <- glmer.nb(eggs_nb ~ size + (1 | loc) + (1 | year), 
                   data = df_eggs,
                   control = cntrl)

• Diagnostics for GLMMs are an art rather than a science

• Use LRTs to test random effects (or just leave them in)

• Use robust inferential methods like bootstrapped confidence intervals