Where We Are

GLMMs Excel When

• Subject heterogeneity drives the association
• Covariate effects are interpreted conditional on latent subject traits
• Predicted subject trajectories or shrinkage estimates are needed
• Interest is at the individual level

GEE vs GLMM: Translation Guide

Key insight: Both answer valid but different scientific questions.

When to Choose Which

Why GLMMs Matter: Asthma Example

Why do they differ? Averaging the nonlinear logit over between-child variability shrinks the population-averaged OR toward 1 (non-collapsibility / Jensen). The subject-specific OR is larger. This is attenuation, not differential benefit.

Non-collapsibility (define it now)

Non-collapsibility: an odds ratio computed within strata can differ from the odds ratio computed in the pooled (marginal) population even when there is no confounding.

• Risk differences and (log) rate ratios are collapsible; odds ratios are not
• Mixing subjects with different baseline risk (\(\mathbf{b}_i\)) shifts the marginal OR toward 1
• This is a property of the odds-ratio scale, not a bias

Model Hierarchy

Level 1 (conditional on \(\mathbf{b}_i \in \mathbb{R}^q\)): \[Y_{ij} \mid \mathbf{b}_i \sim \text{EF}(\mu_{ij}, \phi), \quad g(\mu_{ij}) = \eta_{ij} = X_{ij}\beta + Z_{ij}\mathbf{b}_i\]

Level 2: \[\mathbf{b}_i \sim N(0, G)\]

• EF = exponential-family distribution (binomial, Poisson, …; recall L11)
• Grounding links: logit for binary \(\pi_{ij}\), log for Poisson rates

Jensen’s Inequality and Marginal Means

Marginal mean: \(E(Y_{ij}) \neq g^{-1}(X_{ij}\beta)\) because of Jensen’s inequality:

\[E\!\left[g^{-1}(\eta + b)\right] \neq g^{-1}\!\left(E[\eta + b]\right) \quad \text{for nonlinear } g^{-1}\]

Key: random effects create heterogeneity; averaging over it attenuates effects on the probability scale.

Jensen’s Inequality: Concrete Numerical Example

Setup: logistic model, random intercept \(b_i\) taking three equally likely values \(\{-2, 0, 2\}\).

At \(\eta = 0\):

Plug-in at \(b=0\): 0.500. Average of probabilities: \((0.119+0.500+0.881)/3 = 0.500\). They match here because \(\eta=0\) is symmetric.

Jensen’s Inequality: The Calculation

Now try \(\eta = 1\):

• Plug-in at \(b=0\): \(\text{expit}(1) = 0.731\)
• Average the probabilities: \((0.269+0.731+0.953)/3 = 0.651\)

Difference: 8 percentage points. This is why re.form=NA does NOT give marginal predictions.

Likelihood Challenge

Subject-specific contribution: \[L_i(\beta, G) = \int_{\mathbb{R}^q} \prod_{j=1}^{n_i} f(y_{ij} \mid \mathbf{b}_i;\, \beta)\; f_G(\mathbf{b}_i)\; d\mathbf{b}_i\]

• \(f_G(\mathbf{b}_i) = N(\mathbf{b}_i;\, 0, G)\) is the random-effects density
• No closed form when \(f(\cdot\mid\mathbf{b}_i)\) is non-Gaussian
• This integral is what makes GLMM fitting hard

Estimation: ML via Numerical Integration

Random intercept only, moderate n_i (>=5)?
  -> Laplace approximation (lme4 default, nAGQ = 1): fast, usually adequate
Random slopes, or small/sparse clusters?
  -> Adaptive Gauss-Hermite quadrature (nAGQ = 5-50): more accurate
Complex G, sparse data, full uncertainty wanted?
  -> Bayesian MCMC (brms, MCMCglmm)

R Setup

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(lme4)
  library(emmeans)
  library(broom.mixed)
  library(patchwork)
})
theme_set(theme_minimal(base_size = 16))

Example 1: Binary Asthma Control Trial (simulated)

• 120 children, quarterly visits for 2 years
• Outcome: controlled asthma (1) vs uncontrolled (0)
• Predictors: treatment (steroid vs usual care), visit time, baseline severity
• Random intercept captures latent propensity for control

Example 1: Data Simulation

