16.1 The Main Idea

• Marginal and mixed models differ beyond correlation handling; they answer different questions.

• Linear models mask this difference because \(\operatorname{E}(Y_i)=X_i\beta\) either way.

• Once the link is non-linear, \(\beta\) (population) and \(\beta^*\) (subject) diverge.

• Always start by writing down the estimand before picking software.

16.2 Linear Special Case: Why Equivalence Holds

With identity link, both GLS and LME imply \(\operatorname{E}(Y_i)=X_i\beta\).

Why? Averaging a linear function of \(b_i\) is straightforward: \[\operatorname{E}[Z_i b_i] = Z_i \operatorname{E}(b_i) = 0\]

Random effects only restructure the covariance: \[\boldsymbol\Sigma_i = \mathbf{Z}_i G \mathbf{Z}_i^\top + \sigma^2 \mathbf{I}\]

where \(G\) is the between-subject random-effects covariance and \(\sigma^2 \mathbf{I}\) is the within-subject measurement (error) covariance.

Therefore \(\boldsymbol\beta\) retains a marginal interpretation in linear mixed models: here \(\boldsymbol\beta = \boldsymbol\beta^*\).

16.2 Quick Derivation

Start from the LME conditional mean: \[\operatorname{E}(\mathbf{Y}_i \mid \mathbf{b}_i) = \mathbf{X}_i\boldsymbol\beta + \mathbf{Z}_i \mathbf{b}_i\]

Apply law of total expectation: \[\operatorname{E}(\mathbf{Y}_i) = \operatorname{E}\big[\operatorname{E}(\mathbf{Y}_i \mid \mathbf{b}_i)\big]\]

Substitute: \[\operatorname{E}(\mathbf{Y}_i) = \mathbf{X}_i\boldsymbol\beta + \mathbf{Z}_i \operatorname{E}(\mathbf{b}_i) = \mathbf{X}_i\boldsymbol\beta\]

when \(\operatorname{E}(\mathbf{b}_i) = \mathbf{0}\).

Why Non-Linearity Breaks Equivalence

Mathematical reason: Expectation is a linear operator.

where \(h = g^{-1}\) is the inverse link (for logistic, \(h = \operatorname{expit}\)).

Visual Intuition

Plain English: Averaging curved lines gives a different curve than the average line.

16.3 Link Function Consequences

Define marginal model: \(g(\mu_i) = \mathbf{X}_i \boldsymbol\beta\) with \(\mu_i = \operatorname{E}(Y_i)\).

Define GLMM: \(g\{\operatorname{E}(Y_i\mid \mathbf{b}_i)\} = \mathbf{X}_i \boldsymbol\beta^* + \mathbf{Z}_i \mathbf{b}_i\).

Marginal mean requires integrating the inverse link \(h\) over \(\mathbf{b}_i\): \[\mu_i = \int h(\mathbf{X}_i\boldsymbol\beta^* + \mathbf{Z}_i \mathbf{b}_i)\, f(\mathbf{b}_i)\,d\mathbf{b}_i\]

Result: \(g(\mu_i) \neq \mathbf{X}_i \boldsymbol\beta^*\) in general; the link form is not preserved.

16.3 Attenuation Insight

For a logistic GLMM with a random intercept: \[\text{logit}\{\operatorname{E}(Y_{ij}\mid b_i)\} = X_{ij}\boldsymbol\beta^* + b_i, \qquad b_i \sim N(0,\sigma_b^2)\]

Approximate marginal relationship (Zeger, Liang & Albert, 1988; FLW p. 477): \[\text{logit}\{\operatorname{E}(Y_{ij})\} \approx c \, X_{ij}\boldsymbol\beta^*, \qquad c = \left(1 + c_2^{\,2}\,\sigma_b^2\right)^{-1/2}\]

where the named constant is \(c_2 = \dfrac{16\sqrt{3}}{15\pi} \approx 0.588\), so the multiplier of \(\sigma_b^2\) is \(c_2^{\,2} \approx 0.346\).

Thus \(|\boldsymbol\beta| < |\boldsymbol\beta^*|\) whenever \(\sigma_b^2 > 0\); attenuation grows with heterogeneity.

Important Note on the Attenuation Constant

The multiplier \(c_2^{\,2} \approx 0.346\) is an approximation specific to:

• Logistic GLMMs with random intercepts
• Normally distributed random effects

Source: Zeger, Liang & Albert (1988); \(c_2 = 16\sqrt{3}/(15\pi) \approx 0.588\).

