Cluster-Robust Standard Errors

In plain language

When you analyze data where multiple rows come from the same source — say, five hospital visits from the same patient, or ten patients from the same clinic — those rows are not truly independent of one another. Standard statistical software assumes every row is independent and, as a result, calculates standard errors that are too small, producing confidence intervals that look more precise than the data actually support. Cluster-robust standard errors fix this by grouping all rows from the same source together and recalculating the uncertainty around your estimate to reflect how much truly independent information you have. The result is wider, more honest confidence intervals — the point estimate (for example, a risk difference or odds ratio) does not change, only the uncertainty around it.

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

Any model fit on data where independence holds at the cluster level but not the row level: (1) repeated measures / longitudinal outcomes with multiple records per patient; (2) recurrent-event survival (Andersen-Gill, WLW, PWP), where each patient supplies several at-risk intervals and events; (3) propensity-score matched cohorts, where the matched set is the cluster and within-set outcomes are correlated by construction; (4) multilevel real-world data — patients within clinicians, hospitals, plans, or geographic regions, or siblings within families; (5) pooled / stacked person-time designs (pooled logistic for hazards, case-time-control) where each patient contributes many person-period rows.

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.

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

• Run the regression treating all 12 rows as independent. The model estimates a risk difference of -0.33 (Hospital A has 33 percentage points lower readmission than Hospital B).

• The naive standard error for that estimate comes out to 0.05, giving a 95% confidence interval of roughly -0.43 to -0.23 — it looks very precise.

• But patients within the same hospital are not independent: they share discharge nurses, post-discharge call protocols, and local referral networks. The 6 Hospital A rows are really only 1 independent unit of information about Hospital A, and the 6 Hospital B rows are 1 independent unit about Hospital B.

• Apply cluster-robust standard errors, grouping by hospital_id. The software now sums the score contributions within each hospital before computing the variance, reflecting that we have 2 independent clusters, not 12 independent rows.

• The cluster-robust standard error is 0.09 — nearly twice as large as the naive version — and the 95% confidence interval widens to roughly -0.51 to -0.15.

• The point estimate (-0.33) is identical. Only the uncertainty changed: the naive interval was falsely narrow because it double-counted correlated observations as if they were independent.

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.

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