BIOS 667 - Lecture 12: Marginal Models (Ch. 12)

Lecture Objectives

• Understand why standard GLMs fail for longitudinal data
• Define marginal models and when to use them
• Distinguish marginal vs mixed model interpretations of \(\beta\)
• Describe the three-part marginal model: mean, variance, association
• Preview GEE estimation (Ch. 13) and apply it to a real Case Study

Prerequisites: What You Should Recall

Today we glue these together: a GLM mean model plus a within-subject association structure, fit so that inference stays valid.

Notation reference

Notation (this lecture): GEE working covariance

This is the deck that introduces the GEE working-covariance notation. Reserved symbols (\(R_i\) = within-subject residual covariance) are unchanged.

Motivation

Repeated measures within a subject are correlated.

A naive GLM assumes independence, so its inference is invalid: the standard errors are incorrect (FLW pp. 43-44).

For discrete data, how we model correlation also changes the meaning of \(\beta\).

Two Model Families (FLW Ch. 12)

FLW Ch. 12 contrasts two families: marginal (population-average) vs mixed (subject-specific).

Marginal Model = Mean + Variance + Association

Mean:

\[E(Y_{ij}\mid X_{ij})=\mu_{ij}, \quad g(\mu_{ij})=X_{ij}\beta\]

Variance:

\[\operatorname{Var}(Y_{ij}\mid X_{ij})=\phi\,v(\mu_{ij})\]

Association:

Within-subject dependence via the working covariance \(V_i(\alpha)\) (exchangeable, AR(1), etc.).

Key Insight: Working Correlation

Changing \(V_i(\alpha)\) affects precision, not the interpretation of \(\beta\).

This is the defining feature of marginal models: \(\beta\) targets the same population-average quantity regardless of the working correlation, and robust (sandwich) SEs stay valid even if \(V_i(\alpha)\) is misspecified.

12.2 Notation and Data Structure

• \(N\) subjects; subject \(i\) has \(n_i\) measurements at times \(t_{ij}\)
• Stack responses: \(\mathbf{Y}_i=(Y_{i1},\dots,Y_{in_i})^\top\)
• Across-subject independence
• Within-subject dependence

Covariates

\(X_{ij}\) is the \(1\times p\) covariate row; stack the rows to \(X_i\) (\(n_i\times p\))

Discrete Outcomes: Correlation Is Constrained

For binary outcomes the feasible correlation depends on the marginal means. With \(p_1=E(Y_1)=0.2\) and \(p_2=E(Y_2)=0.8\):

The largest possible \(\Pr(Y_1=1, Y_2=1)\) is \(\min(p_1,p_2)=0.2\) (the Frechet bound):

At that maximum, the correlation is \[\rho_{12}^{\max}=\sqrt{\frac{p_1(1-p_2)}{p_2(1-p_1)}}=\sqrt{\frac{0.2\cdot 0.2}{0.8\cdot 0.8}}=\sqrt{\tfrac{0.04}{0.64}}=0.25.\]

So \(\rho_{12}\le 0.25\): the means alone cap the correlation.

Discrete Outcomes: Odds Ratios Are More Natural

For a pair of binary outcomes, define the four joint cell probabilities:

The pairwise log odds ratio uses these cells directly:

\[ \theta_{jk}=\log\frac{P_{11}\,P_{00}}{P_{10}\,P_{01}} \]

Unlike the correlation, the odds ratio is not mechanically bounded by the marginal means, so it is the more natural association measure for binary data.

Master Simulated Dataset: continuous outcome

A small simulated dataset to illustrate the long format and the continuous fit. The continuous outcome carries subject random effects (\(b_0, b_1\)), so it is genuinely correlated within subject.

suppressPackageStartupMessages({
  library(dplyr); library(tidyr); library(tibble)
  library(ggplot2); library(nlme); library(broom); library(geepack)
})

set.seed(667)
N <- 120
n_i  <- sample(3:6, N, replace = TRUE)
id   <- rep(1:N, times = n_i)
time <- unlist(lapply(n_i, function(n) 0:(n - 1)))
trt  <- rep(rbinom(N, 1, 0.5), times = n_i)

# subject random effects induce within-subject correlation
b0 <- rep(rnorm(N, 0, 1.2), times = n_i)
b1 <- rep(rnorm(N, 0, 0.4), times = n_i)
y  <- 10 + 0.6 * time + 1.0 * trt + 0.5 * (trt * time) +
      b0 + b1 * time + rnorm(length(id), 0, 1)