Does NOT apply to: other link functions or random-slope models.

Clinical Translation

“Individual treatment effects look stronger than population effects because high responders and low responders average out.”

16.3 Attenuation Formula and Severity

\[\boldsymbol\beta_{\text{marginal}} \approx \frac{\boldsymbol\beta^*_{\text{GLMM}}}{\sqrt{1 + 0.346\,\sigma_b^2}}, \qquad c = \frac{1}{\sqrt{1 + 0.346\,\sigma_b^2}}\]

Worked example (\(\boldsymbol\beta^* = 1.5\), \(\sigma_b^2 = 4\)): \(\boldsymbol\beta_{\text{marg}} = 1.5/\sqrt{1 + 0.346(4)} = 1.5/1.544 \approx 0.97\) (factor \(c=0.65\)). Population OR \(\approx 2.6\) vs individual OR \(\approx 4.5\): same direction, weaker population effect.

Check Your Understanding: Attenuation

Question: A researcher reports that a GLMM for diabetes progression yields a treatment OR of 0.5 (protective effect) with random intercept variance \(\sigma_b^2 = 3.0\). They conclude: “Treatment cuts the population odds in half.”

What is wrong with this conclusion, and what is the correct population OR?

16.3 Log Link Exception

For log-link GLMMs: \[\log \{\operatorname{E}(Y_{ij}\mid b_i)\} = X_{ij}\boldsymbol\beta^* + Z_{ij} b_i\]

Marginal mean factorizes: \[\mu_{ij} = \exp(X_{ij}\boldsymbol\beta^*)\, \operatorname{E}\{\exp(Z_{ij} b_i)\}\]

For a single random intercept, \(\log(\mu_{ij}) = X_{ij}\boldsymbol\beta^* + g_{11}/2\): fixed effects not in \(Z_{ij}\) carry over unchanged (only the intercept shifts by \(g_{11}/2\)).

Important Clarification: Marginal Predictions from GLMMs

Common misconception: predict(glmer_fit, re.form = NA) gives marginal predictions.

Reality (FLW textbook, p. 477):

• re.form = NA gives conditional predictions at \(b_i = 0\) (the “typical” subject)
• For non-linear links: \(\operatorname{E}(Y \mid b = 0) \neq \operatorname{E}(Y)\)
• True marginal predictions require integrating over the random effects distribution

Helper Function for True Marginal Predictions

# CAVEAT: this helper assumes a SINGLE random intercept and a BINOMIAL family.
# It does NOT handle random slopes or non-binomial GLMMs.
glmm_marginal_prob <- function(fit, newdata, draws = 2000, seed = 667) {
  stopifnot(family(fit)$family == "binomial")
  set.seed(seed)  # reproducible Monte Carlo (course seed 667)
  # Fixed-part linear predictor at b = 0 (re.form = NA): X*beta-star only
  eta_fix <- predict(fit, newdata = newdata, re.form = NA, type = "link")
  # VarCorr(fit)[[1]] is the covariance of the FIRST grouping factor; its
  # "stddev" attribute is the random-intercept SD (a 1x1 here -> sigma_b).
  sd_b <- as.numeric(attr(VarCorr(fit)[[1]], "stddev"))
  # Monte Carlo integration over b ~ N(0, sigma_b^2): for each of `draws`
  # random intercepts, push eta_fix + b through plogis (= expit = h), then
  # average the probabilities across draws (rowMeans over the draw columns).
  rowMeans(sapply(1:draws, function(...) plogis(eta_fix + rnorm(1, 0, sd_b))))
}

Logistic Attenuation: Simulation Setup

set.seed(667)
# Subject-specific parameters (conditional on random effects)
beta0 <- -1.5        # Log-odds at time=0 for average individual
beta1 <- 0.75        # Subject-specific slope (change per unit time)
sigma_b <- 2         # SD of random intercepts (heterogeneity); variance = 4

# Time grid for plotting
time_grid <- seq(-4, 8, by = 0.25)

# Simulate 25 individuals with varying intercepts
subject_b <- rnorm(25, mean = 0, sd = sigma_b)

# Generate subject-specific curves
subject_curves <- tidyr::expand_grid(time = time_grid, b = subject_b) |>
  mutate(
    linear_predictor = beta0 + beta1 * time + b,
    prob = plogis(linear_predictor),
    curve_id = factor(round(b, 2))
  )

# Generate marginal curve by Monte Carlo integration
marginal_curve <- tibble(
  time = time_grid,
  prob = vapply(
    time_grid,
    function(t) mean(plogis(beta0 + beta1 * t + rnorm(2000, 0, sigma_b))),
    numeric(1)
  )
)

