Marginal Models and Generalized Estimating Equations

12.1 Two Questions Hiding in One Regression

Consider a smoking-cessation program and two apparently identical questions. The first asks how much the risk of relapse would fall if the whole population quit for a month; the second asks how much a particular smoker’s risk would fall if that smoker quit. The first is a population-average or marginal question, comparing average outcomes across states of the population; the second is a subject-specific or conditional question, comparing one individual’s outcome to their own counterfactual. For a continuous outcome analyzed with an identity link the two questions have the same answer, because averaging is linear and the average of individual changes equals the change in the average. For a binary outcome analyzed with a logistic link they do not, and the reason is geometric. Figure 12.1 shows a collection of individual logistic curves, each steep, that differ only in a random intercept, together with their average across the population. The average curve is visibly flatter than any individual curve, because averaging S-shaped curves that are horizontally shifted produces a shallower S, and a flatter curve means a smaller slope. The population-average logistic coefficient is therefore attenuated relative to the subject-specific one.

Note. Each thin curve is one person’s logistic probability of an outcome as a function of a predictor, the curves differing only in a random intercept. Their average across the population, the marginal curve in red, is shallower than any individual curve. The population-average slope is thus attenuated relative to the subject-specific slope; the two estimands differ for nonlinear links.

The attenuation is quantifiable. For a logistic model with a random intercept of variance \(\sigma_u^2\), the population-average slope is approximately the subject-specific slope divided by \(\sqrt{1 + c^2\sigma_u^2}\), where \(c = 16\sqrt{3}/(15\pi) \approx 0.588\) arises from approximating the logistic by a probit. Figure 12.2 confirms the relation by simulation: the ratio of the population-average to the subject-specific slope, estimated by fitting both models to data generated with a known random-intercept variance, tracks the formula closely and falls steadily as the between-person variance grows, from equality when everyone shares the same intercept to roughly one-half when the between-person variance is large. The practical consequence is that a logistic regression coefficient means different things depending on which model produced it, and the two must not be compared or pooled as though they were the same quantity. Table 12.1 states the distinction, and the mapping to research questions is the section’s payload: program evaluation, epidemiological, and policy questions, which concern what happens to a population, call for the population-average estimand, whereas questions about individual process and mechanism, which concern what happens within a person, call for the subject-specific one.

Note. The ratio of the population-average to the subject-specific logistic slope, from simulation (points), against the random-intercept variance, with the approximation \(1/\sqrt{1 + c^2\sigma_u^2}\) (curve). The larger the between-person variance, the more the marginal slope is attenuated below the subject-specific slope. The two coincide only when there is no between-person variation.

12.2 The Machinery of Generalized Estimating Equations

A generalized estimating equation is specified in three independent pieces (Liang & Zeger, 1986). A mean model links a function of the average outcome to the predictors, exactly as in an ordinary generalized linear model, so that for the binary remission outcome of the therapy trial a logistic mean model expresses the log-odds of remission as a function of arm and week. A variance function states how the variance of the outcome depends on its mean, again as in a generalized linear model. And a working correlation \(R(\alpha)\) describes the within-person dependence among the repeated measurements, the piece an ordinary generalized linear model lacks. Figure 12.3 shows the standard choices: independence, which ignores the within-person correlation; exchangeable, which assumes a single common correlation between any two occasions; autoregressive, which lets the correlation decay with the separation in time; and unstructured, which estimates every pairwise correlation freely. The remarkable and defining property of the method, derived in the Foundations box, is that the estimate of the mean-model coefficients is consistent even when the working correlation is wrong, because the estimating equation for the mean is unbiased regardless of the assumed dependence. The working correlation affects only efficiency, not validity, which is why it is called “working”: it is a device to improve precision, not an assumption the conclusions rest on. This is confirmed empirically in the therapy trial, where the estimated arm-by-week effect on remission is essentially the same, near \(0.14\) on the log-odds scale, whether independence, exchangeable, or autoregressive correlation is assumed, while the standard errors differ slightly. Table 12.2 guides the choice by design.

