GEE Population-Average (Marginal) Models

In plain language

GEE (Generalized Estimating Equations) is a statistical method for answering the question: on average across the whole population, how does an outcome differ between two groups, when each person contributes multiple measurements over time? Instead of treating each measurement as if it came from a different person, GEE acknowledges that readings from the same person tend to be similar to each other, and it accounts for that similarity when calculating the result. The answer it gives describes the average shift in the population, not any single patient's trajectory, which is exactly what payers, regulators, and health technology assessors usually need.

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 to use

(decision rules). Population-average comparative effectiveness, safety, or utilization questions on correlated outcomes: repeated lab values (HbA1c), recurrent events / counts (hospitalizations, ED visits) modeled per period, repeated binary status (controlled vs uncontrolled), clustered cross- sections (patients within practices), or paired/eye-level/limb-level data. Use GEE when the estimand is the average effect, the number of independent clusters is reasonably large (rule of thumb ≥ ~40), and missingness is plausibly MCAR or you can build observation weights for MAR. It is the natural outcome-model engine for a population-average treatment contrast after propensity-score weighting.

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).

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

• Identify the clusters: person_id is the cluster. Rows sharing a person_id are correlated because they come from the same person.

• Calculate the raw per-arm average to get an intuition: Drug A patients had (0+0+1+0) = 1 total hospitalization across 4 patient-quarters, an average of 0.25 per quarter. Drug B patients had (1+1+2+1) = 5 total, an average of 1.25 per quarter.

• GEE is told that observations within the same person are likely correlated (exchangeable working correlation: any two quarters from the same person are assumed equally similar). This is the analyst's modeling choice for efficiency, not a requirement for validity.

• GEE solves its estimating equations across all clusters simultaneously, using the sandwich variance so that the standard errors are valid regardless of whether the exchangeable assumption is exactly right.

• The result is a population-average rate ratio: the average quarterly hospitalization rate in Drug A patients divided by the average rate in Drug B patients, interpreted as the contrast you would see if you could intervene on the whole population.

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.

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)
}