Logistic Attenuation: Visualization

logit_plot <- ggplot(subject_curves, aes(x = time, y = prob, group = curve_id)) +
  geom_line(alpha = 0.35, color = "#1f77b4") +
  geom_line(data = marginal_curve, aes(x = time, y = prob),
            inherit.aes = FALSE, color = "#d62728", linewidth = 1.2) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Time",
    y = "Probability of success",
    title = "Subject-specific vs marginal logistic curves",
    subtitle = "Random intercept variance = 4; slopes attenuate after averaging"
  ) +
  theme_minimal(base_size = 16)
print(logit_plot)

Effect of Heterogeneity on Marginal Curves

# Low heterogeneity (variance = 0.5)
sigma_b_low <- sqrt(0.5)
marginal_low <- tibble(
  time = seq(-4, 8, 0.5),
  prob = vapply(seq(-4, 8, 0.5),
                function(t) mean(plogis(beta0 + beta1 * t + rnorm(1000, 0, sigma_b_low))),
                numeric(1)),
  heterogeneity = "Low (var=0.5)"
)

# High heterogeneity (variance = 9)
sigma_b_high <- 3
marginal_high <- tibble(
  time = seq(-4, 8, 0.5),
  prob = vapply(seq(-4, 8, 0.5),
                function(t) mean(plogis(beta0 + beta1 * t + rnorm(1000, 0, sigma_b_high))),
                numeric(1)),
  heterogeneity = "High (var=9)"
)

# Plot comparison
comparison_plot <- ggplot() +
  geom_line(data = marginal_low, aes(x = time, y = prob),
            color = "#2ca02c", linewidth = 1.5) +
  geom_line(data = marginal_high, aes(x = time, y = prob),
            color = "#ff7f0e", linewidth = 1.5) +
  labs(title = "Impact of heterogeneity on marginal curves",
       subtitle = "Green: Low heterogeneity (steep) | Orange: High heterogeneity (flat)",
       x = "Time", y = "Probability") +
  theme_minimal(base_size = 14)
print(comparison_plot)

Logistic Attenuation: Key Observations

16.4 Simple Numerical Illustration

Consider a simple setting: binary responses depend only on time with a random intercept.

GLMM (conditional model): \[\text{logit}\{P(Y_{ij} = 1 \mid b_i)\} = \beta_0^* + \beta_1^* t_{ij} + b_i, \quad b_i \sim N(0, \sigma_b^2)\]

Key insight: The conditional slope \(\beta_1^*\) describes how an individual’s log-odds changes over time.

16.4 Why Marginal Slope Differs

To obtain the marginal (population-average) probability, we must integrate: \[P(Y_{ij} = 1) = \int \text{expit}(\beta_0^* + \beta_1^* t_{ij} + b_i) \, \phi(b_i; 0, \sigma_b^2) \, db_i\]

Critical point: Due to the non-linearity of expit(), this integral does NOT simplify to: \[\text{logit}\{P(Y_{ij} = 1)\} = \beta_0 + \beta_1 t_{ij}\]

with \(\beta_1 = \beta_1^*\). Instead, \(|\beta_1| < |\beta_1^*|\).

16.4 The b=0 Misconception

Common error: Claiming that setting \(b_i = 0\) in a GLMM gives marginal predictions.

Why this is wrong (FLW textbook, p. 477): \[\operatorname{E}[Y \mid b = 0] \neq \operatorname{E}[Y]\]

For non-linear links, the conditional expectation at the mean of \(b\) is NOT the marginal expectation.

16.4 Demonstration: b=0 vs True Marginal

# Parameters
beta0 <- 0
beta1 <- 1
sigma_b <- 2  # variance = 4

# Time points
t_vals <- seq(-3, 3, by = 0.5)

# Prediction at b = 0 (the WRONG "marginal")
prob_b0 <- plogis(beta0 + beta1 * t_vals)

# True marginal (integrate over random effects)
prob_marginal <- sapply(t_vals, function(t) {
  mean(plogis(beta0 + beta1 * t + rnorm(5000, 0, sigma_b)))
})

# Compare
comparison_df <- data.frame(
  time = t_vals,
  prob_b0 = prob_b0,
  prob_marginal = prob_marginal,
  difference = prob_b0 - prob_marginal
)

cat("Comparison of P(Y=1|b=0) vs true E[Y]:\n")

16.4 Visual Comparison

