GEE Population-Average (Marginal) Models

Generalized estimating equations (GEE)

extend generalized linear models to correlated data — repeated measures within a patient, encounters within a provider, or members within a household — by modeling the marginal mean of the outcome as a function of covariates while treating the within-cluster correlation as a nuisance. You specify

1. a link and variance function for the marginal mean (e.g., logit + binomial, log + Poisson, identity + Gaussian),
2. a working correlation structure (independence, exchangeable, AR(1), unstructured), and
3. the robust ("sandwich") variance estimator.

The estimating equations are solved without a full likelihood, so GEE is semiparametric: the point estimate of the regression coefficients is consistent for the marginal mean model even if the working correlation is wrong, and the sandwich variance delivers valid standard errors as the number of clusters grows. This makes GEE the workhorse for population-average questions in claims, EHR, and registry data, where each beneficiary contributes many correlated observations and you want a single contrast that describes the average effect across the population.

Core conceptual distinction (marginal vs conditional estimand)

This is the line every reviewer checks first. GEE estimates a population-average (marginal) effect: how the average outcome in the population shifts when a covariate changes, averaging over all clusters. A generalized linear mixed model (GLMM, e.g., a random-intercept logistic model) estimates a subject-specific (conditional) effect: how an individual cluster's outcome shifts, holding its random effect fixed.

For linear (identity) and log links these two targets coincide (the link is collapsible), so a GEE rate ratio and a mixed-model rate ratio estimate the same number. For the logit link they do not: because the logistic function is non-collapsible, the population-average log-odds-ratio is attenuated toward the null relative to the subject-specific one by an approximate shrinkage factor of (c^2 σ_b^2 + 1)^(-1/2) with c = 16√3/(15π) ≈ 0.588, so c^2 ≈ 0.346 (Zeger, Liang & Albert 1988). A 0.7 subject-specific log-OR with a large random-intercept variance can become a materially smaller marginal log-OR — same data, different estimand, both correct for their question.

Three further pillars: (i) consistency vs efficiency — the working correlation choice does not bias the marginal coefficients, it only affects efficiency; the sandwich is what protects you, so an honest default is "independence working correlation + sandwich SE" unless you have a strong, exogenous reason to model the correlation. (ii) Missing-data regime — ordinary GEE is consistent only under MCAR (missing completely at random); under MAR you must use weighted GEE (inverse-probability-of-observation weights) or move to a likelihood-based mixed model, which is valid under MAR.

(iii) Time-varying covariate feedback — the Pepe–Anderson (1994) caveat: when a time-varying covariate is endogenous (affected by prior outcome, i.e., there is feedback), GEE with a non-independence working correlation is biased for the marginal effect; the independence working correlation with the sandwich variance is the safe choice, and for time-varying confounders on the causal pathway you should leave GEE entirely for marginal structural models / g-methods.

Interpreting the output

Consider a Poisson GEE fit to quarterly hospitalization counts (SGLT2 vs DPP-4 initiators), with exchangeable working correlation and the sandwich variance. The model returns an adjusted hospitalization rate ratio of 0.82 (95% CI 0.71–0.94) for the SGLT2 arm.

Formal interpretation: The GEE coefficient exp(beta_arm) = 0.82 is a population-average (marginal) rate ratio — it compares the average hospitalization rate across all SGLT2 initiators to the average rate across all DPP-4 initiators, averaging over the distribution of random effects in the population.

This is explicitly not the same question as a mixed model: if a GLMM were fit to the same data, its conditional log rate ratio would be larger in magnitude (further from 1) than this marginal 0.82, because the logistic/log link is non-collapsible and the subject-specific effect for any given patient is more extreme than the population average. Both numbers are correct for their respective estimands; they are not interchangeable.

Practical interpretation: On average across the study population, SGLT2 initiators were hospitalized at 18% lower rate per quarter than DPP-4 initiators after covariate adjustment. The sandwich variance is what makes the 95% CI honest despite the working correlation structure being only approximately correct. For policy, budget-impact, and HTA purposes — "what is the effect on hospitalization burden if we switch this population?" — the marginal GEE rate ratio is the directly relevant quantity. For predicting an individual patient's trajectory, a mixed model is the appropriate tool.

Pros, cons, and trade-offs

(specific and comparative).