Master Simulated Dataset: discrete outcomes

# Discrete outcomes generated INDEPENDENT within subject (illustration only)
lp_bin  <- -1 + 0.25 * time + 0.7 * trt
y_bin   <- rbinom(length(id), 1, plogis(lp_bin))

mu_cnt  <- exp(1.2 + 0.15 * time + 0.3 * trt)
y_count <- rpois(length(id), mu_cnt)

cuts <- quantile(y, probs = c(.25, .5, .75))
ordinal_y <- cut(y, breaks = c(-Inf, cuts, Inf),
                 labels = c("L", "M", "H", "VH"), ordered_result = TRUE)

dat <- tibble(
  id   = factor(id),
  time = as.integer(time),
  trt  = factor(trt, levels = c(0, 1), labels = c("Ctl", "Tx")),
  y = y, y_bin = as.integer(y_bin),
  y_count = as.integer(y_count), ordinal_y = ordered(ordinal_y)
)

Long Format: One Subject’s Rows

head(dat, 8)

Visualizing Imbalance

ggplot(dat, aes(time, y, group = id, color = trt)) +
  geom_line(alpha = 0.15) +
  stat_summary(aes(group = trt), fun = mean, geom = "line", linewidth = 1.2) +
  labs(title = "Unbalanced Trajectories by Treatment (continuous y)",
       y = "Continuous outcome y") +
  theme_minimal()

Implicit Mean Assumption

The marginal mean depends only on the current covariate:

\[ E(Y_{ij}\mid X_{i1},\dots,X_{in_i})=E(Y_{ij}\mid X_{ij}) \]

This holds for design-fixed covariates (fixed by the study design, e.g. randomized treatment, planned visit time).

   time ->        X_1  ->  Y_1     X_2  ->  Y_2     X_3  ->  Y_3
                            |                |
   FORBIDDEN feedback:      +----> X_2        +----> X_3   (Y_(j-1) -> X_j)

It fails when a future covariate responds to a past outcome (feedback), e.g. a clinician raising the dose after a poor prior response.

Outcome-Type Specifications (one table)

All four use the same three-part recipe; only the link and variance function change.

Association in every row: a working covariance \(V_i(\alpha)\) (exchangeable, AR(1), unstructured, …).

Ordinal Specification (Proportional Odds)

Genuinely new structure: cumulative dichotomizations with shared slopes.

\[\Pr(Y_{ij}\le k\mid X_{ij})=F_{ijk}, \qquad \log\frac{F_{ijk}}{1-F_{ijk}}=\kappa_k - X_{ij}\beta\]

• \(\kappa_k\): category-specific thresholds (cut-points), \(\kappa_1<\dots<\kappa_{K-1}\) (FLW writes these \(\alpha_k\); we use \(\kappa_k\) to keep \(\alpha\) for the GEE working correlation)
• \(\beta\): population-average log cumulative-odds ratios, shared across cut-points (the proportional-odds assumption); the minus-sign form matches L11 and R’s polr, so a positive \(\beta\) pushes toward higher categories

GEE bridge: how to read geeglm (estimation is L13)

Continuous Outcome: Fit (gls)

dat_cont <- dat |> dplyr::select(id, time, trt, y)

fit_cont <- gls(
  y ~ time * trt,
  data = dat_cont,
  correlation = corAR1(form = ~ time | id),
  method = "REML"
)

coef(summary(fit_cont))

Continuous: Plain-Language Interpretation

In a sentence: both groups increase over time, but the treatment group starts higher and increases faster, so the separation grows.

Case Study: Epilepsy Seizure Counts

Case Study: Load and Reshape (wide -> long)

epi_url  <- "https://content.sph.harvard.edu/fitzmaur/ala2e/epilepsy.txt"
epi_file <- "../../data/epilepsy.txt"
if (!file.exists(epi_file)) {
  dir.create(dirname(epi_file), showWarnings = FALSE, recursive = TRUE)
  download.file(epi_url, epi_file)
}

# FLW .txt has a comment header; data rows start at the first "<int> <int> ..."
raw <- readLines(epi_file)
data_start <- grep("^\\s*[0-9]+\\s+[0-9]", raw)[1]
epi_wide <- read.table(text = raw[data_start:length(raw)])
names(epi_wide) <- c("id", "trt", "age", "baseline", "c1", "c2", "c3", "c4")