ggplot(comparison_df, aes(x = time)) +
  geom_line(aes(y = prob_b0), color = "#1f77b4", linetype = "dashed", linewidth = 1.2) +
  geom_line(aes(y = prob_marginal), color = "#d62728", linewidth = 1.2) +
  annotate("text", x = 2, y = 0.9, label = "P(Y=1 | b=0)", color = "#1f77b4", size = 5) +
  annotate("text", x = 2, y = 0.65, label = "True E[Y]", color = "#d62728", size = 5) +
  labs(
    x = "Time",
    y = "Probability",
    title = "Conditional at b=0 vs True Marginal Probability",
    subtitle = expression(paste("GLMM with ", sigma[b], " = 2; b=0 overestimates probabilities near 0.5"))
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  theme_minimal(base_size = 14)

Why Is the Difference Largest Near p = 0.5?

Difference Near p = 0.5: When It Matters

Takeaway: the marginal-vs-conditional gap is largest for steep regions of the link and large \(\sigma_b^2\).

16.4 Key Takeaway

The difference matters most when:

• \(\sigma_b^2\) is large (high heterogeneity)
• Predictions are near probability 0.5 (steepest part of logistic curve)

This bridges to the ECG crossover case study where we compare GEE vs GLMM estimates.

16.5 ECG Crossover Overview

Study design (FLW Table 16.1; Jones & Kenward 1989):

• 67 patients randomized to sequences: Placebo\(\to\)Active (34) or Active\(\to\)Placebo (33)
• Binary response each period: abnormal ECG (1) vs normal (0)
• Sampling zero in Sequence 1, pattern (1,0) - sparse cell issue

Goal: Compare treatment effect estimates from marginal vs mixed models.

ECG Crossover: Key Questions

16.5 Crossover Design: Sequence x Period -> Treatment

                Period 1        Period 2
            +-------------+   +-------------+
 Seq P->A   |  Placebo    |   |  Active     |   treatment = 0, then 1
 (34 pts)   |  (trt = 0)  |   |  (trt = 1)  |
            +-------------+   +-------------+
 Seq A->P   |  Active     |   |  Placebo    |   treatment = 1, then 0
 (33 pts)   |  (trt = 1)  |   |  (trt = 0)  |
            +-------------+   +-------------+

Each subject contributes two binary ECG responses; treatment is the active/placebo indicator for that period, decoded from (sequence, period).

16.5 Data Reconstruction

ecg_counts <- tribble(
  ~sequence, ~resp_p1, ~resp_p2, ~count,
  "P-A", 1, 1, 6,
  "P-A", 1, 0, 0,  # Sampling zero
  "P-A", 0, 1, 6,
  "P-A", 0, 0, 22,
  "A-P", 1, 1, 9,
  "A-P", 1, 0, 4,
  "A-P", 0, 1, 2,
  "A-P", 0, 0, 18
)

ecg_long <- ecg_counts |>
  tidyr::uncount(count) |>
  # one UNIQUE subject id per expanded row (67 patients, 2 rows each).
  # Do NOT reuse uncount's within-cell .id: that repeats ids across cells.
  mutate(id = factor(row_number())) |>
  tidyr::pivot_longer(
    cols = c(resp_p1, resp_p2),
    names_to = "period",
    values_to = "response"
  ) |>
  mutate(
    period = if_else(period == "resp_p1", 1L, 2L),
    # decode active-drug indicator from (sequence, period); see schematic
    treatment = case_when(
      sequence == "P-A" & period == 1L ~ 0L,
      sequence == "P-A" & period == 2L ~ 1L,
      sequence == "A-P" & period == 1L ~ 1L,
      TRUE ~ 0L
    )
  ) |>
  arrange(id, period)

cat("Subjects:", length(unique(ecg_long$id)),
    "| rows:", nrow(ecg_long), "\n\n")

16.5 Model Fits: GEE

# Marginal model (GEE), exchangeable working correlation
gee_fit <- geepack::geeglm(
  response ~ treatment + period,
  id = id,
  family = binomial,
  corstr = "exchangeable",
  data = ecg_long
)
# Note: only 2 obs/subject, so robust SEs can be biased downward
# (small-sample corrections shown later).
summary(gee_fit)

16.5 Model Fits: GLMM

# Mixed model (GLMM): 100-point adaptive Gaussian quadrature (FLW Table 16.3).
# Laplace (nAGQ = 1) is inadequate here because Var(b) is very large (~24).
glmm_fit <- lme4::glmer(
  response ~ treatment + period + (1 | id),
  family = binomial,
  data = ecg_long,
  nAGQ = 100,
  control = lme4::glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)
summary(glmm_fit)

16.5 The Same Two Fits in SAS

proc genmod data=ecg descending;
  class id period;
  model response = treatment period / dist=bin link=logit;
  repeated subject=id / type=exch;   /* exchangeable working correlation */
run;

proc glimmix data=ecg method=quad(qpoints=100);  /* adaptive Gaussian quadrature */
  class id;
  model response = treatment period / dist=bin link=logit solution;
  random intercept / subject=id;
run;

This slide is reference only (not executed). qpoints=100 mirrors FLW Table 16.3; type=exch mirrors the working-correlation GEE.

16.5 Extract and Compare Treatment Effects

# Extract treatment effects
gee_coef <- broom.mixed::tidy(gee_fit) |>
  filter(term == "treatment") |>
  mutate(model = "Marginal (GEE)", estimand = "Population odds ratio")

glmm_coef <- broom.mixed::tidy(glmm_fit, effects = "fixed") |>
  filter(term == "treatment") |>
  mutate(model = "Mixed (GLMM)", estimand = "Subject-specific odds ratio")

# Combine results
contrast_tbl <- bind_rows(gee_coef, glmm_coef) |>
  transmute(
    model, estimand,
    log_odds = estimate,
    se = std.error,
    odds_ratio = exp(estimate),
    lower = exp(estimate - 1.96 * se),
    upper = exp(estimate + 1.96 * se)
  )
print(contrast_tbl)

16.5 Demonstrating True Marginal Predictions

# Create prediction data
nd <- data.frame(treatment = c(0L, 1L), period = 1L)

# Conditional prediction at b_i = 0 (typical subject)
nd$prob_conditional <- predict(glmm_fit, newdata = nd, re.form = NA, type = "response")

# True marginal prediction (averaging over random effects)
nd$prob_marginal <- glmm_marginal_prob(glmm_fit, nd)

# GEE marginal prediction for comparison
nd$prob_gee <- predict(gee_fit, newdata = nd, type = "response")

cat("Predictions for Period 1:\n")

16.5 Treatment Effect Comparison

Results interpretation (from computed values):

Why the difference? Marginal model averages over heterogeneity; mixed model conditions on \(b_i\). With \(\widehat{\operatorname{Var}}(b_i)\) very large and only 2 obs/subject, the GLMM OR is a point estimate from a weakly identified \(\sigma_b\), so read it as illustrative of direction and order of magnitude.

Interpreting the ECG Results: A Template

For the FDA (population-average perspective):

“Active drug increased the population odds of abnormal ECG by approximately 77% (OR = 1.77, 95% CI: 1.12 to 2.79). This population-level estimate accounts for the average effect across all patients regardless of their baseline cardiac status.”

For the clinician (subject-specific perspective):

“For an individual patient, the active drug raises that patient’s personal odds of an abnormal ECG by a substantial factor (OR \(\approx\) 6.44, 95% CI: 1.05 to 39.64), holding their baseline cardiac propensity constant. The wide interval reflects very large, imprecisely estimated between-patient variability.”

Key point: Neither interpretation is “wrong” – they answer different questions! Report the conditional OR with its uncertainty, not as a precise percentage.

Non-collapsibility in Numbers (no confounding)

Two subgroups of equal size, conditional (within-subgroup) OR = \(9\) in both (so nothing is confounded):

Marginal OR \(= 5.44 \neq 9 =\) conditional OR, with no confounder. The odds ratio is non-collapsible\(^\dagger\).

Which to Report?

Both are correct – they answer different questions!

16.5 Small-Sample GEE Variance Adjustments

With only 2 observations per subject, standard GEE robust SEs can be biased downward.

# Option 1: Jackknife standard errors
gee_jack <- geepack::geeglm(
  response ~ treatment + period, id = id,
  family = binomial, corstr = "exchangeable",
  data = ecg_long, std.err = "jack"
)

# Option 2: CR2 correction via clubSandwich
library(clubSandwich)
gee_cr2 <- coef_test(gee_fit, vcov = "CR2", cluster = ecg_long$id)

# Option 3: Fay correction
gee_fay <- geepack::geeglm(
  response ~ treatment + period, id = id,
  family = binomial, corstr = "exchangeable",
  data = ecg_long, std.err = "fay"
)

Note: These corrections typically increase SEs by 10-30% with small clusters.

16.5 Structured Comparison

16.5 OR Contrast Visualization

or_plot <- ggplot(contrast_tbl, aes(x = model, y = odds_ratio)) +
  geom_pointrange(aes(ymin = lower, ymax = upper), size = 0.9, color = "#1f65b7") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray40") +
  scale_y_log10(breaks = c(0.5, 1, 2, 4, 8, 16)) +
  coord_flip() +
  labs(
    x = NULL,
    y = "Treatment odds ratio (Active vs Placebo)",
    title = "Population vs subject-specific treatment effects",
    subtitle = "Same data, different estimands -- both valid for their intended use"
  ) +
  theme_minimal(base_size = 16)
print(or_plot)

16.5 The Sampling Zero Problem

Observation: ECG data has ZERO in cell (Sequence P-A, response pattern 1,0).

Why this matters:

FLW’s sensitivity refits (add \(1/4\) to the (P-A, 1,0) cell): marginal Treatment \(\approx 0.55\), GLMM Treatment \(\approx 1.64\), \(\widehat{g}_{11} \approx 20.7\) - substantive conclusions unchanged.

Alternatives to Continuity Correction

Key point: The 0.25 correction is a diagnostic probe, not a solution.

16.5 Robustness Check with Conditional Logistic

# Conditional logistic (fixed-effects / within-subject), stratified by subject.
# In a 2-period crossover, treatment and period are confounded WITHIN subject,
# so we condition on subject and estimate the treatment effect alone; adding
# period to the conditional model is not separately identified here.
clogit_fit <- survival::clogit(
  response ~ treatment + strata(id),
  data = ecg_long
)

clogit_or <- exp(coef(clogit_fit)["treatment"])
cat("Conditional logistic OR for treatment:", round(clogit_or, 2), "\n")

16.5 Additional Diagnostics

For sparse data with small \(n_i\):

Key Messages for Sparse Data

• Distributional assumptions on \(b_i\) matter most when \(n_i\) is small
• With only 2 observations/subject, both approaches have limitations
• Report sensitivity analyses alongside main results

16.6 Conclusion Highlights

Key insight: Mixed models do NOT automatically provide marginal summaries without extra integration.

16.6 No Grand Unified Model

FLW’s verdict: do not chase one model that answers every question. Different scientific questions demand different models.

“The megalomaniacal strategy of fitting a grand unified model… is dangerous, unnecessary and counterproductive.” – Drum & McCullagh (1993), quoted approvingly by FLW (p. 486)

One size does not fit all. Choose marginal or mixed by the question, and report both when the audience needs both.

Why GLMMs Do Not Solve Everything

Bottom line: One size does not fit all.

When to Choose Marginal vs Mixed Models

Can you fit both? YES! Report both when stakeholders need both estimands.

Common Misconceptions

More Misconceptions

Common Software Pitfalls

Common Mistakes When Contrasting Models

Advanced (beyond FLW): Marginalized Mixed Models

The challenge: want BOTH random effects (heterogeneity) AND a direct marginal interpretation. FLW notes this is possible only with a non-standard conditional-mean specification (footnote, p. 478); the constructions themselves are beyond Ch. 16.

Trade-off: best of both worlds in theory, but complex, specialized software, computationally intensive.

Recipe Card: Decision Table

Step 1 - Define the estimand, then read across:

Recipe Card: Decision Flowchart

                    START
                       |
            What is your estimand?
           /                      \
    Population effect          Individual effect
          |                           |
      Use GEE                     Use GLMM
          |                           |
   Report population OR/RR     Report conditional OR/RR
          |                           |
   Consider: robust SE,        Consider: nAGQ, random-
   working correlation, MCAR   effects structure, MAR
          |                           |
          +----------+----------------+
                     |
        Report BOTH if the audience includes
        clinicians AND regulators

Methods-section template: “We used [GEE/GLMM] because our primary interest was [population-average/subject-specific] effects.”

16.7 Further Reading

Most accessible:

Comprehensive coverage:

Further Reading: Technical Foundations

Summary Checklist

Before analyzing longitudinal/clustered data, ask:

• What is my primary research question?
• Who is my audience (regulator, clinician, patient)?
• Do I need population or individual inference?
• Is my outcome linear or non-linear?
• How much heterogeneity exists (\(\sigma_b^2\))?
• Do I have sparse cells or small clusters?
• Am I computing marginal predictions correctly?
• Should I report both marginal and conditional estimates?

Final Thought

Understanding the distinction between marginal and conditional models is essential for correct inference in longitudinal data analysis.