Cluster-Robust Standard Errors

Cluster-robust standard errors (CRSE)

replace the model-based variance of a regression coefficient with a sandwich estimator that is consistent under arbitrary within-cluster correlation. The "bread" is the usual inverse-information (or inverse-Hessian) matrix; the "meat" is built from the sum of cluster-level score contributions rather than individual-observation contributions. The point estimate (hazard ratio, odds ratio, rate ratio, risk difference) does not change — only its standard error, and therefore its Wald confidence interval and p-value.

CRSE is the variance counterpart to the design fact that, in real-world data, rows are not independent: the same `person_id` contributes multiple visits or recurrent events, propensity-score matching creates correlated matched sets, and patients nest within clinicians, facilities, and health plans.

Core conceptual distinction — what CRSE does and does not estimate

CRSE targets the marginal (population-average) variance of a coefficient under a working model. It is agnostic about the source and form of the correlation: you do not specify a within-cluster correlation structure correctly to get valid inference, you only need clusters to be independent of one another and reasonably numerous.

This is the key difference from a random-effects / mixed model, which models the correlation explicitly (a random intercept per cluster) and, in nonlinear models (logistic, Poisson), changes the estimand itself from a marginal to a conditional (subject-specific) effect. A GEE with an independence working correlation plus the robust sandwich is numerically identical to fitting the naive model and applying CRSE; richer working correlations (exchangeable, AR-1) only buy efficiency, not validity.

So the decision is not "CRSE vs GEE" — CRSE is the variance engine of GEE — but "marginal effect with robust variance" (CRSE/GEE) vs "conditional effect with a random effect" (mixed model). Pick the estimand first; the variance method follows.

Pros, cons, and trade-offs

• vs naive (model-based / independence) standard errors: CRSE corrects the near-universal under-estimation of SEs that occurs when correlated rows are treated as independent, which otherwise yields anticonservative CIs and inflated type-I error (a recurrent-event or repeated-visit analysis run without clustering will routinely report CIs ~30-60% too narrow). Cost: essentially none at the point-estimate level and trivial computationally; the only real risk is using CRSE when you have too few clusters (see below). Always prefer CRSE over naive SEs whenever the unit of analysis is finer than the unit of independence.
• vs random-effects / mixed-effects models: CRSE makes no distributional assumption about the cluster effects and keeps a clean marginal interpretation, which is usually what a regulator or payer wants ("the average effect in the population"). Cost: it discards the efficiency gains a correct random-effects model provides, gives no cluster-level prediction (no BLUPs), and is downward-biased with few clusters, whereas a well-specified mixed model can be more efficient and handles small numbers of large clusters better. Prefer CRSE/GEE when the target is a marginal effect and the number of clusters is large; prefer a mixed model when you want subject-specific effects, cluster-level prediction, or have few clusters with many members.
• vs heteroskedasticity-robust (HC0/HC1, "White") SEs: ordinary robust SEs correct for heteroskedasticity but still assume independent rows; they are the special case of CRSE where every cluster has size one. Using HC instead of CR on clustered data leaves the SEs as wrong as the naive estimator. Prefer CRSE whenever any cluster has more than one observation.
• vs cluster bootstrap: the bootstrap (resampling whole clusters with replacement) is a robust alternative that often performs better with moderate cluster counts and for non-smooth statistics. Cost: computation and seed management. Prefer the cluster bootstrap or a small-sample correction (CR2/CR3) over the standard CRSE when clusters are few.

When NOT to use — and when CRSE is actively misleading