set.seed(667)
n_subj <- 120; visits <- 8
asthma_dat <- tibble(
  id    = rep(1:n_subj, each = visits),
  visit = rep(0:(visits - 1), times = n_subj),
  weeks = rep(0:(visits - 1), times = n_subj) * 12,
  treat = factor(rep(sample(c("Steroid","Usual"), n_subj, replace = TRUE), each = visits)),
  severity = rep(sample(c("Mild","Moderate","Severe"), n_subj, replace = TRUE,
                        prob = c(0.5, 0.35, 0.15)), each = visits)
) %>%
  mutate(
    sev_score = case_when(severity == "Mild" ~ 0, severity == "Moderate" ~ 0.6, TRUE ~ 1.2),
    treat_num = if_else(treat == "Steroid", 0.6, 0),   # TRUTH: steroid log-OR = 0.6
    time_eff  = 0.05 * weeks / 12,                      # TRUTH: time log-OR per yr = 0.05
    rand_int  = rep(rnorm(n_subj, sd = 0.9), each = visits),  # TRUTH: RE intercept SD = 0.9
    linpred   = -1 + treat_num + time_eff - 0.8 * sev_score + # TRUTH: severity log-OR = -0.8
                0.15 * treat_num * (weeks / 12) + rand_int,    # TRUTH: treat x time = 0.15
    prob      = plogis(linpred),
    control   = rbinom(n_subj * visits, size = 1, prob = prob)
  )
asthma_dat$treat <- relevel(asthma_dat$treat, ref = "Usual")

Fit the Random-Intercept Logistic GLMM

glmm_fit <- glmer(
  control ~ treat * weeks + severity + (1 | id),
  data = asthma_dat,
  family = binomial(link = "logit")
)
broom.mixed::tidy(glmm_fit, effects = "fixed", conf.int = TRUE) %>%
  dplyr::select(term, estimate, std.error, conf.low, conf.high)

Interpreting the Fixed Effects (conditional)

• \(\exp(\beta)\) is an odds ratio conditional on the random intercept (within a given child)
• Positive steroid effect = better control for the same child
• The interaction tests diverging within-child trajectories
• Critical: these are NOT population-averaged effects

Conditional vs Marginal: Why They Differ

Converting Conditional to Marginal (Zeger approx.)

For random-intercept logistic models (Zeger, Liang & Albert 1988; FLW Ch. 16 / p. 477):

\[\beta^{\text{marg}} \approx \frac{\beta^{\text{cond}}}{\sqrt{1 + c_2^{2}\,\sigma_b^2}}, \qquad c_2^{2} \approx 0.346\]

Applying the Zeger Approximation

beta_cond <- fixef(glmm_fit)["treatSteroid"]
var_u     <- as.numeric(attr(VarCorr(glmm_fit)$id, "stddev"))^2
atten     <- sqrt(1 + 0.346 * var_u)      # 0.346 = c_2^2 multiplier
beta_marg <- beta_cond / atten

cat("Conditional beta (GLMM):", round(beta_cond, 3),
    " -> OR =", round(exp(beta_cond), 2), "\n")

Side-by-Side: Conditional vs Marginal

Check Your Understanding: Zeger

A GLMM for smoking cessation gives conditional log-OR \(\beta^{\text{cond}} = 1.1\) (OR = 3.0) and \(\sigma_b^2 = 2.0\). Approximate marginal OR?

Recap of Part I

• GLMM: \(g(E[Y_{ij}\mid\mathbf{b}_i]) = X_{ij}\beta + Z_{ij}\mathbf{b}_i\), \(\mathbf{b}_i\sim N(0,G)\)
• \(\exp(\beta)\) is a conditional (within-subject) odds/rate ratio
• Marginal effect is attenuated by Jensen / non-collapsibility (Zeger approx., Ch. 16)
• Fit by ML (Laplace / adaptive quadrature); PQL/MQL are Ch. 15
• GLMM valid under MAR; standard GEE needs MCAR

How to Report GLMM Results

Correct:

• “Within a given child, steroid increased odds of control by OR = 2.73 (95% CI …)”
• “At the population level (marginal), the OR is approximately 2.36”

Incorrect:

• “Steroid increased the odds of control” (conditional or marginal?)

Always state the frame of reference.

Three Types of Predictions

For nonlinear links, “at \(\mathbf{b}=0\)” \(\neq\) “marginal mean.”

Common Misconception: re.form = NA

# WRONG to call this a marginal prediction.
# It is E(Y | b = 0): conditional at the median latent subject.
pred_at_b0 <- predict(glmm_fit,
  newdata = data.frame(treat = "Steroid", weeks = 48, severity = "Mild"),
  re.form = NA, type = "response")
cat("Prediction at b=0:", round(pred_at_b0, 4), "  (this is E(Y|b=0), NOT E(Y))\n")