Note. Four working-correlation structures as heatmaps: independence (no within-person correlation), exchangeable (one common correlation), autoregressive (decaying with time separation), and unstructured (every pair free). The coefficient estimates are consistent under any of these; the choice affects only efficiency, and the structure should match the design.

The robustness of the sandwich estimator is not free, and its cost is a small-sample fragility: the meat of the sandwich is an average of residual products across clusters, and when the number of clusters is small that average is a poor and downward-biased estimate of the true variance, so the standard errors are too small and confidence intervals too narrow. Figure 12.4 shows the consequence in a coverage simulation: with a hundred clusters the ninety-five-percent interval covers at nearly its nominal rate, but with ten or fifteen clusters it covers only ninety-two percent of the time, an inflation of the true error rate that would make many published small-sample generalized-estimating-equation analyses anticonservative. The remedies are the bias-corrected sandwich estimators, which adjust the residuals to remove the downward bias, and the practical guidance is to apply such a correction whenever the number of clusters is below roughly forty to fifty (Mancl & DeRouen, 2001; Kauermann & Carroll, 2001). Table 12.3 lists the corrections and when each is required.

Note. Empirical coverage of the nominal ninety-five-percent confidence interval for a time slope, from a simulation, against the number of clusters. With few clusters the robust standard error is downward biased and the interval covers below its nominal rate; the shortfall is negligible by about fifty clusters. Small-sample corrections are needed below roughly forty.

12.3 Model Evaluation and Selection

Because a generalized estimating equation is not fitted by maximizing a likelihood, the familiar likelihood-based tools are unavailable, and their closest replacement is the quasi-likelihood information criterion (QIC) of Pan (2001), which adapts the Akaike criterion to the quasi-likelihood setting and can be used to compare working-correlation structures and, in a variant, mean models. The criterion must be used judiciously and never worshipped, because it is not a true likelihood and its behavior is less well understood than the Akaike criterion it imitates; it is a rough guide to the working correlation, not an arbiter of substantive models. The absence of a likelihood also means there is no natural coefficient of determination, and the honest way to communicate a generalized-estimating-equation fit is through predicted quantities on the scale of interest, the population-average probabilities or means implied by the model, rather than through a summary statistic that the method does not naturally provide.

12.4 Missing Data: The Method’s Soft Spot

The generalized estimating equation carries a liability that is easy to overlook and directly opposed to a common misconception. Because it makes no full distributional assumption, it is often supposed to handle missing data gracefully, but the truth is the reverse: the unweighted estimating equation is valid only when the data are missing completely at random, a far stronger requirement than the missing-at-random condition under which the likelihood methods of Chapter 6 remain valid. Under the missing-at-random dropout typical of longitudinal studies, where those doing worse are more likely to leave, the unweighted method is biased, because it effectively averages over the selected sample of those who remain. Figure 12.5 is the chapter’s signature exhibit, a known-truth simulation in which a symptom declines over time and dropout depends on the previous observed value, so that sicker patients leave and the survivors are systematically healthier. The unweighted generalized estimating equation, tracking the observed survivors, reports a mean trajectory that falls well below the truth and so overstates improvement, whereas the mixed model estimated by full-information maximum likelihood recovers the true trajectory almost exactly, because likelihood is valid under missing-at-random. The remedy within the marginal framework is the inverse-probability-weighted generalized estimating equation (Robins et al., 1995), which weights each observed record by the inverse of its estimated probability of still being in the study, upweighting the kinds of people who tend to drop out so that the observed sample is made to represent the full population. Figure 12.5 shows that this weighting repairs part of the bias, but only part, and the reason is instructive: the weighting uses only the marginal observed data, discarding the within-person information about where a departed person’s trajectory was heading that the likelihood method exploits in full. The exhibit therefore delivers two lessons at once, that the marginal method is fragile under realistic missingness, and that even its remedy is dominated by the likelihood approach that the rest of the book favors.

Note. A known-truth simulation in which a symptom declines and sicker patients drop out (missing at random, depending on the previous observed value). The unweighted generalized estimating equation (orange) falls below the true cohort mean (dashed) and overstates improvement; the mixed model by maximum likelihood (blue) recovers the truth; the inverse-probability-weighted estimating equation (green) repairs part of the bias but not all, because it discards the within-person information likelihood uses.