epi <- epi_wide |>
  tidyr::pivot_longer(c(c1, c2, c3, c4),
                      names_to = "visit", values_to = "seizures") |>
  dplyr::mutate(
    visit = as.integer(gsub("c", "", visit)),
    trt   = factor(trt, levels = c(0, 1),
                   labels = c("Placebo", "Progabide")),
    lbase = log(baseline / 4)   # 8-wk baseline -> per-2-wk rate scale
  )

head(epi, 6)

Case Study: Marginal Poisson GEE

fit_epi <- geeglm(
  seizures ~ trt + lbase,
  id      = id,
  data    = epi,
  family  = poisson(link = "log"),
  corstr  = "exchangeable"
)

summary(fit_epi)

Case Study: Rate Ratios and Interpretation

est <- coef(summary(fit_epi))
rr  <- exp(cbind(RR = est[, "Estimate"],
                 Lower = est[, "Estimate"] - 1.96 * est[, "Std.err"],
                 Upper = est[, "Estimate"] + 1.96 * est[, "Std.err"]))
round(rr, 3)

Interpretation (rate ratios, population-average):

• Progabide vs placebo: \(\widehat{RR}\approx 0.90\) (95% CI crosses 1) - no significant reduction in seizure rate once baseline is accounted for.
• Log baseline: \(\widehat{RR}\approx 3.2\) per unit - the dominant predictor; higher baseline rate strongly predicts more post-randomization seizures.

Case Study: Model-Based vs Robust SEs

# Naive Poisson GLM ignores clustering: model-based SEs
fit_naive <- glm(seizures ~ trt + lbase, data = epi,
                 family = poisson(link = "log"))

data.frame(
  term       = c("trtProgabide", "lbase"),
  naive_SE   = summary(fit_naive)$coef[c("trtProgabide", "lbase"), "Std. Error"],
  robust_SE  = coef(summary(fit_epi))[c("trtProgabide", "lbase"), "Std.err"]
) |> dplyr::mutate(ratio = round(robust_SE / naive_SE, 2))

The Same Model in SAS

Static reference (not run here): the marginal Poisson GEE in PROC GENMOD with a REPEATED statement.

proc genmod data=epi;
   class id trt;
   model seizures = trt lbase / dist=poisson link=log;
   repeated subject=id / type=exch corrw;
run;

Why Full Likelihood Is Hard

needed_terms <- function(m) 2^m - m - 1
dfc <- data.frame(m = 4:12, extra_assoc = sapply(4:12, needed_terms))

ggplot(dfc, aes(m, extra_assoc)) +
  geom_point() + geom_line() +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Explosion of Association Parameters",
       x = "Number of repeats (m)",
       y = "Parameters beyond means (2^m - m - 1)") +
  theme_minimal()

Full Likelihood vs Marginal + GEE

GEE wins for practicality and interpretability.

Practical Decisions: Model Selection

Looking Ahead (beyond Ch. 12): Missing Data

Common Mistakes in Marginal Models

Key Takeaways

• Marginal model = mean + variance + association (\(g\), \(v(\mu)\), \(V_i(\alpha)\))
• \(\beta\) = population-average effect (log rate/odds ratio for count/binary), robust to the working correlation
• The working correlation \(V_i(\alpha)\) affects precision, not interpretation; robust SEs stay valid
• Case Study (epilepsy): progabide \(\widehat{RR}\approx 0.90\) (not significant); robust SEs ~4x the naive ones
• For discrete data, full joint likelihoods are intractable; use GEE (L13)
• Standard GEE assumes MCAR; missing-data extensions in L17

Further Reading

• FLW (2011) Ch. 12, “Marginal Models: Introduction and Overview” (pp. 341-352), and its Further Reading section.
• Liang, K.-Y. & Zeger, S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73, 13-22. (The original GEE paper; estimation in L13.)
• Robins, J.M., Rotnitzky, A. & Zhao, L.P. (1995). Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. JASA 90, 106-121. (IPW-GEE; missing data in L17.)
• Thall, P.F. & Vail, S.C. (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46, 657-671. (The epilepsy data.)

Next lecture (L13): GEE estimation - the estimating equations, the sandwich variance, and choosing the working correlation.