• vs GLMM / mixed-effects models (random intercepts/slopes): GEE targets the marginal estimand, makes no distributional assumption on cluster effects, is robust to correlation misspecification, and is computationally trivial. Cost: it gives no cluster-level (subject-specific) inference, is less efficient than a correctly specified GLMM, and — critically — is only MCAR-valid, whereas a mixed model is MAR-valid and recovers individual trajectories. Prefer GEE when the policy/HTA question is "what is the average effect in the population" (e.g., average reduction in hospitalizations under drug A vs B); prefer a GLMM when you need subject-specific prediction, have informative dropout (MAR), or want to model heterogeneity in slopes.
• vs MMRM (mixed model for repeated measures): MMRM is a Gaussian GLMM with an unstructured time covariance, the regulatory default for continuous longitudinal endpoints in RCTs, and MAR-valid. GEE is more flexible across non-Gaussian outcomes (binary, count) and gives a marginal interpretation, but is anticonservative with few clusters and MCAR-limited. Prefer MMRM for a continuous endpoint with monotone MAR dropout; prefer GEE for binary/count marginal effects.
• vs cluster-robust ("sandwich") SE on an ordinary GLM: A GLM fit independence + cluster-robust SE is algebraically the special case of GEE with an independence working correlation. GEE adds value only when you model a non-trivial working correlation for efficiency; otherwise the two are the same estimator. Prefer the explicit GEE machinery when you want exchangeable/AR(1) efficiency gains or weighted (IPW) GEE for MAR.
• vs marginal structural models (MSM) / g-methods: Both target marginal effects, but MSMs handle time-varying confounding affected by prior treatment via inverse-probability-of-treatment weighting — exactly the setting where plain GEE is biased. Never use plain GEE for a sustained treatment effect with feedback; use an IPTW-weighted GEE (which is the MSM fitting step) instead.

When NOT to use — and when it is actively misleading or dangerous

