statminds

Marginal GLM (Correlated-Outcome Model)Parametric Reference12-stage workflow

GEE (Generalized Estimating Equations)

The engine for Population-Averaged Discovery. GEE audits clustered and longitudinal data by focusing on the 'Global Average' effect while utilizing robust standard errors to neutralize internal correlations.

Model familyMarginal GLM (Correlated-Outcome Model)

Hypothesistwo-tailed (can be one-tailed for directional hypotheses)

AliasesPopulation-Averaged Model · Marginal Model · Quasi-Likelihood Estimation

Population-Average Audit

Determine how predictors shift the entire population mean, regardless of individual subject trajectories.

Correlation Structure Forensics

Model the 'Internal Web' of data (e.g., Autoregressive or Exchangeable) to ensure efficient precision.

Robust Stability Discovery

Achieve consistent parameter estimates even if the exact distribution of random effects is unknown or messy.

Visual Overview Dashboard

What is it?

GEE (Generalized Estimating Equations) is designed to analyze clustered, longitudinal, or repeated measures data by modeling both population trends and correlation structures.

Goals & Indications

Population-Average Audit: Determine how predictors shift the entire population mean, regardless of individual subject trajectories.
Correlation Structure Forensics: Model the 'Internal Web' of data (e.g., Autoregressive or Exchangeable) to ensure efficient precision.
Robust Stability Discovery: Achieve consistent parameter estimates even if the exact distribution of random effects is unknown or messy.

Core Idea Diagram

Hypotheses

H₀: H₀: β₁ = 0 (no population-average effect of predictor on outcome)

Hₐ: Hₐ: β₁ ≠ 0 (predictor has population-average effect on outcome)

How it works

Specify the link function and marginal population-average model.
Select a working correlation structure (Independent, Exchangeable, AR1).
Solve Generalized Estimating Equations to obtain regression coefficients.
Apply the Huber-White 'sandwich' estimator to obtain robust standard errors.

Assumptions

Outcome variable follows a distribution from the exponential family: Outcome has appropriate marginal distribution

Correct specification of link function relating linear predictor to mean response: Appropriate link function for outcome type

Data are missing completely at random: Missing data unrelated to outcomes or predictors

Important Note

GEE estimates population-average (marginal) effects, not subject-specific effects. GEE uses sandwich/robust SEs that are valid even if correlation structure is misspecified. Unlike GLMMs, GEE does not model random effects; it accounts for correlation via working correlation matrix. Interpretation: average effect across population, not conditional on cluster membership.

Worked Example

Variable	Model SE	Robust SE
Treatment Effect	0.124	0.145 (Robust)
AR1 Corr (ρ)	0.528	Working Correlation

Interactive Sandbox

GEE Working Correlation Matrix

Switch between correlation structures (Independent, Exchangeable, AR1) and adjust time correlation strength ρ.

Correlation Strength (ρ)0.50

Repeated Measures (T)4

Corr(Y_t, Y_s) = ρ (Constant correlation over time)

The 12-Stage Precision Workflow

01Marginal Influence

Hypotheses

We test how the global 'Average Responder' moves across categories or time—the definitive metric for public health discovery.

02Cluster Independence

Assumptions

Ensuring that while data *within* a cluster is correlated, the clusters themselves are independent—the foundational GEE mandate.

03Working Correlation

Diagnostics

Selecting the 'Structure' (Independence, Exchangeable, AR1, Unstructured) that best captures the physical reality of your data's correlation.

04focus

Auditing the average FlowMotion success rate across 50 cities—focusing on the 'City-Level' shift rather than the individual person.

05LMM / GLMM Pivot

Alternatives

Knowing when to switch to Mixed Models if you need 'Subject-Specific' predictions rather than just population averages.

06The Robust Score

Significance

Utilizing Sandwich Estimators (Huber-White) to provide p-values that remain valid even if your correlation structure is slightly 'Wrong'.

07Population Odds/Rates

Effect Size

Interpreting parameter estimates as the shift in the global baseline—essential for summarizing broad treatment impacts.

08The 30-Cluster Shield

Sample Size

Determining the N required—ideally seeking at least 30 clusters to ensure the 'Robust' standard errors actually stabilize.

09The Marginal Narrative

Reporting

Reporting the 'Population-Averaged' results clearly to avoid the 'Ecological Fallacy' of applying global averages to specific individuals.

10geepack / gee Logic

Software

Executing 'geeglm()' commands, ensuring the 'corstr' argument aligns with your forensic audit of the data's temporal or spatial structure.

11focus

Identifying the 'MCAR Requirement'—GEE is sensitive to non-random dropouts, which can dangerously bias the marginal average.

12focus

Tracing the model back to Liang and Zeger (1986) and the foundational shift from maximum likelihood to estimating equations.

01Hypothesis test logic

Hypotheses

Pragmatic null and alternative hypotheses defined in mathematical notation.

A hypothesis is a question sharpened to a point. Ambiguity is the enemy of inference.

Logic Core

Null · H₀

H₀: β₁ = 0 (no population-average effect of predictor on outcome)

Alternative · Hₐ

Hₐ: β₁ ≠ 0 (predictor has population-average effect on outcome)

Why it matters two-tailed (can be one-tailed for directional hypotheses)

02Model diagnostics

Assumptions

The core mathematical criteria needed to ensure that statistical testing remains unbiased and valid.

Build your analysis on rock, not sand. Verify the mathematical foundation before building the model.

Integrity Shield

Assumptions

Critical / High Severity