• Too few clusters. The standard CRSE is consistent only as the number of clusters grows; with fewer than roughly 40-50 clusters it is downward-biased, producing CIs that are too narrow and type-I error well above nominal. A multi-site study with 8 hospitals clustered at the hospital level is the classic trap. Use a bias-reduced estimator (CR2/CR3, the Mancl-DeRouen or Kauermann-Carroll corrections) with Satterthwaite degrees of freedom, or a cluster wild bootstrap. Clustering at a finer level with many units (e.g., patient rather than site) does not have this problem.
• You actually want a subject-specific effect. If the question is conditional ("for a given patient/clinician, what is the effect"), CRSE on a marginal model answers a different question; fit a mixed model.
• The clustering variable is endogenous or post-treatment. Clustering on something affected by exposure (e.g., the treating specialist chosen because of the drug) does not fix bias and can distort variance; clusters must be exogenous, independent groupings.
• Singleton-heavy or one-observation-per-cluster data. If essentially every cluster has size one, CRSE collapses to HC robust SEs and adds nothing — and small numbers of large clusters mixed with many singletons can behave erratically.
• CRSE does not cure confounding or model misspecification of the mean. It fixes the variance, not the point estimate; a biased coefficient with a correct CRSE is still biased, now with a confidently wrong interval.

Data-source operational depth

• Claims (FFS or commercial): The natural cluster is `person_id`, because a single enrollee generates many pharmacy and medical claims and, in recurrent-event or pooled-logistic designs, many analysis rows. Cluster on `person_id`, not on the claim. Higher-level clustering (plan, PBM, provider TIN, geographic region) matters for utilization and cost outcomes that are correlated by formulary, network, and local practice. Failure modes: plan-switching moves a person across plan clusters mid-follow-up — decide a priori whether the cluster is the person (usual choice) or the plan-spell; Medicare Advantage-only person-time lacks fee-for-service claims, so apparent "independent" intervals are really missingness, and clustering cannot repair an undercounted denominator; few large plans (a database dominated by 5-10 payers) recreates the few-clusters problem if you cluster at the plan level — cluster at the person level or use a small-cluster correction.
• EHR: Encounter-driven capture makes the patient the primary cluster, but clinic/site is a strong second-level cluster because workflow, coding, and care patterns are shared. With a handful of contributing sites, site-level CRSE is unreliable (few clusters) — prefer patient-level clustering, a mixed model with a site random effect, or a wild bootstrap. External-care leakage means some within-patient correlation is unobserved; CRSE only accounts for correlation among the rows you actually have.
• Registry: Patients nest within enrolling centers; center is a meaningful cluster for adjudicated outcomes and for quality-of-care contrasts. Registries often have moderate center counts (10-30), which sits squarely in the few-clusters danger zone — report a CR2/Satterthwaite or bootstrap variant and state the cluster count.
• Linked claims-EHR-vital-records: Multiple plausible cluster levels coexist (person, site, plan). Pre-specify the clustering level in the SAP; nesting (patient within site within region) usually calls for clustering at the highest level at which independence is credible, or for two-way / multiway clustering when neither dimension nests cleanly (e.g., patients seen across several facilities).

Worked claims example — two scenarios where the choice changes inference

Scenario A (within-person recurrent events). Question: rate of recurrent COPD exacerbations on inhaled corticosteroid + LABA vs LABA alone in a commercial + Medicare FFS cohort. Each `person_id` is followed from the index fill and contributes multiple exacerbation events and at-risk intervals, so rows within a person are correlated.

Fit an Andersen-Gill Cox model on `(start, stop, event)` person-interval rows with the treatment covariate; the naive SE treats every interval as independent and understates the true SE. Cluster on `person_id` (`COVS(AGGREGATE)` / `cluster(id)`): the HR is unchanged, the SE widens appropriately, and the 95% CI now reflects that 10,000 rows came from 3,200 people. Sensitivity: compare to a shared-frailty (random-effect) model to confirm the marginal and conditional stories agree. Scenario B (1:k PS-matched cohort).

Question: 1-year all-cause mortality, drug A vs drug B, after 1:1 propensity- score matching among adults with `>=`365 days of continuous A/B/D enrollment and a drug-free washout. Matching induces within-pair correlation (matched patients share covariate values by construction). Fit the outcome model on the matched cohort and cluster on the matched-set id, not the person, so the variance accounts for the matched design (Austin 2014).

Failure mode to avoid: ignoring the matched-set clustering yields anticonservative SEs; over- clustering (e.g., on plan when the design unit is the pair) targets the wrong correlation. Report the cluster count (number of matched sets) and, if few large clusters appear at any candidate level, switch to CR2/Satterthwaite or a cluster wild bootstrap before trusting the interval.