True Marginal Predictions via Monte Carlo

sigma_b   <- as.numeric(attr(VarCorr(glmm_fit)$id, "stddev"))
nd        <- data.frame(treat = "Steroid", weeks = 48, severity = "Mild")
eta_fixed <- predict(glmm_fit, newdata = nd, re.form = NA, type = "link")

set.seed(667)
b_sim       <- rnorm(5000, mean = 0, sd = sigma_b)   # draw from N(0, G)
marginal_pred <- mean(plogis(eta_fixed + b_sim))      # average expit over b

cat("At b=0 (NOT marginal):           ", round(pred_at_b0, 4), "\n")

Why the Difference Matters

A Real FLW Case Study: Seizure Counts

• Outcome: seizure count per period, Poisson with offset \(\log(T_{ij})\) (\(T = 8\) wk baseline, \(2\) wk follow-up)
• Random intercept and slope on Time (baseline vs follow-up): \(q = 2\)
• Question: does progabide reduce the subject-specific seizure rate?

Load and Reshape the Seizure Data

sz_url  <- "https://content.sph.harvard.edu/fitzmaur/ala2e/epilepsy.txt"
sz_file <- "../../data/epilepsy-data.txt"     # local copy (header stripped)
if (!file.exists(sz_file)) {
  raw <- readLines(url(sz_url))               # FLW file has a comment header
  dat <- raw[grepl("^\\s*[0-9]", raw)]        # keep numeric rows only
  writeLines(dat, sz_file)
}
sz_wide <- read.table(sz_file, header = FALSE,
  col.names = c("id", "trt", "age", "baseline", "y1", "y2", "y3", "y4"))

sz_long <- sz_wide %>%
  pivot_longer(c(baseline, y1, y2, y3, y4),
               names_to = "period_lab", values_to = "count") %>%
  mutate(
    period = recode(period_lab, baseline = 0L, y1 = 1L, y2 = 2L, y3 = 3L, y4 = 4L),
    Time   = if_else(period == 0L, 0L, 1L),         # 0 = baseline, 1 = follow-up
    Tij    = if_else(period == 0L, 8, 2),           # exposure weeks
    logT   = log(Tij),
    Trt    = trt                                    # 0 = placebo, 1 = progabide
  ) %>% arrange(id, period)
cat("Rows:", nrow(sz_long), " Subjects:", dplyr::n_distinct(sz_long$id), "\n")

Fit the Seizure GLMM (FLW Table 14.6)

sz_fit <- glmer(
  count ~ Trt * Time + (1 + Time | id),
  offset = logT, data = sz_long,
  family = poisson(link = "log"),
  control = glmerControl(optimizer = "bobyqa")
)
broom.mixed::tidy(sz_fit, effects = "fixed", conf.int = TRUE) %>%
  dplyr::select(term, estimate, std.error, conf.low, conf.high)

Interpreting the Seizure Fit (subject-specific)

• All covariates are dichotomous, so \(\exp(\beta)\) are subject-specific (log) rate ratios (FLW Table 14.5)
• \(\exp(\beta_4)\) (Trt\(\times\)Time) is a ratio of rate ratios: progabide vs placebo change from baseline
• \(\widehat{\beta}_4 = -0.306 < 0\): progabide reduces the seizure rate from baseline
• For a “typical” progabide patient (\(b_{2i}=0\)): rate ratio \(e^{\widehat\beta_2 + \widehat\beta_4} \approx 0.74\) (about a 26% drop)

Seizure Variance Components

print(VarCorr(sz_fit), comp = c("Variance", "Std.Dev."))

The Same GLMM in SAS (reference)

Static reference (not run here): the seizure log-linear GLMM in PROC GLIMMIX, mirroring FLW Table 14.11 (ML via adaptive quadrature).

proc glimmix data=seizure method=quad(qpoints=50);
   class id trt;
   model count = trt time trt*time / dist=poisson link=log solution
                 offset=logT;
   random intercept time / subject=id type=un;
run;

Overdispersion Check

pr  <- residuals(sz_fit, type = "pearson")
od  <- sum(pr^2) / df.residual(sz_fit)
cat("Overdispersion ratio (Pearson X^2 / df):", round(od, 2), "\n")

Overdispersion Remedies

Option 1 - observation-level random effect (OLRE):

sz_long$obs <- factor(seq_len(nrow(sz_long)))
glmer(count ~ Trt * Time + (1 + Time | id) + (1 | obs),
      offset = logT, data = sz_long, family = poisson)