How to check

Quick

For binary outcomes: verify 0/1 coding. For count outcomes: check non-negative integers, plot histogram to assess mean-variance relationship. For continuous outcomes: check normality of residuals. GEE is more robust to distributional misspecification than GLMMs due to robust SEs

Rigorous

For count data: compare mean to variance across clusters; if variance >> mean, consider negative binomial family or use quasi-Poisson with scale parameter. For binary: check prevalence; if rare (<5%) or very common (>95%), logistic link may have poor properties. For continuous: Q-Q plots and histogram of Pearson residuals

If violated

If count data overdispersed: specify family='NegBin' (negative binomial) or use quasi-Poisson with scale parameter estimated from data. If binary data shows extreme prevalence: use log-binomial (family=poisson, link='log') for relative risk, or use modified Poisson with robust SEs. If continuous data non-normal: GEE with robust SEs is relatively robust; consider transformations or use Gaussian family with identity link (GEE does not assume normality, only correct mean structure)

gee

How to check

Quick

Use standard links: identity for Gaussian (continuous), logit for binomial (binary), log for Poisson (count). Plot residuals vs fitted values to check for systematic patterns suggesting link misspecification

Rigorous

Compare models with different link functions using QIC (quasi-likelihood information criterion). For binomial: test logit vs probit vs cloglog. Use component+residual plots (partial residual plots) to assess linearity on link scale. Check if predicted values fall within valid range (e.g., probabilities in [0,1])

If violated

For binomial: if asymmetric response or rare events, try cloglog link. To estimate relative risk instead of odds ratio: use log link (log-binomial model). For count: if non-linear relationship on log scale, add polynomial or spline terms. For continuous: try log or square-root transformations if relationship non-linear. Use GAM to explore non-linearities, then incorporate into GEE

gee

How to check

Quick

Calculate missingness proportions by group and time. If missingness differs systematically (e.g., higher dropout in treatment vs control, or among those with worse outcomes at baseline), MCAR likely violated

Rigorous

Little's MCAR test (tests if missingness pattern is random). Compare baseline characteristics between completers vs those with missing data. Plot missingness indicators vs observed outcomes. Logistic regression: predict missingness from baseline covariates and observed outcomes (significant predictors suggest MAR/MNAR, not MCAR)

If violated

GEE with standard methods assumes MCAR. If missing at random (MAR) - missingness depends on observed covariates: (1) Use inverse probability weighting (IPW) GEE to account for differential dropout; (2) Include baseline covariates predicting missingness as covariates in model; (3) Use multiple imputation then GEE (MI-GEE). If missing not at random (MNAR) - missingness depends on unobserved values: GEE biased; consider selection models, pattern-mixture models, or sensitivity analyses. GLMMs with maximum likelihood handle MAR better than GEE (but assume random effects distribution correct)

gee

How to check

Quick

Review study design: were clusters (subjects, schools, clinics) independently sampled? Check for spatial/temporal clustering. Ensure cluster IDs are unique and observations don't belong to multiple clusters

Rigorous

If clusters have spatial locations, test for spatial autocorrelation using Moran's I on cluster-level residuals. If clusters have temporal ordering, plot residuals by time to detect trends. Check for hierarchical clustering (e.g., subjects nested in sites nested in regions) that should be modeled

If violated

If spatial correlation between clusters: (1) Include spatial covariates (region, latitude/longitude); (2) Use spatial correlation structure (not standard in GEE but available in some packages); (3) Use hierarchical models with cluster-level random effects for higher-level units. If temporal correlation between clusters: include time-period fixed effects or use GEE with nested correlation structure. If higher-level clustering: use multilevel GEE or GLMM with crossed/nested random effects

gee

How to check

Quick

Plot outcome vs continuous predictors to check linearity on link scale. Test for interactions if theoretically plausible. Add quadratic terms if relationship appears curved

Rigorous

Use generalized additive models (GAMs) with smoothers to explore non-linear relationships, then incorporate findings into GEE (e.g., add polynomial terms, splines, or categorize). Component+residual plots for each predictor. Test interactions using Wald tests or QIC model comparison

If violated

If non-linear relationships: (1) Add polynomial terms (quadratic, cubic); (2) Use restricted cubic splines (rcs package in R, then use spline basis in GEE); (3) Transform predictors (log, square root); (4) Categorize continuous predictors (use with caution - loses information). If missing interactions: test theoretically motivated interactions using Wald tests. Misspecified mean structure leads to biased estimates even with robust SEs

gee

How to check

Quick

Count number of clusters (k) and calculate mean, median, min, max cluster size. Rule of thumb: k ≥ 30-40 clusters for valid sandwich SE inference, at least 3-5 observations per cluster on average. For binary outcomes: ensure sufficient events per cluster

Rigorous

Simulation studies show sandwich SEs are anticonservative (too small) with k < 30-40, leading to inflated Type I error. Check if any clusters have < 2 observations (uninformative for within-cluster correlation). For rare outcomes: ensure at least 10 events per predictor (EPV rule)

If violated

If few clusters (k < 30): (1) Use small-sample SE corrections: Mancl-DeRouen, Fay-Graubard, or Kauermann-Carroll adjustments (available in geepack, geesmv packages in R); (2) Use bootstrap SEs with cluster resampling; (3) Use GLMM instead (better small-sample properties if random effects normally distributed); (4) Report results cautiously, noting limited precision. If small cluster sizes but many clusters: GEE still valid, but correlation structure harder to estimate; use independence or exchangeable. If very unbalanced cluster sizes: standard GEE valid but efficiency reduced