• Few clusters (< ~30–40). The sandwich variance is consistent only as the number of clusters grows; with few clusters it is anticonservative (SEs too small, type-I error inflated). Apply a bias-corrected sandwich (Mancl–DeRouen 2001 or Kauermann–Carroll) or a small-sample df correction — do not report the naive robust SE.
• Endogenous time-varying covariate with feedback. If a covariate is affected by prior outcome (e.g., dose titrated in response to last period's value), GEE with exchangeable/AR(1) is biased (Pepe–Anderson). Use independence + sandwich at minimum, or move to an MSM/g-method if the covariate is a time-varying confounder on the causal pathway.
• Informative cluster size. When the number of observations per cluster is associated with the outcome (e.g., sicker patients generate more encounters), standard GEE silently up-weights large clusters and estimates a cluster-size-weighted marginal effect, not the per-subject one; cluster- weighted GEE (CWGEE) is required to recover the intended estimand.
• Non-ignorable (MNAR) dropout. Neither GEE nor weighted GEE rescues you; sensitivity analysis under explicit MNAR assumptions is needed.
• The actual question is subject-specific. If you want an individual patient's predicted trajectory or to quantify between-patient heterogeneity, a marginal model answers the wrong question — fit a GLMM.

Data-source operational depth

(claims vs EHR vs registry vs linked).

• Claims (FFS vs MA): Repeated observations are naturally generated by period-level aggregation (e.g., quarterly hospitalization counts per `person_id`). The dominant failure mode is that cluster size is itself a function of observability: Medicare Advantage person-time lacks adjudicated FFS claims, so an MA-only quarter shows zero events not because none occurred but because they are unobserved — pool MA and FFS person-time and you bias the marginal rate downward differentially by plan. Restrict the analytic window to FFS-observable person-time (Parts A/B, or commercial medical+pharmacy) and treat plan-switching as censoring, not a true zero. Continuous-enrollment and washout rules must be applied per analysis period, and the correlation across a patient's periods is exactly what GEE handles — but only if the periods are real (an unenrolled gap is missing, not a correlated zero). Watch immortal time when the cluster is built around a procedure: person-time before the qualifying procedure is event-free by construction.
• EHR: Repeated measures are encounter-driven, so observation times are themselves informative — sicker patients are measured more often, inducing informative cluster size and MAR/MNAR observation. Model the visit process (IPW-GEE with weights for being observed) or you get a measured- when-sick bias. External-care leakage means a patient's "improvement" may be a visit elsewhere; site workflow drives systematic missingness by clinic.
• Registry: Scheduled visits give cleaner balanced repeated measures and adjudicated outcomes, but enrollment is selective and follow-up completeness varies by site — report per-wave completeness and consider weighting to the source population for transportability.
• Linked claims–EHR–registry: Best substrate (EHR/registry severity + claims completeness for the event count denominator), but linkage selects the linkable subset and the period boundaries from the three sources must be reconciled before clusters are defined; date discrepancies create spurious within-cluster correlation if a single real event is counted in two adjacent periods.

Worked claims example

Question: does initiating an SGLT2 inhibitor vs a DPP-4 inhibitor change the population-average rate of all-cause hospitalization over the first year among new initiators with type 2 diabetes in a 100% Medicare FFS sample?

1. Build the cohort with an active-comparator, new-user design (365-day continuous A/B enrollment washout, first qualifying fill = `index_date`, arm from the dispensed NDC).
2. Split each patient's first year into four 90-day periods, keeping only periods fully covered by FFS enrollment (drop MA-only and post-disenrollment periods so a zero is a true zero, not unobserved).
3. For each `person_id` × period, count inpatient stays from facility claims by revenue/bill-type, and record `offset = log(person-days observed in the period)` to handle partial periods.
4. Fit a marginal Poisson GEE: `hosp_count ~ arm + period + arm:period + age + sex + baseline_comorbidity`, link = log, distribution = Poisson, subject = person_id, working correlation = exchangeable (periods within a patient are positively correlated), with the sandwich variance. The `exp(arm)` coefficient is the population-average hospitalization rate ratio; the `arm:period` terms test whether the marginal effect changes over the year.
5. Because periods within a patient with informative dropout violate MCAR, refit with inverse-probability-of-observation weights (weighted GEE) for the disenrollment process and compare.
6. Because the cluster count is large but cluster size (number of observed periods) is associated with frailty, run a cluster-weighted GEE sensitivity analysis; and because the working-correlation choice is a nuisance, confirm the exchangeable and independence fits give the same point estimate (they should — only the SE differs).

Decision diagram

flowchart TD
  Q[Correlated / clustered outcome:<br/>repeated measures, recurrent events, clustered cross-section] --> E{What is the estimand?}
  E -->|Population-average<br/>marginal effect| F{Time-varying covariate<br/>affected by prior outcome?}
  E -->|Subject-specific effect,<br/>individual trajectory, heterogeneity| GLMM[Generalized linear mixed model<br/>PROC GLIMMIX / lme4 / MixedLM]
  F -->|No feedback| GEE[GEE marginal model<br/>working corr + sandwich SE]
  F -->|Yes: time-varying confounding| MSM[Marginal structural model<br/>IPTW-weighted GEE / g-methods]
  GEE --> M{Missingness regime?}
  M -->|MCAR| OK[Unweighted GEE valid]
  M -->|MAR| IPW[Weighted GEE<br/>inverse-prob-of-observation]
  M -->|MNAR| SENS[MNAR sensitivity analysis]
  GEE --> C{Number of clusters?}
  C -->|Large| ROB[Naive sandwich SE OK]
  C -->|Small < ~40| BC[Bias-corrected sandwich<br/>Mancl-DeRouen / Kauermann-Carroll]

flowchart LR
  subgraph IND[Independence]
    i1["1 0 0<br/>0 1 0<br/>0 0 1"]
  end
  subgraph EXC[Exchangeable]
    e1["1 r r<br/>r 1 r<br/>r r 1"]
  end
  subgraph AR1["AR(1)"]
    a1["1 r r^2<br/>r 1 r<br/>r^2 r 1"]
  end
  subgraph UNS[Unstructured]
    u1["1 r12 r13<br/>r12 1 r23<br/>r13 r23 1"]
  end
  IND -->|safe default;<br/>robust to feedback| USE[Working correlation choice<br/>affects efficiency, not consistency]
  EXC -->|unordered clusters,<br/>members within site| USE
  AR1 -->|equally spaced<br/>repeated measures| USE
  UNS -->|few balanced<br/>time points| USE

Worked example

Scenario

A researcher wants to know whether patients with type 2 diabetes who start Drug A have fewer hospitalizations per quarter than those who start Drug B, on average across the population. Four patients each contribute up to four quarterly observations. Because the same patient appears in multiple rows, those rows are correlated and ordinary regression would understate the uncertainty. GEE accounts for the within-person correlation and returns a single population-average rate comparison.

Dataset

Quarterly hospitalization counts per patient. Each row is one patient-quarter. Two patients per arm, two quarters each (a minimal illustration).

Steps

Result

GEE population-average rate ratio: Drug A vs Drug B = 0.25 / 1.25 = 0.20. Interpreted: on average across the population, patients on Drug A have about one-fifth the quarterly hospitalization rate of patients on Drug B. This is a population-level statement, not a prediction for any specific patient.

Trade-offs

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.genmod.cov_struct import Exchangeable, Independence

def fit_pa_poisson_gee(panel: pd.DataFrame):
    panel = panel.sort_values(["person_id", "period"]).copy()
    panel["log_obs_days"] = np.log(panel["obs_days"])  # offset for partial/uneven periods

    # Marginal Poisson GEE: subject = person_id (the cluster), exchangeable working correlation.
    # The robust (sandwich) covariance is statsmodels' default for GEE.
    model = smf.gee(
        formula="hosp_count ~ C(arm, Treatment('DPP4')) + C(period) + age + sex + baseline_comorbidity",
        groups="person_id",
        data=panel,
        offset=panel["log_obs_days"].values,
        cov_struct=Exchangeable(),
        family=sm.families.Poisson(),
    )
    res = model.fit()  # res.cov_type is 'robust' (sandwich) by default

    # Population-average rate ratio for the active arm vs the comparator, with sandwich CI.
    params, ci = res.params, res.conf_int()
    arm_term = [t for t in params.index if t.startswith("C(arm")][0]
    rr = np.exp(params[arm_term])
    rr_lo, rr_hi = np.exp(ci.loc[arm_term, 0]), np.exp(ci.loc[arm_term, 1])

    # Sensitivity: independence working correlation should give the SAME point estimate (only SE differs).
    rr_indep = np.exp(model.__class__(
        model.endog, model.exog, groups=panel["person_id"],
        offset=panel["log_obs_days"].values,
        cov_struct=Independence(), family=sm.families.Poisson(),
    ).fit().params[list(params.index).index(arm_term)])

    return {"rate_ratio": rr, "ci95": (rr_lo, rr_hi),
            "rate_ratio_independence": rr_indep, "summary": res.summary()}

library(geepack)

fit_pa_poisson_gee <- function(panel) {
  panel <- panel[order(panel$person_id, panel$period), ]
  panel$arm <- relevel(factor(panel$arm), ref = "DPP4")

  # Marginal Poisson GEE: id = cluster, exchangeable working correlation, log offset for exposure time.
  fit <- geeglm(
    hosp_count ~ arm + factor(period) + age + sex + baseline_comorbidity,
    id          = person_id,
    data        = panel,
    family      = poisson(link = "log"),
    corstr      = "exchangeable",
    offset      = log(panel$obs_days),
    std.err     = "san.se"            # robust sandwich SE (the GEE default)
  )

  # Population-average rate ratio for the active arm and its Wald 95% CI from the sandwich SE.
  s   <- summary(fit)$coefficients
  arm <- grep("^arm", rownames(s), value = TRUE)[1]
  est <- s[arm, "Estimate"]; se <- s[arm, "Std.err"]
  rr  <- exp(est); ci <- exp(est + c(-1.96, 1.96) * se)

  # Robustness check: independence working correlation gives the same point estimate.
  fit_indep <- update(fit, corstr = "independence")
  list(rate_ratio = rr, ci95 = ci,
       rate_ratio_independence = exp(coef(fit_indep)[arm]), fit = fit)
}

/* Log person-time offset for partial/uneven periods. */
data work.panel;
  set work.panel;
  log_obs_days = log(obs_days);
run;

/* GENMOD must read the data sorted by the subject (cluster) id for the REPEATED statement. */
proc sort data=work.panel; by person_id period; run;

/* Marginal (population-average) Poisson GEE. The REPEATED statement is what makes GENMOD fit GEE. */
proc genmod data=work.panel;
  class person_id arm(ref='DPP4') sex period / param=ref;
  model hosp_count = arm period age sex baseline_comorbidity
        / dist=poisson link=log offset=log_obs_days;
  repeated subject=person_id / type=exch corrw modelse;  /* exchangeable working corr; print empirical + model SE */
  estimate 'PA rate ratio SGLT2 vs DPP4' arm 1 / exp;     /* exp() gives the population-average rate ratio */
run;

/* Robustness: independence working correlation should reproduce the point estimate (only SE changes). */
proc genmod data=work.panel;
  class person_id arm(ref='DPP4') sex period / param=ref;
  model hosp_count = arm period age sex baseline_comorbidity
        / dist=poisson link=log offset=log_obs_days;
  repeated subject=person_id / type=ind;
  estimate 'PA rate ratio (independence)' arm 1 / exp;
run;

Citations

1. [1]Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13-22.
2. [2]Zeger SL, Liang KY, Albert PS. Models for longitudinal data: a generalized estimating equation approach. Biometrics. 1988;44(4):1049-1060.
3. [3]Hubbard AE, Ahern J, Fleischer NL, Van der Laan M, Lippman SA, Jewell N, Bruckner T, Satariano WA. To GEE or not to GEE: comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology. 2010;21(4):467-474.

1. [4]Pepe MS, Anderson GL. A cautionary note on inference for marginal regression models with longitudinal data and general correlated response data. Communications in Statistics - Simulation and Computation. 1994;23(4):939-951.

1. [5]Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126-134.