Option 2 - negative binomial (FLW Table 14.8):

glmmTMB::glmmTMB(count ~ Trt * Time + (1 + Time | id),
                 offset = logT, data = sz_long, family = nbinom2)

Hypothesis Testing: Fixed Effects

• Wald tests (from summary()) are asymptotic; small clusters \(\rightarrow\) prefer profile likelihood

confint(sz_fit, parm = "beta_", method = "profile")  # profile CIs for fixed effects

Hypothesis Testing: Variance Components

Testing \(H_0\): random slope dropped ((1+Time|id) \(\to\) (1|id)) removes two parameters (slope variance and intercept-slope covariance).

Parametric-Bootstrap LRT for the Random Slope

null_fit <- glmer(count ~ Trt * Time + (1 | id), offset = logT,
                  data = sz_long, family = poisson,
                  control = glmerControl(optimizer = "bobyqa"))
lr_obs <- as.numeric(2 * (logLik(sz_fit) - logLik(null_fit)))

set.seed(667)
B <- 100
sim_y   <- simulate(null_fit, nsim = B)          # simulate under the null
lr_boot <- vapply(seq_len(B), function(b) {
  f0 <- refit(null_fit, newresp = sim_y[[b]])
  f1 <- refit(sz_fit,  newresp = sim_y[[b]])
  as.numeric(2 * (logLik(f1) - logLik(f0)))
}, numeric(1))

p_boot <- (1 + sum(lr_boot >= lr_obs)) / (1 + B)
p_mix  <- 0.5 * pchisq(lr_obs, 1, lower.tail = FALSE) +
          0.5 * pchisq(lr_obs, 2, lower.tail = FALSE)
cat("LR statistic:        ", round(lr_obs, 2), "\n")

Post-Hoc Contrasts with emmeans

emm_contrasts <- emmeans(glmm_fit, specs = pairwise ~ treat | weeks,
                         at = list(weeks = c(0, 48)), type = "response")
summary(emm_contrasts$contrasts)

Interpreting emmeans Contrasts

These are subject-specific odds ratios evaluated at \(\mathbf{b}_i = 0\) (the median subject on the logit scale), NOT the population-average OR.

• At \(\mathbf{b}_i = 0\) = conditional on the median latent subject
• The population-average OR is attenuated (Jensen); see the Three Types of Predictions table
• Do not read these as “the average subject in the population”

Random-Effects Heterogeneity (seizure slopes)

Diagnostics: Residuals vs Fitted

Diagnostics: Random-Effects Q-Q

Convergence Troubleshooting Recipe

Rule of thumb: if a random-slope variance is \(< 0.01\) and convergence struggles, drop it.

Convergence: A Worked Fix

# Rescaling the continuous predictor often resolves convergence trouble.
asthma_dat$weeks_z <- as.numeric(scale(asthma_dat$weeks))
fit_scaled <- glmer(control ~ treat * weeks_z + severity + (weeks_z | id),
                    data = asthma_dat, family = binomial,
                    control = glmerControl(optimizer = "bobyqa"))
cat("Converged:", is.null(fit_scaled@optinfo$conv$lme4$messages), "\n")

Common Pitfalls and Remedies

Reporting Checklist

• State “conditional on the random effects” (or “within-subject”) when interpreting \(\beta\)
• Report the estimand: conditional OR/RR, and the marginal version if relevant
• Report variance components with CIs
• State the estimation method (Laplace, nAGQ = X, or MCMC)
• For variance-component tests, state the mixture / bootstrap used

Extensions and Software

Wrap-Up

• GLMM = GLM plus latent subject effects; inference is subject-specific (conditional)
• The marginal effect is attenuated (Jensen / non-collapsibility); Zeger approximates it (Ch. 16)
• Fit by ML (Laplace / adaptive quadrature); valid under MAR
• re.form=NA is conditional at \(\mathbf{b}=0\), not marginal; integrate for true marginal
• Test variance components with the right mixture or a parametric bootstrap (never exactRLRT for a GLMM)
• Next chapter (Ch. 15): approximate likelihood methods (PQL/MQL) and their tradeoffs

Reading

Required: Fitzmaurice et al. Ch. 14, sections 14.1-14.6, plus the Appendix on adaptive Gaussian quadrature.

Optional:

• Zeger, Liang & Albert (1988), marginal vs subject-specific approximation (the 0.346 factor; see FLW Ch. 16)
• Bolker et al. (2009), GLMMs: a practical guide

Session Info

sessionInfo()