gee

How to check

Quick

For longitudinal data: try exchangeable (constant correlation), AR(1) (autoregressive), or unstructured. For clustered cross-sectional: exchangeable is typical. Use QIC (quasi-likelihood information criterion) to compare structures (lower = better). Note: misspecification affects efficiency (power), not validity (Type I error), due to robust SEs

Rigorous

Compare QIC, QICu (corrected QIC) across multiple correlation structures. Plot within-cluster correlations over time/distance to visualize structure. Examine estimated correlation parameters (e.g., AR(1) rho, exchangeable alpha). Use model-based vs robust SEs: if similar, correlation well-specified; if different, misspecification present (but robust SEs still valid)

If violated

GEE with robust SEs gives valid inference even if correlation structure misspecified (key advantage over likelihood-based methods). However, efficiency reduced (wider CIs, less power) if misspecified. To improve efficiency: (1) Try multiple structures (independence, exchangeable, AR(1), unstructured) and select using QIC; (2) For longitudinal data with equal spacing: AR(1) often reasonable; (3) For unequal spacing: use continuous-time AR; (4) If uncertain: exchangeable is conservative and robust; (5) Independence structure with robust SEs is always valid but least efficient. Sensitivity analysis: report results under multiple structures to assess robustness

03Residual Forensics

Diagnostics

Checking residual plots and indices to examine model deviations and ensure standard error integrity.

Trust, but verify. The outliers often hold more truth than the averages.

System Health

Essential checks

Check convergence (algorithm converged, no warnings)
Examine Pearson or standardized residuals vs fitted values
Check for outliers and influential clusters (Cook's distance, DFBETAS)
Compare model-based vs robust SEs (large discrepancy suggests correlation misspecification)
Plot observed vs predicted values by cluster
Calculate proportion of missing data and patterns

Recommended checks

QIC and QICu for model selection (lower = better fit)
Test multiple working correlation structures, compare QIC
Check residuals by cluster (look for systematic patterns)
Sensitivity analysis: compare results under different correlation structures
Small-sample SE corrections if k < 40 (Mancl-DeRouen, Fay-Graubard)
Bootstrap SEs with cluster resampling for confirmation
Multicollinearity check (VIF) for predictors

04Live Instances

Applied Minds

Review concrete study examples, data layout guidelines, and copy executable syntax scripts.

Theory is the map. Practice is the terrain. Simulation bridges the gap.

Applied Wisdom

Example 01

Blood Pressure Reduction in Hypertension (Continuous Longitudinal Outcome)

Research question: Does a lifestyle intervention reduce systolic blood pressure (SBP) compared to usual care over 12 months in adults with hypertension? Design: RCT with 120 participants (60 per group), measured at baseline, 3, 6, 9, 12 months (5 timepoints). Outcome: Systolic BP (mmHg, continuous). GEE with exchangeable correlation structure estimates population-average effect.

Designlongitudinal RCT

Outcome Scalesystolic BP (mmHg, continuous)

# GEE: Continuous outcome (blood pressure) with longitudinal data
# Population-average effects with exchangeable correlation

library(geepack)     # For geeglm()
library(doBy)        # For esticon() for contrasts
library(emmeans)     # For estimated marginal means
library(ggplot2)

set.seed(2025)

# Simulate realistic data
n_subjects <- 120
n_timepoints <- 5
subject_id <- rep(1:n_subjects, each=n_timepoints)
treatment <- rep(rep(c("Control", "Intervention"), each=60), each=n_timepoints)
month <- rep(c(0, 3, 6, 9, 12), times=n_subjects)

# Subject-level random effects (persistent individual differences)
subject_effects <- rnorm(n_subjects, mean=0, sd=8)
subject_effects_expanded <- rep(subject_effects, each=n_timepoints)

# Generate SBP with treatment effect
treatment_numeric <- as.numeric(treatment == "Intervention")
sbp <- 145 +                                    # Baseline SBP
       -0.8 * month +                           # Time trend in control
       -2.0 * treatment_numeric +               # Treatment main effect
       -0.25 * month * treatment_numeric +      # Treatment × time interaction
       subject_effects_expanded +               # Subject-level effect
       rnorm(length(subject_id), 0, 6)          # Measurement error

data <- data.frame(
  subject_id = factor(subject_id),
  treatment = factor(treatment, levels=c("Control", "Intervention")),
  month = month,
  sbp = sbp
)

# === STEP 1: Descriptive Statistics ===
library(dplyr)
data %>%
  group_by(treatment, month) %>%
  summarise(n = n(),
            mean_sbp = mean(sbp),
            sd_sbp = sd(sbp),
            se = sd_sbp / sqrt(n))

# === STEP 2: Fit GEE with Exchangeable Correlation ===
# Exchangeable: assumes constant correlation within subject over time
gee_model <- geeglm(sbp ~ month * treatment,
                    id = subject_id,
                    data = data,
                    family = gaussian(link="identity"),
                    corstr = "exchangeable")

summary(gee_model)

# Output:
# Coefficients:
#                        Estimate Std.err  Wald Pr(>|W|)    
# (Intercept)            145.12    1.08  18012  < 2e-16 ***
# month                   -0.79    0.08     95  < 2e-16 ***
# treatmentIntervention   -2.04    1.53      1.8   0.18    
# month:treatmentIntervention -0.24  0.12   4.0   0.046 *
#
# Estimated correlation: alpha = 0.52 (moderate within-subject correlation)

# === STEP 3: Compare Correlation Structures ===
# Try different structures to assess sensitivity

gee_indep <- geeglm(sbp ~ month * treatment, id=subject_id, data=data,
                    family=gaussian, corstr="independence")

gee_ar1 <- geeglm(sbp ~ month * treatment, id=subject_id, data=data,
                  family=gaussian, corstr="ar1")

gee_unstr <- geeglm(sbp ~ month * treatment, id=subject_id, data=data,
                    family=gaussian, corstr="unstructured")

# Compare QIC (lower = better)
QIC(gee_indep)[1]
QIC(gee_model)[1]   # Exchangeable
QIC(gee_ar1)[1]
QIC(gee_unstr)[1]
# Result: AR(1) or exchangeable typically best for equally-spaced longitudinal data

# === STEP 4: Diagnostics ===

# 4a. Check convergence
if(gee_model$converged) {
  cat("Model converged successfully\n")
}

# 4b. Residual plots
data$residuals <- residuals(gee_model, type="pearson")
data$fitted <- fitted(gee_model)

ggplot(data, aes(x=fitted, y=residuals)) +
  geom_point(alpha=0.4) +
  geom_hline(yintercept=0, linetype="dashed", color="red") +
  geom_smooth(se=FALSE, color="blue") +
  labs(title="Residuals vs Fitted Values",
       x="Fitted SBP(mmHg)", y="Pearson Residuals") +
  theme_minimal()

# 4c. Compare model-based vs robust SEs
summary(gee_model)$coefficients[, "Std.err"]      # Robust (sandwich) SEs
summary(gee_model)$coefficients[, "Std.err"]      # Model-based not directly available
# If robust >> model-based, correlation misspecified (but inference still valid)

# === STEP 5: Estimated Marginal Means ===
library(emmeans)
emm <- emmeans(gee_model, ~ treatment | month, at=list(month=c(0, 12)))
summary(emm)

# Baseline (month 0):
#   Control: 145.1 mmHg, Intervention: 143.1 mmHg (diff = -2.0 mmHg)
# Month 12:
#   Control: 135.6 mmHg, Intervention: 130.7 mmHg (diff = -4.9 mmHg)

# Contrast: treatment difference at 12 months
contrast(emmeans(gee_model, ~ treatment | month, at=list(month=12)),
         method="pairwise")
# Estimate: -4.9 mmHg (95% CI [-7.8, -2.0]), p = .001

# === STEP 6: Visualizations ===

# Plot 1: Observed means over time
data_summary <- data %>%
  group_by(treatment, month) %>%
  summarise(mean_sbp = mean(sbp),
            se = sd(sbp) / sqrt(n()))

ggplot(data_summary, aes(x=month, y=mean_sbp, color=treatment)) +
  geom_line(size=1.2) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=mean_sbp - 1.96*se, ymax=mean_sbp + 1.96*se), 
                width=0.5) +
  labs(title="Systolic Blood Pressure Over Time by Treatment",
       x="Month", y="Mean SBP ± 95% CI(mmHg)",
       color="Treatment") +
  theme_minimal()