Interpreting the output

Consider a two-hospital study: 6 patients each in Hospital A (all treated) and Hospital B (all comparators). The observed risk difference is −0.33. Under a naive independence assumption the SE is 0.05, yielding a 95% CI of approximately −0.43 to −0.23. Clustering on hospital raises the SE to 0.09, yielding a 95% CI of approximately −0.51 to −0.15.

(1) Formal statistical interpretation. The point estimate of −0.33 is identical under both SE approaches — clustering corrects the uncertainty, not the central estimate. The naive SE of 0.05 underestimates true variability because it treats the 12 patients as 12 independent draws; within-hospital patients share unmeasured environmental and practice factors, so the effective information is closer to two hospitals than twelve patients. The cluster-robust sandwich estimator targets the between-cluster variance, producing a wider CI calibrated to the actual clustering structure.

With only two clusters, however, the sandwich estimator itself is unreliable; CR2/Satterthwaite or a cluster wild bootstrap is preferred.

(2) Practical interpretation for a decision-maker. The risk difference is the same (−0.33), but the honest uncertainty range is wider: −0.51 to −0.15 rather than −0.43 to −0.23. In a two-hospital study, a result that looks statistically significant under naive SEs may not survive correct variance estimation. Report the cluster count alongside the cluster-robust interval so readers can judge whether the precision claim is credible.

Decision diagram

flowchart TD
  Data[Analysis table<br/>rows finer than the unit of independence] --> Q{Where does independence hold?}
  Q -->|Same patient, many rows<br/>visits / recurrent events / pooled person-time| Pat[Cluster on person_id]
  Q -->|PS-matched cohort| Match[Cluster on matched-set id]
  Q -->|Patients within<br/>clinician / site / plan / region| High[Cluster on the higher level]
  Pat --> N{Number of clusters}
  Match --> N
  High --> N
  N -->|Many >= 40-50| Std[Standard cluster sandwich<br/>point estimate unchanged, SE widened]
  N -->|Few < 40-50| Small[CR2/CR3 + Satterthwaite df<br/>or cluster wild bootstrap]
  Std --> Report[Report robust CI / p-value<br/>+ cluster count + cluster level]
  Small --> Report

flowchart LR
  subgraph Naive[Naive variance]
    I1[Row 1] --> M1[Sum over independent rows]
    I2[Row 2] --> M1
    I3[Row 3] --> M1
  end
  subgraph Cluster[Cluster-robust variance]
    C1[Patient A rows<br/>score sum] --> M2[Sum over independent CLUSTERS]
    C2[Patient B rows<br/>score sum] --> M2
  end
  M1 --> B1[Bread x Meat x Bread<br/>too small -> CI too narrow]
  M2 --> B2[Bread x Meat x Bread<br/>valid -> CI correct]

Worked example

Scenario

A researcher is studying whether patients admitted to Hospital A have a lower 30-day readmission rate than patients admitted to Hospital B. The dataset has 12 patients — 6 from each hospital. Because patients at the same hospital share the same doctors, protocols, and discharge practices, their outcomes are correlated. Fitting a regression that ignores this clustering produces a standard error that is too small and a confidence interval that is falsely narrow.

Dataset

One row per patient. The cluster variable is hospital_id, which groups patients who share a care environment.

Steps

Result

Naive SE: 0.05, 95% CI approximately -0.43 to -0.23. Cluster-robust SE: 0.09, 95% CI approximately -0.51 to -0.15. The lesson: ignoring clustering made the result look twice as precise as it really was. With only 2 hospitals, even the cluster-robust SE has limitations (very few clusters), but it is far closer to the truth than the naive version.

Trade-offs

Runnable example

import pandas as pd
import statsmodels.formula.api as smf
from lifelines import CoxPHFitter

# --- Repeated-measures / pooled logistic: cluster on person_id ---
# Naive SEs assume independent person-period rows and are too narrow.
model = smf.logit("y ~ treat + age + sex + comorbidity_score", data=df)
naive = model.fit(disp=0)                                  # model-based (anticonservative) SEs
robust = model.fit(disp=0, cov_type="cluster",            # cluster sandwich
                   cov_kwds={"groups": df["person_id"]})