12.5 Generalized Estimating Equations in Practice

The primary worked example fits the population-average probability of remission in the therapy trial as a function of arm and week. The model reports how the average remission rate in each arm evolves, and Figure 12.6 displays those model-implied population-average probabilities: both arms improve, and the treatment arm’s remission probability rises faster, the arm-by-week interaction being the population-average treatment effect. This is exactly the quantity a program evaluation wants, the change in the population rate, and it is reported as a probability rather than as a log-odds coefficient precisely because the population-average probability is interpretable while the attenuated coefficient invites the subject-specific misreading warned against above. A second example, a day-level binary outcome in the experience-sampling data modeled on within-person predictors, would illustrate the honest limit of the approach: when the interest is genuinely in the within-person process, whether a person is more likely to act a certain way on days when they are more stressed than usual, the population-average estimand answers a subtly different question, and the subject-specific mixed model of Chapters 15 and 23 is the better tool. Table 12.4 is the reporting checklist.

Note. Model-implied population-average probability of remission in the therapy trial for each arm across weeks, from a logistic generalized estimating equation. Both arms improve and the treatment arm improves faster; the arm-by-week difference is the population-average treatment effect, the quantity a program evaluation asks about.

12.6 Choosing Between Marginal and Mixed Models

The choice between a generalized estimating equation and a mixed model is not a matter of statistical taste but of estimand and assumptions, and Figure 12.7 lays out the decision. The first and governing question is which estimand the research wants: a population-average answer points to the marginal model, a subject-specific answer to the mixed model, and for a linear outcome the two coincide so either will do. The second question is the missingness mechanism: data missing at random, the common case, are handled validly by the likelihood-based mixed model but require weighting for the marginal model. The third concerns what the analysis must deliver: only the mixed model provides estimates of the random-effect variances, the individual differences in change that are often the substantive quantity of interest, and their absence from the marginal model is both a simplification and a genuine loss. The remaining considerations, the number of clusters and computational robustness, favor the marginal model when clusters are few and the mixed model’s random-effect structure is hard to estimate. Table 12.5 summarizes the comparison. The two methods are coordinated, not opposed: the generalized linear mixed model of Chapter 15 is the subject-specific counterpart of the marginal model developed here, and the attenuation relation of Section 12.1 is exactly the bridge between their coefficients.

Note. The estimand governs the choice: a population-average question leads to the generalized estimating equation, a subject-specific question to the generalized linear mixed model. Missingness beyond missing-completely-at-random requires weighting or a likelihood method; the need for random-effect variances requires the mixed model. For a linear outcome the two estimands coincide.

12.7 Running Generalized Estimating Equations in R

The geepack package fits generalized estimating equations, taking the cluster identifier through its id argument and the working correlation through corstr, and reporting robust standard errors by default. The population-average remission model and the comparison across working correlations are a few lines.

library(geepack); library(dplyr)
rct <- readRDS("Examples/data/therapy_rct.rds") |>
  mutate(arm = factor(arm, levels = c("Control","Treatment")),
         remit = as.integer(hdrs <= 7)) |>
  filter(!is.na(hdrs)) |> arrange(patient_id, week)

# --- Population-average logistic model; robust SEs are the default ---
fit <- geeglm(remit ~ week * arm, id = patient_id, data = rct,
              family = binomial, corstr = "exchangeable")
summary(fit)                                   # week:armTreatment is the PA effect

# --- Coefficients are consistent across working correlations ---
sapply(c("independence","exchangeable","ar1"), function(cs)
  coef(geeglm(remit ~ week * arm, id = patient_id, data = rct,
              family = binomial, corstr = cs))["week:armTreatment"])

The population-average probabilities for the margins plot come from predict on the fitted model, and the inverse-probability-weighted analysis passes weights to the same fitting function. The complete code, including the attenuation and coverage simulations and the missing-at-random fragility exhibit with its cumulative dropout weights, is the shipped script ch12_analysis_V01.R, with figures drawn by ch12_figures_V01.R.

# --- Population-average predicted probabilities by arm and week ---
newd <- expand.grid(week = 0:11, arm = factor(c("Control","Treatment")))
newd$p <- predict(fit, newdata = newd, type = "response")