# Plot 2: Predicted values from GEE (population-average)
newdata <- expand.grid(
  month = seq(0, 12, by=1),
  treatment = c("Control", "Intervention")
)
newdata$pred_sbp <- predict(gee_model, newdata=newdata, type="response")

ggplot() +
  geom_line(data=newdata, aes(x=month, y=pred_sbp, color=treatment), 
            size=1.2) +
  geom_point(data=data_summary, aes(x=month, y=mean_sbp, color=treatment),
             size=3, alpha=0.6) +
  labs(title="GEE Predicted SBP vs Observed Means",
       x="Month", y="Systolic BP(mmHg)",
       color="Treatment") +
  theme_minimal()

# Plot 3: Individual trajectories (sample)
sample_ids <- sample(unique(data$subject_id), 20)
data_sample <- data[data$subject_id %in% sample_ids, ]

ggplot(data_sample, aes(x=month, y=sbp, group=subject_id, color=treatment)) +
  geom_line(alpha=0.5) +
  facet_wrap(~ treatment) +
  stat_summary(aes(group=1), fun=mean, geom="line", 
               color="black", size=2) +
  labs(title="Individual SBP Trajectories(20 subjects per group)",
       x="Month", y="Systolic BP(mmHg)") +
  theme_minimal()

# === STEP 7: APA-Style Reporting ===
cat("
=== APA-Style Report ===
Generalized estimating equations(GEE) with exchangeable correlation structure
and robust standard errors were used to analyze systolic blood pressure(SBP)
over 12 months in a randomized trial comparing lifestyle intervention to usual care
(n=120, 60 per group, 5 timepoints).

Results showed a significant time × treatment interaction(Wald χ² = 4.0, p = .046),
indicating that the intervention produced greater SBP reduction over time compared
to usual care. At baseline, groups did not differ significantly(intervention: 143.1
mmHg, control: 145.1 mmHg, difference = -2.0 mmHg, p = .18). By 12 months, the
intervention group had significantly lower SBP(130.7 mmHg) compared to control
(135.6 mmHg), with a between-group difference of -4.9 mmHg(95% CI [-7.8, -2.0],
p = .001).

The estimated within-subject correlation was 0.52 (exchangeable structure),
indicating moderate correlation of repeated SBP measurements. Sensitivity analyses
using AR(1) and unstructured correlation structures produced similar results
(effect estimates within 0.3 mmHg), supporting robustness of findings.

These results support lifestyle intervention as an effective population-level
strategy for blood pressure reduction in adults with hypertension, with clinically
meaningful effects emerging by 12 months.
")

Interpretation Blueprint

Time × treatment interaction (Wald χ² = 4.0, p = .046) shows lifestyle intervention produces greater SBP reduction over time. At 12 months, intervention reduced SBP by 4.9 mmHg more than control (95% CI [-7.8, -2.0], p = .001). Exchangeable correlation (α = 0.52) indicates moderate within-subject correlation. Results robust across correlation structures (QIC similar for exchangeable and AR1). Findings align with Appel et al. (2003) JAMA showing lifestyle intervention reduces SBP by ~4 mmHg. Clinical significance: 5 mmHg reduction translates to ~10-20% reduction in cardiovascular events.

05Tactical Pivots

Alternatives

Structured fallback pathways for choosing alternative tests when normality or slopes requirements fail.

When the path is blocked, pivot. Rigor is not rigidity; it is the intelligent adaptation to reality.

Adaptive Strategy

Measurement Precision Ladder Ideal · Clustered Nominal / Continuous

Ratio

Maintain GEE. Ideal for population-average continuous effects when the exact random distribution is unknown.

Standard Signal

Binary / Count

Maintain GEE logic. The most robust engine for marginal probability and rate discovery.

Peak Precision

Ordinal

Pivot to Ordinal GEE to preserve the cumulative rank story across clusters.

Logic Leak

Temporal Trajectory Audit Longitudinal Marginal Snapshot

Population Mean

Global trajectory.

Stay with GEE. Audit how the average responder shifts while accounting for internal correlation.

Subject-Specific

Trajectory growth.

Pivot to Mixed Models (LMM/GLMM) to model individual recovery curves.

Adaptive Technical Safeguards · adaptive safeguards

non random missingness

Multiple Imputation GEE — Account for 'Missing Not At Random' data—GEE is highly sensitive to dropout bias.
Weighted GEE (WGEE) — Apply inverse-probability weights to neutralize attrition signals.

wrong correlation choice

QIC Selection Strike — Utilize the Quasi-Likelihood Information Criterion to select the optimal 'Working Correlation' structure.
Robust Standard Errors — Maintain p-value integrity even if the internal web is misspecified.

06Adjusted Comparisons

Post-hoc

Group mean comparisons and correction controls (e.g. Tukey HSD, Bonferroni) to protect against Family-Wise Error Rates.

The omnibus test opens the door; post-hoc analysis explores the room.

Forensic Detail

Adjusted Comparisons

Post-hoc pairwise tests defined for this model.

Interpretation Guidelines

All comparisons on link scale (log-odds, log-rate), then exponentiate for OR/RR. Use robust SEs for all post-hoc tests

07Standardized scale impact

Effect Size

Understanding effect sizes (e.g., Cohen's d, Partial Eta-Squared) and clinical impact benchmarks.

Significance is noise. Magnitude is the signal. Measure the impact, not just the probability.

Impact Magnitude

Population-average OR (not conditional). OR = 1.5 means 50% higher odds on average across population. For common outcomes (>10%), OR overestimates relative risk. Convert to approximate RR: RR ≈ OR / [(1 - p0) + (p0 × OR)] where p0 is baseline risk

RR = 1.3 means 30% higher rate/risk in treatment vs control on average. More interpretable than OR. RR = 1 means no effect. For protective effects: RR = 0.7 means 30% reduction

Alpha (exchangeable): constant within-cluster correlation. Rho (AR(1)): correlation between adjacent timepoints. Higher values indicate more clustering/autocorrelation, justifying GEE over independence

RD = p1 - p0 (absolute difference in proportions). Most interpretable for public health: RD = 0.05 means 5% absolute increase. Number needed to treat: NNT = 1/RD

Recommended Metric: For binomial: report both OR and predicted probabilities (or RD) for interpretability. For counts: RR with 95% CI. For continuous: mean difference. Always report correlation parameter (alpha/rho) to quantify clustering and justify GEE approach

Small
0.2

Medium
0.5

Large
0.8

0.50

For binomial: report both OR and predicted probabilities (or RD) for interpretability. For counts: RR with 95% CI. For continuous: mean difference. Always report correlation parameter (alpha/rho) to quantify clustering and justify GEE approach

Recommended Measure

Available Metrics

ReportUse For binomial: report both OR and predicted probabilities (or RD) for interpretability. For counts: RR with 95% CI. For continuous: mean difference. Always report correlation parameter (alpha/rho) to quantify clustering and justify GEE approach to represent clinical impact magnitude.

08Statistical Power

Sample Size

Guidelines for minimum sample requirements and power analysis parameters.

An underpowered study is an ethical failure. Respect the data by collecting enough of it.

Power Protocol

Floor Requirements

The '30-Cluster Shield': A minimum of 30 independent clusters is required to ensure that the 'Robust' standard errors (Sandwich Estimators) stabilize and provide valid p-values.

Effect Size	Parameters	Required n
Small Effect	Odds Ratio = 1.5 (Small)	n ≈ 45 clusters
Medium Effect	Odds Ratio = 2.5 (Medium)	n ≈ 25 clusters
Large Effect	Odds Ratio = 4.0 (Large)	n ≈ 12 clusters

Key considerations

The 'ICC Penalty': Higher internal correlation (ρ) within clusters *reduces* your effective sample size. If subjects within a site are highly similar, you must recruit 20-40% more clusters to achieve the same discovery authority.

G*Power StrategyBenchmark: Z-tests → GEE (Population-Averaged). Parameters: 2 groups, 3 timepoints, ICC = .30, α = .05, Power = .80. Note: GEE power depends on the 'Design Effect'—calculated as 1 + (m-1)ρ.

09APA narrative blueprint

Reporting

How to compile statistical results into publication prose matching APA and journal style guides.

Data does not speak for itself. It requires a translator. Be clear, be precise, be honest.

Narrative Arc

Worked APA paragraph example

Generalized estimating equations (GEE) with binomial family (logit link), exchangeable correlation structure, and robust standard errors were used to analyze 6-month smoking abstinence in a cluster-randomized trial (n=1000 patients across 40 clinics). The model adjusted for age and baseline smoking. Proactive telephone counseling significantly increased abstinence compared to standard care (OR = 1.40, 95% CI [1.03, 1.91], Wald χ² = 4.5, p = .034). At average age (45 years) and baseline smoking (18 cigarettes/day), predicted abstinence was 16.3% for counseling versus 12.2% for standard care, representing a 4.1 percentage point increase (NNT = 24). The estimated within-clinic correlation was 0.08, indicating 8% of variance was due to clinic-level factors. Residual diagnostics showed adequate model fit. These findings support telephone counseling as an effective population-level intervention for smoking cessation in primary care.

Reusable template

Generalized estimating equations (GEE) with family, e.g., 'Gaussian/binomial/Poisson' family (link function), correlation structure, e.g., 'exchangeable/AR(1)/unstructured' correlation structure, and robust (sandwich) standard errors were used to analyze outcome over time/clustered by X. Sample description: n participants, k clusters, repeated measures. The model included list fixed effects predictors. If applicable: Correlation structure selection: 'AR(1) was selected over exchangeable based on lower QIC (QIC_AR1 = X vs QIC_exch = Y).' OR 'Sensitivity analyses using multiple correlation structures produced similar results (estimates within X units).'. Results showed describe main finding: significant/non-significant effect, Wald χ² = X.X, df = X, p = .XXX. For binomial: The population-average odds ratio was OR = X.XX (95% CI X.XX, X.XX), indicating interpretation. Predicted probabilities for interpretability: 'Predicted abstinence was X% for treatment A vs Y% for treatment B, representing a Z percentage point difference.'. For counts: The rate ratio was RR = X.XX (95% CI X.XX, X.XX), indicating interpretation. For continuous: The mean difference was X.XX units (95% CI X.XX, X.XX). The estimated correlation parameter, e.g., 'within-cluster correlation (alpha)' or 'AR(1) parameter (rho)' was X.XX, indicating interpretation of clustering/autocorrelation. Model diagnostics: 'Residual plots showed no systematic patterns. Comparison of model-based and robust SEs suggested [adequate/inadequate correlation structure specification.']. Conclude with substantive interpretation in research context.

Essential statistics to report

Model specification (family, link, correlation structure)
Sample size (n observations, k clusters, cluster sizes)
Fixed effects estimates (β), robust SEs, Wald statistics, p-values
Exponentiated coefficients (OR or RR) with 95% CIs for non-linear links
Correlation parameter estimate (alpha, rho) with interpretation
Model fit: QIC for correlation structure comparison
Diagnostics: convergence, residual plots, comparison of model-based vs robust SEs
For binomial outcomes with common events: predicted probabilities and/or risk differences

10Exhibit Builder

Manuscript Lab

Copy standard summary tables and forensic reporting grids to outline analysis details.

Table 1: GEE for Population-Averaged Binary Outcomes

Predictor	B (Logit)	Robust SE	Wald χ²	p	OR	95% CI
(Intercept)	0.12	0.08	2.25	.134	1.12	[0.96, 1.31]
Time	0.45	0.10	20.25	< .001	1.57	[1.30, 1.90]
Intervention	0.85	0.22	14.90	< .001	2.34	[1.52, 3.60]

Note. Outcome: Clinical Success (1/0). Correlation Structure: Exchangeable. N = 200.

Intervention OR = 2.34Identifies the Population Benefit. On average, the odds of success across the entire population increase by 134% when deploying the intervention.

Header glossary

The 'Consistency Guard'. Adjusted for the fact that observations within the same person are correlated.

The Correlation Assumption. Assumes that any two measurements from the same person are equally correlated, regardless of time interval.

11Algorithmic Logic

Command Center

Syntax libraries and function parameters for executing calculations in stats packages.

Code is the modern laboratory. Clean execution ensures reproducible discovery.

Execution Engine

# 1. Fit GEE Model
model <- geepack::geeglm(success ~ time + treatment, data = df, 
                         id = subject_id, family = 'binomial', 
                         corstr = 'exchangeable')
summary(model)

Library stack

geepackmultcomp

Python

statsmodels

Elite Forensic Strike

GEE is the 'Public Health' model. Use it when you want to make statements about the 'Average American' rather than how 'Individual A' changed relative to themselves.

# Audit QIC (Equivalent of AIC for GEE) to choose correlation structure
MESS::QIC(model)

12The Over-adjustment Trap

Common Mistakes

Analytical caveats and corrections to maintain modeling integrity.

Wisdom is learning from the failures of others. Anticipate the error before it occurs.

Defensive Logic

Why it's wrong

Sandwich (robust) SEs are asymptotically valid, requiring many clusters. With few clusters, sandwich SEs are anticonservative (too small), leading to inflated Type I error and falsely narrow confidence intervals. Simulation studies show coverage of 95% CIs drops to 85-90% with k < 20. The problem is worse with small cluster sizes

The correction

If k < 30: (1) Use small-sample SE corrections: Mancl-DeRouen (MD), Fay-Graubard (FG), or Kauermann-Carroll (KC) adjustments (available in R geesmv package, SAS PROC GENMOD with EMPIRICAL=MCALL); (2) Use GLMM instead (better small-sample properties if random effects normally distributed); (3) Use bootstrap SEs with cluster resampling; (4) Report results with caution, acknowledge limited precision. Minimum recommended: k ≥ 30-40 for GEE with sandwich SEs

Why it's wrong

GEE estimates population-average effects: average effect across entire population, averaging over between-cluster variability. GLMM estimates subject-specific effects: effect for a given subject, holding random effects constant. Due to non-linearity in logistic/Poisson models, these differ numerically. Example: GEE OR = 2.0 (population-average) might correspond to GLMM OR = 2.5 (subject-specific). Misinterpretation leads to incorrect inference and policy recommendations

The correction

Always clarify interpretation: (1) GEE results: 'The population-average odds ratio is X, meaning on average across the population...'; (2) GLMM results: 'The subject-specific odds ratio is X, meaning for a given individual...'. Choose method based on research question: GEE for population-level policy (e.g., public health interventions), GLMM for individual-level prediction (e.g., clinical decision-making). Both are valid but answer different questions. For continuous outcomes with identity link, marginal = conditional effects, so GEE and LMM give similar estimates

Why it's wrong

GEE with standard methods assumes MCAR: missingness is unrelated to observed or unobserved outcomes. If data are MAR (missing at random - missingness depends on observed covariates) or MNAR (missing not at random - missingness depends on unobserved outcomes), standard GEE produces biased estimates. Common scenario: participants with worse outcomes drop out (MNAR), leading to underestimation of treatment effects. Even with MAR, GEE can be biased if missing data mechanism not addressed

The correction

Investigate missingness: (1) Compare baseline characteristics between completers vs those with missing data; (2) Use Little's MCAR test or logistic regression predicting missingness from covariates/outcomes. If MAR: (3) Use inverse probability weighting (IPW-GEE): weight observations by inverse of probability of being observed, estimated from logistic regression; (4) Use multiple imputation then GEE (MI-GEE); (5) Include baseline predictors of missingness as covariates in model. If MNAR suspected: (6) Sensitivity analyses with different missingness assumptions; (7) Consider selection models or pattern-mixture models; (8) Use GLMM with maximum likelihood (handles MAR better under correct distributional assumptions). Always report proportion and pattern of missing data

Why it's wrong

While QIC (quasi-likelihood information criterion) is useful for comparing correlation structures (lower = better), (1) QIC differences of <2-3 are often negligible; (2) QIC can be unreliable with small samples; (3) Even if correlation structure is misspecified, robust SEs ensure valid inference (key advantage of GEE). Overreliance on single 'best' QIC ignores robustness of GEE. Results should be similar across reasonable correlation structures; if not, raises concerns about model specification or data issues

The correction

Use QIC as guide but not sole criterion: (1) Fit multiple plausible correlation structures (independence, exchangeable, AR(1), unstructured); (2) Compare effect estimates (βs), not just QIC: if estimates similar (e.g., within 10-20%), results are robust; if estimates differ substantially, investigate why (outliers? mean structure misspecified?); (3) Choose structure based on: (a) Study design (longitudinal → AR(1), cross-sectional clusters → exchangeable), (b) QIC, (c) Interpretability; (4) Report results from primary structure + sensitivity analysis: 'Using AR(1) (QIC=450) vs exchangeable (QIC=455), estimates differed by <5%, supporting robustness'. Remember: even with misspecified correlation, sandwich SEs remain valid (though may be less efficient)

Why it's wrong

Uncentered predictors make intercept uninterpretable (outcome when all predictors = 0, often meaningless, e.g., age = 0). With interactions, centering prevents multicollinearity between main effects and interaction terms, improving estimation stability. Centering also aids convergence in iterative algorithms. Without centering, correlation structure estimation can be unstable due to scale differences

The correction

Grand-mean center continuous predictors before fitting GEE: predictor_c = predictor - mean(predictor). Benefits: (1) Intercept = expected outcome at mean predictor values (interpretable); (2) Reduces multicollinearity in interaction models; (3) Improves numerical stability and convergence; (4) Makes effect estimates more interpretable (effect of 1-unit change from mean, not from zero). For interactions: center all components before creating interaction term. Standardizing (z-scores) additionally allows comparing effect sizes across predictors (1-unit change = 1 SD)

Why it's wrong

GEE and GLMM are complementary, not competing methods. They answer different questions and have different strengths/weaknesses. GEE: estimates population-average effects, robust to correlation misspecification (sandwich SEs), requires fewer assumptions, better with many small clusters, but assumes MCAR and needs large k. GLMM: estimates subject-specific effects, handles MAR missingness (via ML), allows prediction for new clusters, models data-generating process, but requires random effects normality and correlation structure correctly specified. Neither is universally 'better'

The correction

Choose method based on: (1) Research question: Population-level policy (public health, epidemiology) → GEE; Individual-level prediction (clinical decision-making, personalized medicine) → GLMM; (2) Sample size: k < 30 → GLMM may be better; k ≥ 40 → GEE valid; (3) Missing data: MAR → GLMM better (ML handles MAR); MCAR → either works; MNAR → both problematic, use sensitivity analyses; (4) Cluster structure: Many small clusters → GEE; Fewer larger clusters → GLMM; Prediction for new clusters → GLMM. Report method choice rationale in methods section. Consider fitting both and comparing results to assess robustness (expect similar results if assumptions met)

Why it's wrong

The correlation parameter quantifies within-cluster dependence and justifies use of GEE over simpler methods (GLM with independence). If correlation ≈ 0, GEE is unnecessary (reduces power vs GLM). Reviewers expect correlation parameter reported. Correlation magnitude informs design of future studies (sample size calculations require ICC/correlation). Failing to report suggests analyst didn't verify clustering was present

The correction

Always extract and report estimated correlation: (1) For exchangeable: 'The estimated within-cluster correlation (alpha) was X.XX, indicating X% of variance was due to clustering'; (2) For AR(1): 'The estimated AR(1) parameter (rho) was X.XX, indicating moderate/strong autocorrelation between adjacent timepoints'; (3) For unstructured: report matrix or range. Interpret: alpha/rho = 0.05 (weak clustering), 0.10-0.20 (moderate), >0.20 (strong). If correlation near zero: acknowledge GEE may be conservative; consider simpler GLM. High correlation (>0.5): strong clustering, definitely need GEE/GLMM. Correlation informs: (a) Whether clustering matters, (b) Design effect for power analyses, (c) Validity of independence assumption

Why it's wrong

While robust SEs ensure valid inference even with misspecified correlation (Type I error control), using independence when strong correlation exists is very inefficient: much wider CIs, reduced power. Can miss real effects. Independence assumes no within-cluster correlation, ignoring data structure entirely. Although technically valid with sandwich SEs, it's a poor choice when better structures available (exchangeable, AR(1)). With strong clustering (ICC > 0.2), independence can lose 30-50% power vs correctly specified structure

The correction

Always try multiple correlation structures: (1) Start with structure matching design: longitudinal → AR(1) or exchangeable; cross-sectional clusters → exchangeable; (2) Compare QIC across structures; (3) Examine estimated correlation: if α/ρ > 0.10, clustering is substantial and appropriate structure improves efficiency; (4) Report results from best structure + sensitivity under independence to show robustness. Use independence only when: (a) Correlation truly near zero (α < 0.05), or (b) As sensitivity analysis to confirm robustness, not as primary analysis. Remember: goal is efficiency (power) not just validity. Correct correlation improves power to detect real effects while maintaining Type I error control via robust SEs

13Academic Lineage

References

Scholarly lineage and citation keys grounding the statistical framework.

We stand on the shoulders of giants. Honor the source of the method.

Academic Lineage

[1]

Liang, K. Y., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13-22.

Seminal paper introducing GEE methodology. Essential theoretical foundation for population-average models with correlated data

doi: 10.1093/biomet/73.1.13

[2]

Zeger, S. L., Liang, K. Y., & Albert, P. S. (1988). Models for longitudinal data: a generalized estimating equation approach. Biometrics, 44(4), 1049-1060.

Extension of GEE to longitudinal data with time-dependent covariates. Discusses correlation structures and sandwich estimators

doi: 10.2307/2531734

[3]

Hanley, J. A., Negassa, A., Edwardes, M. D. D., & Forrester, J. E. (2003). Statistical analysis of correlated data using generalized estimating equations: An orientation. American Journal of Epidemiology, 157(4), 364-375.

Accessible tutorial on GEE for epidemiologists. Covers practical issues: correlation structures, interpretation, comparison with GLMM

doi: 10.1093/aje/kwf215

[4]

Hubbard, A. E., Ahern, J., Fleischer, N. L., Van der Laan, M., Lippman, S. A., Jewell, N., ... & Satariano, W. A. (2010). To GEE or not to GEE: comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology, 21(4), 467-474.

Comparison of GEE vs GLMM: when to use each, marginal vs conditional interpretation, simulation study results

doi: 10.1097/EDE.0b013e3181caeb90

[5]

McNeish, D. (2019). Small sample methods for multilevel modeling: A colloquial elucidation of REML and the Kenward-Roger correction. Multivariate Behavioral Research, 54(5), 661-670.

Discussion of small-sample issues in clustered data, including why GEE problematic with k < 30 and small-sample corrections

doi: 10.1080/00273171.2019.1585672

[6]

Stead, L. F., Hartmann-Boyce, J., Perera, R., & Lancaster, T. (2013). Telephone counselling for smoking cessation. Cochrane Database of Systematic Reviews, 8, CD002850.

Meta-analysis showing telephone counseling increases cessation (OR ≈ 1.4). Basis for Example 2

doi: 10.1002/14651858.CD002850.pub3

GEE sees the forest, not the trees. Use it when the global health of the population is more important than the trajectory of any single individual.

The Interpretive Rigor Directive

statminds · GEEMind reference · v2.2 · updated 2026-01-1715 of 15 sections

GEE (Generalized Estimating Equations)

What is it?

Goals & Indications

Core Idea Diagram

Hypotheses

How it works

Assumptions

Important Note

Worked Example

GEE Working Correlation Matrix

Hypotheses

Assumptions

Diagnostics

Applied Minds

Blood Pressure Reduction in Hypertension (Continuous Longitudinal Outcome)

Alternatives

Post-hoc

Effect Size

Sample Size

Reporting

Manuscript Lab

Command Center

Common Mistakes

References

Related Tests

simpler

more complex

conditional approach

for subject specific effects

with random slopes

continuous outcome

non normal distributions

small samples