# Coefficients are identical; compare the standard errors:
print(naive.bse["treat"], robust.bse["treat"])
print(robust.summary())                                    # report the robust CI/p-value

# --- Recurrent-event (Andersen-Gill) Cox: cluster on person_id ---
# surv_df is in counting-process (start, stop, event) form, multiple rows per person.
cph = CoxPHFitter()
cph.fit(surv_df, duration_col="stop", entry_col="start",
        event_col="event", cluster_col="person_id",        # robust SE for within-person correlation
        formula="treat + age + sex + comorbidity_score")
cph.print_summary()                                         # 'robust' SE column reflects clustering

# --- PS-matched cohort: cluster on the matched-set id, not person ---
matched_fit = smf.logit("y ~ treat", data=matched).fit(
    disp=0, cov_type="cluster", cov_kwds={"groups": matched["matchid"]})

library(sandwich); library(lmtest); library(geepack); library(survival)

# --- GLM + cluster-robust (CR) variance on person_id ---
fit <- glm(y ~ treat + age + sex + comorbidity_score, family = binomial, data = df)
coeftest(fit, vcov = vcovCL(fit, cluster = ~ person_id))            # CR sandwich
# Few clusters: bias-reduced CR2 variant
coeftest(fit, vcov = vcovCL(fit, cluster = ~ person_id, type = "HC2"))

# --- Equivalent population-average GEE (independence working corr -> same robust SE) ---
df <- df[order(df$person_id), ]
gee <- geeglm(y ~ treat + age + sex + comorbidity_score, family = binomial,
              id = person_id, corstr = "independence", data = df)
summary(gee)                                                        # 'Std.err' is the robust (sandwich) SE

# --- Recurrent-event Cox (Andersen-Gill) with robust clustered variance ---
ag <- coxph(Surv(start, stop, event) ~ treat + age + sex + comorbidity_score +
              cluster(person_id), data = surv_df)                   # robust SE for within-person events
summary(ag)                                                         # 'robust se' column

# --- PS-matched cohort: cluster on matchid ---
mfit <- glm(y ~ treat, family = binomial, data = matched)
coeftest(mfit, vcov = vcovCL(mfit, cluster = ~ matchid))

/* --- GEE/GENMOD: empirical (cluster-robust) SE under an independence working correlation --- */
proc sort data=work.df; by person_id; run;
proc genmod data=work.df descending;
  class person_id;
  model y = treat age sex comorbidity_score / dist=binomial link=logit;
  repeated subject=person_id / type=ind;     /* type=ind -> robust SE = cluster sandwich, point est unchanged */
run;                                          /* read the 'Empirical Standard Error Estimates' table */

/* --- Recurrent-event Cox (Andersen-Gill) with robust clustered variance --- */
proc phreg data=work.surv covsandwich(aggregate);
  model (start, stop)*event(0) = treat age sex comorbidity_score;
  id person_id;                              /* aggregate score residuals within person -> robust SE */
run;                                          /* compare 'StdErrRatio' = robust/model-based SE */

/* --- Clustered binary outcome via survey procedures (alternative) --- */
proc surveylogistic data=work.df;
  cluster person_id;                          /* design-based cluster-robust variance */
  model y(event='1') = treat age sex comorbidity_score;
run;

/* --- PS-matched cohort: cluster on the matched-set id --- */
proc surveylogistic data=work.matched;
  cluster matchid;                            /* matched-set is the cluster (Austin 2014) */
  model y(event='1') = treat;
run;

Citations

1. [1]Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13-22.
2. [2]White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838.
3. [3]Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126-134.

1. [4]Petersen MA. Estimating standard errors in finance panel data sets: comparing approaches. Review of Financial Studies. 2009;22(1):435-480.

1. [5]Austin PC. The use of propensity score methods with survival or time-to-event outcomes: reporting measures of effect similar to those used in randomized experiments. Statistics in Medicine. 2014;33(7):1242-1258.