# --- Inverse-probability weighting to relax MCAR toward MAR ---
# weights w = 1 / (cumulative probability of remaining), from a dropout model;
# pass to geeglm and use an independence working correlation
gee_ipw <- geeglm(y ~ factor(t), id = id, data = observed, weights = w,
                  corstr = "independence")

12.8 Interpreting and Reporting the Results

A generalized-estimating-equation result is reported by naming the estimand, describing the three model pieces, and interpreting effects on the scale of interest. A model paragraph for the therapy trial reads: “The population-average probability of remission was modeled with a logistic generalized estimating equation with an exchangeable working correlation and robust standard errors, treating patients as clusters (\(N = 240\)). The arm-by-week interaction, the population-average treatment effect, was positive on the log-odds scale (\(\hat\beta = 0.13\), robust \(SE = 0.09\)), and the model-implied population-average remission probability rose faster in the treatment arm, reaching a higher rate by week 12 (Figure 12.6). Estimates were essentially unchanged across independence, exchangeable, and autoregressive working correlations, as expected given the consistency of the coefficients. Because dropout was plausibly missing at random rather than completely at random, the primary inference relied on the mixed model of Chapter 13, with the generalized estimating equation reported as the population-average complement.” The estimand is named, the robust standard errors and working correlation are stated, and the missingness caveat is honest.

12.9 Common Misconceptions

Several beliefs about generalized estimating equations mislead. The first is that a generalized estimating equation is just a mixed model fitted by different software; the two answer different questions, the marginal and the conditional, which for nonlinear links are different quantities, and are valid under different missingness assumptions. The second is that robust standard errors fix everything; the sandwich estimator protects against a misspecified working correlation but not against a misspecified mean model, and it is itself fragile when clusters are few. The third is that exchangeable correlation is a safe default for longitudinal data; longitudinal correlations decay with time, so an autoregressive structure is usually more plausible, though the choice affects only efficiency. The fourth, and the most consequential, is that the method handles missing data well because it makes no distributional assumptions; the reverse is true, as the unweighted method requires the strong missing-completely-at-random condition and is biased under the missing-at-random dropout that likelihood methods handle. A recurring question, why a reviewer might demand one method or the other, is answered by returning to the estimand: name the quantity the study is about, and the method follows.

Where to Go Next

The marginal model is one of two roads out of the classical methods, and the mixed model is the other, developed across Part IV. Chapter 13 introduces the linear mixed model whose random intercept and slope restore the individual differences the marginal model omits, and whose likelihood estimation is valid under missing at random; the robust standard errors learned here reappear there as an option for clustered inference. Chapter 15 completes the pair begun in this chapter, developing the generalized linear mixed model as the subject-specific counterpart to the marginal model, where the attenuation relation of Section 12.1 becomes the explicit translation between the two sets of coefficients. The inverse-probability-weighting logic returns in the survival models of Chapter 29, and the estimand-first discipline, naming the quantity before choosing the model, recurs wherever a causal or descriptive target must be made precise, most pointedly in the cross-lagged panel debates of Chapter 21. The lesson carried forward is that the question chooses the estimand, the estimand chooses the model, and no amount of robustness substitutes for asking the question correctly.

Exercises

1. 12.1 Verify the attenuation. Simulate person-specific logistic data with a known random-intercept variance, fit both a generalized linear mixed model and a generalized estimating equation, and confirm that the ratio of their slopes matches the attenuation formula.
2. 12.2 Working correlations. Refit the remission example under independence, exchangeable, and autoregressive working correlations, and report how the coefficient and its robust standard error change; explain the stability of the coefficient.
3. 12.3 Build the fragility. Construct a missing-at-random dropout mechanism on a simulated dataset, show the bias of the unweighted generalized estimating equation against the mixed model, and repair it with inverse-probability weighting.
4. 12.4 Small-cluster coverage. Simulate a study with fifteen clusters and estimate the empirical coverage of the naive and a bias-corrected sandwich standard error.
5. 12.5 Estimand essay. For three research vignettes, decide whether the population-average or the subject-specific estimand is wanted, and defend each choice in a short paragraph.

References