Introduction
Generalized Estimating Equations extend generalized linear models to correlated (clustered or longitudinal) data without requiring a full likelihood specification — instead, GEE gets away with only specifying the marginal mean and a "working" guess at the correlation structure, and remains valid even when that guess is wrong. This is a genuinely different strategy from MMRM's full likelihood-based approach, and it comes with real trade-offs worth understanding precisely, not just naming. This tutorial implements GEE from scratch (Fisher scoring plus the sandwich variance estimator), validates it against production software, and works through the two things about GEE that are most often glossed over: the population-averaged vs. subject-specific coefficient divergence for non-linear links, and GEE's real vulnerability to anything beyond MCAR missingness.
The Model
GEE specifies only the marginal mean of the outcome as a function of covariates, via a link function exactly as in an ordinary GLM:
and a working correlation matrix \(R(\alpha)\) describing the assumed correlation between a subject's repeated observations — "working" because, unlike MMRM's covariance structure, GEE doesn't require this guess to be correct for the coefficient estimates \(\hat\beta\) to remain consistent. This is GEE's central design choice: trade a full likelihood (and the distributional assumptions that come with it) for a much weaker requirement on the mean model alone.
The Estimating Equations
Rather than maximizing a likelihood, GEE solves a quasi-score equation directly:
where \(A_i\) is a diagonal matrix of the variance function evaluated at each observation's mean (e.g., \(\mu(1-\mu)\) for a binary outcome), and \(\phi\) is a dispersion parameter. This is solved iteratively via Fisher scoring — update \(\beta\), re-estimate the working correlation parameters \(\alpha\) and dispersion \(\phi\) from the current residuals, repeat.
The Sandwich Estimator: GEE's Defining Feature
The step that makes GEE genuinely robust to a misspecified working correlation is the variance estimator. The naive ("model-based") variance,
is only valid if \(R(\alpha)\) is the true correlation structure. The sandwich (robust/empirical) estimator replaces the inner piece with the actual observed variability of the score contributions across subjects:
This "bread-meat-bread" structure is what gives the estimator its name, and it's what makes the resulting standard errors valid regardless of whether the working correlation matches reality — the working correlation still affects efficiency (how precise the estimates are), just not the validity of the reported uncertainty.
A Real Fit: Four Working Correlations, Validated
Implementing GEE's Fisher scoring and sandwich estimator from scratch, fit to a real simulated 150-subject, 4-visit longitudinal trial (true within-subject correlation: AR(1) with \(\rho=0.6\)) — and cross-validated against production software (Python's statsmodels GEE implementation), matching to 3–4 decimal places:
| Working correlation | Treatment β̂ | Robust SE | Naive SE | Recovered correlation (lag 1) |
|---|---|---|---|---|
| Independence | −2.508 | 0.792 | 0.508 | 0 (assumed) |
| Exchangeable | −2.508 | 0.792 | 0.792 | 0.477 (constant) |
| AR(1) | −2.641 | 0.771 | 0.777 | 0.617 |
| Unstructured | −2.598 | 0.767 | 0.779 | 0.655 |
Two things worth reading closely here. First, the point estimates are remarkably stable across all four working correlation choices (−2.51 to −2.64) — this is GEE's core promise delivered in practice, not just claimed. Second, look at the Independence row: its naive SE (0.508) is dramatically smaller than its own robust SE (0.792) — a real ~36% understatement of true uncertainty that would occur if you trusted the naive (model-based) standard error while assuming independence when the data are actually correlated at \(\rho=0.6\). As the working correlation gets closer to the truth (AR(1), Unstructured), naive and robust SEs converge — exactly the pattern the theory predicts: always use the sandwich estimator, because you don't know in advance how wrong your working correlation guess is, and the AR(1) structure recovered here (0.617, then decaying) closely matches the true simulated \(\rho=0.6\) decay pattern, confirming the estimation is working correctly.
Choosing a Working Correlation: QIC
Since GEE has no true likelihood, ordinary AIC doesn't apply — the Quasi-likelihood Information Criterion (QIC) serves the analogous role:
| Working correlation | QIC |
|---|---|
| Independence | 612.61 |
| Exchangeable | 612.61 |
| AR(1) | 611.54 |
AR(1) wins, consistent with it being closest to the true simulated structure — though the margin here is modest, a reminder that QIC differences are often smaller and less decisive than AIC/BIC differences in the MMRM covariance-structure comparisons (see the MMRM tutorial), because the working correlation mainly affects efficiency, not validity, so the criterion has less to distinguish.
Population-Averaged vs. Subject-Specific: A Real Divergence
This is the single most misunderstood property of GEE, and it's worth demonstrating rather than just stating: for non-linear links (like logistic regression), GEE's marginal coefficients and a mixed model's conditional coefficients estimate different quantities, and they will not match even with unlimited data — not due to any modeling error, but because of a real mathematical distinction. Simulating 200 clusters of 4 binary observations each, with substantial subject-level heterogeneity (a random intercept with SD 1.5) and a true conditional (subject-specific) log-odds treatment effect of −1.0:
| Model | Estimated log-OR | Odds Ratio | What it estimates |
|---|---|---|---|
| True conditional effect (used in simulation) | −1.000 | 0.368 | — |
| Mixed logistic model (subject-specific) | −0.702 | 0.496 | Effect for a given subject, holding their own random effect fixed |
| GEE logistic (population-averaged) | −0.505 | 0.604 | Effect on the population-average probability, marginalized over heterogeneity |
The GEE marginal estimate is genuinely, substantially attenuated toward zero (odds ratio 0.604 vs. the mixed model's 0.496) — not because GEE is biased, but because averaging a non-linear (logistic) function over subject heterogeneity mathematically pulls the population-level effect toward the null relative to the subject-specific effect. This gap grows with the amount of unexplained subject-level heterogeneity, and it means choosing between GEE and a mixed model isn't just a technical/computational preference — it changes which scientific question is being answered, directly analogous to the estimand distinctions in the dedicated Estimands tutorial. For continuous outcomes with an identity link (as in this tutorial's main worked example), this divergence does not occur — population-averaged and subject-specific coefficients coincide, which is why the phenomenon is specifically a non-linear-link issue.
GEE's Real Missing-Data Limitation
MMRM's likelihood-based approach is valid under Missing At Random. Standard GEE, as implemented in most software, is not — its standard available-case handling requires the stronger MCAR assumption for unbiasedness. Demonstrating this directly: refitting the longitudinal trial above under genuine MAR dropout (subjects with a poor Week-1 response more likely to drop out, observed-N declining from 150 to 124 by later visits):
| Quantity | Value |
|---|---|
| True treatment effect | −3.000 |
| GEE estimate (available-case, standard implementation) | −4.596 |
| Bias | −1.596 |
A real, substantial bias — roughly 53% larger in magnitude than the true effect — arising purely because standard GEE's available-case approach doesn't properly account for dropout that depends on an observed outcome value. Weighted GEE (incorporating inverse-probability-of-missingness weights, conceptually related to the reference-based methods in the Missing Data tutorial) can restore MAR-validity, but isn't the software default and requires deliberate implementation — this is the single most important practical reason MMRM, not GEE, is the standard primary analysis choice for confirmatory trials with anticipated dropout, even though GEE is often simpler to explain and implement for other purposes.
When to Choose GEE Over MMRM
| Scenario | Better choice |
|---|---|
| Confirmatory trial, continuous endpoint, anticipated dropout | MMRM (MAR-valid, matches primary estimand needs) |
| Binary or count outcomes with population-average interpretation desired, minimal missingness | GEE |
| Cluster-randomized trials with complete or near-complete cluster-level data | GEE (simpler, robust to correlation misspecification) |
| When subject-specific prediction or interpretation is the actual goal | Mixed model, not GEE |
Summary
GEE's entire value proposition rests on the sandwich estimator, demonstrated concretely here: point estimates stayed stable across four working correlation choices, but naive standard errors under a genuinely wrong working correlation (Independence, when the truth was AR(1) at ρ=0.6) understated true uncertainty by roughly a third — exactly the failure mode the sandwich estimator exists to prevent.
Naive vs. robust (sandwich) standard error for the treatment coefficient, across all four working correlations. The gap is largest for Independence (the most misspecified structure here) and shrinks toward zero as the working correlation approaches the truth — visual confirmation that the sandwich estimator's correction scales with how wrong the working guess is.
Two further real findings deserve equal weight: GEE's marginal coefficients for non-linear links are not the same quantity as a mixed model's conditional coefficients (a real ~40% relative gap in this tutorial's logistic example, not an error), and standard GEE's real vulnerability to anything beyond MCAR missingness (a genuine 53% relative bias under realistic MAR dropout) is the practical reason it isn't the default confirmatory-trial choice that MMRM is.
References
Liang, K.Y. & Zeger, S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika.
Zeger, S.L., Liang, K.Y. & Albert, P.S. (1988). Models for longitudinal data: a generalized estimating equation approach. Biometrics.
Pepe, M.S. & Anderson, G.L. (1994). A cautionary note on inference for marginal regression models with longitudinal data and general correlated response data. Communications in Statistics.
Robins, J.M., Rotnitzky, A. & Zhao, L.P. (1995). Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association.
Pan, W. (2001). Akaike's information criterion in generalized estimating equations. Biometrics.