Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
65,873 characters · 33 sections · 55 citation commands
Design-Based Variance Estimation for Modern Heterogeneity-Robust Difference-in-Differences Estimators
Keywords: difference-in-differences, survey sampling, Taylor series linearization, influence functions, heterogeneity-robust estimation, complex survey design
Policy evaluations frequently rely on nationally representative surveys --- NHANES for health outcomes, ACS for demographics and housing, CPS for labor force dynamics, and MEPS for medical expenditure. These surveys employ stratified multi-stage cluster sampling to achieve national coverage at manageable cost. The resulting data carry two features that invalidate iid-based standard errors: observations within the same primary sampling unit (PSU) are correlated, and stratification constrains sampling variability. The ratio of design-based variance to the iid HC1 variance --- the design effect (DEFF) --- commonly ranges from 2 to 5 in health and social surveys.
This matters especially for difference-in-differences (DiD) estimation. Treatment is often assigned at geographic levels --- state policies, county programs, school-district mandates --- that encompass multiple survey PSUs, so within-PSU homogeneity in treatment status intensifies the design effect. DiD estimands involve contrasts across groups and time periods, amplifying distortions in variance estimation. Incorrect standard errors can flip significance conclusions for policy-relevant effect sizes, undermining the credibility of program evaluations.
The modern heterogeneity-robust DiD literature derives estimators and their asymptotic properties under sampling assumptions that are incompatible with complex survey designs. The foundational papers either assume iid sampling explicitly or adopt frameworks that do not incorporate strata, PSU clustering, or finite population corrections:
The most comprehensive recent synthesis --- roth2023 --- discusses design-based inference in the treatment-assignment sense atheyimbens2022, where randomness comes from which units receive treatment. It does not address survey sampling design, survey weights, or strata/PSU/FPC-based variance estimation.
A note on terminology is warranted. The recent DiD literature uses “design-based” to refer to treatment-assignment design atheyimbens2022, where uncertainty arises from which units receive treatment. Throughout this paper, “design-based” refers to survey sampling design binder1983, where uncertainty arises from which units are sampled. Same term, different referent.
Existing software implementations reflect this theoretical gap. R's did package callawaysantanna2021 accepts a weightsname parameter for point estimation and supports cluster-level multiplier bootstrap, but does not account for stratification or finite population corrections. Stata's csdid accepts pweight for point estimation but does not support the svy: prefix --- there is no mechanism for strata or FPC. Other packages --- did_multiplegt_dyn, \texttt{eventstudyinteract}, \texttt{didimputation} --- likewise support weights for point estimation and allow cluster-robust inference, but none provides full survey-design variance estimation that jointly accounts for strata, PSU clustering, and finite population corrections.
The survey statistics literature has developed design-based variance theory for smooth functionals binder1983, demnati2004, lumley2004, and recent work has extended this to causal inference --- but primarily for cross-sectional estimators or simple two-period designs, not for modern staggered DiD.
dugoff2014 provide practical guidance on combining propensity score methods with complex surveys, but address cross-sectional treatment effects, not DiD. zeng2025 derive sandwich variance for survey-weighted propensity score estimators using influence functions --- the closest work to the bridge we describe --- but for cross-sectional IPW/augmented weighting, not staggered DiD. ye2025 study DiD with repeated cross-sectional survey data, combining propensity scores with survey weights. However, their estimator is limited to two periods and two groups, uses bootstrap-only variance (no analytical design-based derivation), and does not address the modern heterogeneity-robust estimators considered here.
To our knowledge, prior work has not explicitly applied survey-sampling Taylor linearization --- in the strata/PSU/FPC sense --- to the influence-function representations of modern heterogeneity-robust DiD estimators.
This paper bridges the two literatures. The core argument (Section (ref)) is a verification step: we show that the influence-function representations established in the original papers for each modern DiD estimator satisfy Binder's smoothness conditions, so that the standard survey-statistics Taylor linearization applies. Binder's theorem then implies that the stratified-cluster variance formula applied to these influence-function values produces a design-consistent variance estimator. The contribution is not the proposition itself --- which is a direct corollary of binder1983 --- but the explicit identification of these estimators as covered by it, and the empirical demonstration that the resulting standard errors are well-calibrated in finite samples.
We demonstrate the practical consequences through a Monte Carlo study (Section (ref)) showing that HC1 confidence interval coverage can fall as low as 34% at the nominal 95% level, with design effects of 2--17 times the HC1 variance in the baseline scenario and exceeding 100 times under informative sampling. The simulation also decomposes where the correction has to come from. Two ingredients matter for valid coverage: the survey-weighted point estimate (which corrects estimand bias under informative sampling) and PSU-level clustering (which corrects within-cluster correlation). Combining the two gives near-nominal coverage in every scenario we examine. Adding strata and finite-population corrections is real but incremental --- a precision gain on top of an already-valid specification, not a coverage requirement. An empirical illustration using NHANES data on the ACA's dependent coverage provision (Section (ref)) shows that ignoring the survey design changes both the point estimate (by 48%) and the significance conclusion. Section (ref) describes the diff-diff Python package (v3.3.2; diffdiff2026), which implements this connection and is, to our knowledge, the first open-source implementation that jointly supports strata, PSU clustering, finite population corrections, and replicate-weight methods for these modern DiD estimators. Section (ref) discusses limitations and directions for future work.
Consider a finite population $\mathcal{U} = \{1, \ldots, N\}$. The population is partitioned into $H$ non-overlapping strata. Within stratum $h$, there are $N_h$ PSUs in the population, of which $n_h$ are sampled. This describes the standard stratified cluster design used by most federal statistical agencies. The variance estimator in Proposition (ref) is a first-stage formula that treats within-PSU totals as the unit of analysis; when observations are sub-sampled within PSUs (as in multi-stage designs), this is the standard “ultimate cluster” approximation lumley2004, which is exact when PSUs are fully enumerated. When observations are sub-sampled within PSUs, including the first-stage FPC $(1 - f_h)$ without a second-stage variance component can be slightly anti-conservative; omitting the FPC (i.e., assuming with-replacement sampling at the first stage) restores the conservative property.
Each sampled observation $i$ carries a raw sampling weight $w_i^* = 1/\pi_i$, where $\pi_i$ is the inclusion probability. The raw weights estimate the finite-population size: $\hat{N} = \sum_i w_i^* = n \bar{w}^*$, where $\bar{w}^* = n^{-1} \sum_i w_i^*$. For computation, we work with normalized weights $w_i = w_i^* / \bar{w}^*$, so that $\sum_i w_i = n$. Throughout, $w_i$ denotes the normalized weight, and we write $\hat{W} \equiv \sum_i w_i = n$ for the normalized-weight total (distinct from the population estimate $\hat{N}$). Normalization preserves relative representativeness and avoids scale dependence in regression coefficients; the ratio $w_i / \hat{W} = w_i^* / \hat{N}$ is invariant to normalization, so weighted averages take the same value either way.
The sampling fraction in stratum $h$ is $f_h = n_h / N_h$. When $f_h$ is close to 1, most of the finite population has been observed and sampling variability is reduced. The finite population correction factor $(1 - f_h)$ enters the variance formula to account for this.
We formalize the estimand as $\theta = T(F)$, where $T$ is a functional mapping a distribution to a real number (or vector) and $F$ is the finite-population distribution of potential outcomes, treatment status, and covariates. For DiD, $T$ extracts average treatment effects on the treated (ATTs) from this joint distribution.
Throughout this paper, uncertainty arises exclusively from the survey sampling mechanism: the population values $\{Y_i(0), Y_i(1), D_i, X_i\}$ are treated as fixed, and randomness enters through which units are sampled. This is the standard design-based framework of binder1983, and it is the reason that the variance estimator includes a finite population correction $(1 - f_h)$ that drives variance to zero as the sampling fraction approaches one. This framework is distinct from the treatment-assignment uncertainty considered by atheyimbens2022, and from the superpopulation framework in which parallel trends is typically stated. We adopt the convention that the parallel trends assumption holds in the finite population; the design-based variance then captures uncertainty about which members of that population we observe. Survey weights affect the estimand (they determine which population the ATT targets) as well as the variance (they determine how precisely we estimate it).
Under the survey design, the survey-weighted empirical distribution is the H\'{a}jek (self-normalized) form:
where the sum is over sampled observations and $\delta_{x_i}$ is the point mass at $x_i$. When $T$ is a smooth functional, the plug-in estimator $\hat{\theta} = T(\hat{F}_w)$ is design-consistent for $\theta = T(F)$: as the sample size grows within the finite-population asymptotic framework, $\hat{\theta}$ converges in probability to $\theta$.
The abstraction of $\theta$ as a functional of $F$ is what enables the bridge between survey statistics and DiD: both literatures reason about functionals, from different perspectives.
For regression-based estimators, the point estimates solve weighted least squares (WLS), minimizing $\sum_i w_i (Y_i - X_i'\beta)^2$. For influence-function-based estimators, point estimates are constructed from survey-weighted sample moments, replacing simple sample averages $(1/n)\sum_i$ with weighted averages $(\sum_i w_i)^{-1} \sum_i w_i$. For doubly robust and IPW variants, the same principle applies to propensity score estimation (via survey-weighted logistic regression) and outcome regression.
Under the design-based perspective adopted here, survey weights are needed to ensure that treatment effect estimates correspond to the finite population, not just the sample solon2015. Without weights, ATT estimates reflect the sample composition, which may over-represent certain strata due to the sampling design.
Table (ref) summarizes the notation used throughout the paper.
This section presents the paper's core contribution: a formal argument that design-based variance estimation is valid for modern heterogeneity-robust DiD estimators. The argument proceeds by showing that these estimators satisfy the conditions of binder1983, so that standard stratified-cluster variance formulas applied to their influence function values produce design-consistent variance estimators.
The influence function (IF) of a functional $T$ at distribution $F$ is the Gateaux derivative:
This is a property of the map $T$ and the distribution $F$. Crucially, the IF does not depend on how the sample was drawn --- the same functional $T$ has the same IF regardless of whether the data come from simple random sampling, stratified sampling, or cluster sampling. The IF characterizes each observation's first-order contribution to the estimator.
Each modern heterogeneity-robust DiD estimator can be written as $\theta = T(F)$ for a smooth functional $T$ that admits an IF representation. Table (ref) classifies the key estimators and the source of their smoothness.
We now state the main result formally. The argument combines two established bodies of theory --- influence functions for smooth functionals and design-based variance for complex surveys --- in a setting where they have not previously been connected.
Proposition (ref) provides the theoretical foundation. In practice, the variance is computed through one of two operational paths, depending on the estimator's structure.
\paragraph{Regression-based TSL sandwich.} For regression-based estimators (TWFE, Sun--Abraham, Stacked DiD), the variance-covariance matrix is the stratified-cluster sandwich binder1983:
where $W = \mathrm{diag}(w_1, \ldots, w_n)$ and the stratum-level meat is
with $T_{hj} = \sum_{i \in \text{PSU}_{hj}} w_i X_i u_i$ the PSU-level weighted score total ($u_i = Y_i - X_i'\hat{\beta}$ is the residual).
\paragraph{IF-based TSL.} For IF-based estimators (Callaway--Sant'Anna, imputation DiD, two-stage DiD), the variance is computed directly from the weighted linearized variables $\psi_i$ without the bread matrix, using equation (ref). The estimator first computes per-unit influence function values $\mathrm{IF}_i$ for each group-time cell, forms the weighted linearized variables $\psi_i = w_i \mathrm{IF}_i / \hat{W}$, aggregates across cells, and then applies the stratified-cluster formula.
\paragraph{Degrees of freedom.} Inference uses the $t$-distribution with survey degrees of freedom $\text{df} = \sum_h n_h - H$, i.e., the total number of sampled PSUs minus the number of strata. For designs with replicate weights, the degrees of freedom are computed from the rank of the replicate weight matrix, matching the convention of R's survey package lumley2004.
\paragraph{Replicate weight variance.} An alternative to TSL is replicate weight variance, which perturbs the weights and observes the resulting variation in estimates rather than linearizing the estimator. The general form is $V_{\text{rep}} = c \sum_r s_r (\hat{\theta}_r - \hat{\theta}_{\text{center}})^2$, where $\hat{\theta}_r$ is the estimate under replicate $r$ and $c$ and $s_r$ are method-specific constants. This approach is useful when pre-computed replicate weights are provided with public-use survey files (e.g., ACS PUMS, CPS, BRFSS, RECS); for surveys without packaged replicate weights --- including NHANES, where the CDC documentation recommends Taylor series linearization with the masked strata and PSU variables --- TSL is the canonical path. Five standard methods --- balanced repeated replication (BRR), Fay's BRR, jackknife (JK1, JKn), and successive difference replication (SDR) --- are detailed in Appendix (ref).
Before proceeding to empirical validation, we clarify the boundaries of the design-based framework. Survey weighting addresses the sampling design. It does not resolve every threat to valid causal inference with DiD. Practitioners should be aware of the following limitations.
\paragraph{Parallel trends.} Survey weighting ensures that treatment effect estimates target the correct population. It does not validate the parallel trends assumption. Under the design-based framework, parallel trends must hold for the population, not just the sample. If parallel trends fail in the population, survey-weighted estimates remain biased --- with correctly estimated standard errors around the wrong estimand.
\paragraph{Small-cluster asymptotics.} TSL variance requires at least 2 PSUs per stratum ($n_h \geq 2$). With few PSUs per stratum --- common in some state-based surveys --- the $t$-distribution approximation with $\text{df} = \sum_h n_h - H$ may be anti-conservative. Practitioners should assess the survey degrees of freedom directly.
\paragraph{Estimand and weights.} Survey weighting is necessary for the estimator to target the finite-population ATT rather than the sample ATT --- this is precisely the correction demonstrated in Scenario 2 of the simulation study. However, when different weighting schemes are available (e.g., person-level vs.\ household-level weights, or calibrated vs.\ design weights), the choice of weights determines which population parameter is targeted solon2015. The design-based variance estimator is consistent for the variance of whatever the weighted estimator targets, but practitioners should ensure the weighting scheme corresponds to the intended estimand.
\paragraph{Weight variability.} Highly variable weights reduce effective sample size. The Kish design effect due to unequal weighting, $\mathrm{DEFF}_w = n \sum w_i^2 / (\sum w_i)^2$, measures this: when $\mathrm{DEFF}_w \gg 1$, estimates are less precise than the nominal sample size suggests.
\paragraph{Model misspecification.} For doubly robust and IPW estimators, the IF corrections for propensity score and outcome regression uncertainty assume correct specification of at least one nuisance model. Survey weighting does not rescue a badly specified propensity score or outcome model.
We assess the practical consequences of ignoring survey design in DiD estimation through a Monte Carlo study with four scenarios, each targeting a distinct aspect of the design-based variance framework developed in Section (ref).
\paragraph{Data-generating process.} We generate staggered DiD panel data with stratified multi-stage survey structure using a purpose-built DGP. The outcome for unit $i$ in period $t$ is:
where $\alpha_i \sim N(0, \sigma^2_\alpha)$ is a unit fixed effect, $\gamma_{p(i)} \sim N(0, \sigma^2_\gamma)$ is a PSU random effect for the PSU $p(i)$ to which unit $i$ belongs, $\eta_{p(i),t} \sim N(0, \sigma^2_\eta)$ is a PSU-by-period shock, $x_{1i} \sim N(0,1)$ and $x_{2i} \sim \text{Bernoulli}(0.5)$ are covariates with coefficients $\beta_1 = 0.5$ and $\beta_2 = 0.3$, $\tau_{it}$ is the treatment effect (zero for untreated observations), and $\varepsilon_{it} \sim N(0, \sigma^2_\varepsilon)$ is idiosyncratic noise. The PSU random effects $\gamma_{p(i)}$ create intra-cluster correlation in outcome levels (relevant for cross-sectional estimators and the within-transformation); the PSU-period shocks $\eta_{p(i),t}$ are the primary driver of design effects in panel DiD, as they survive time-differencing and inflate design-based SEs relative to HC1 SEs.
The survey structure consists of $H = 5$ strata with $n_h = 8$ sampled PSUs per stratum and finite population correction $N_h = 200$. Sampling weights are stratum-based with dispersion controlled by the coefficient of variation (CV). Treatment is staggered across two cohorts (first treated at periods 3 and 5) with a never-treated comparison group.
\paragraph{Estimators.} We compare four estimators spanning the two variance paths from Section (ref):
For each estimator and replication, we compute the ATT and standard error under three inference approaches: (i) HC1, using heteroskedasticity-robust standard errors with no clustering, no weighting, and no survey design --- the default in most DiD software; (ii) PSU cluster, applying the survey-weighted point estimate with PSU-level clustering only (no strata, no FPC) --- the practitioner heuristic of cluster=psu combined with survey weights; and (iii) design-based, using TSL with strata, PSU, FPC, and survey weights as developed in Section (ref). Approaches (ii) and (iii) share the same survey-weighted point estimate; they differ only in the variance estimator.
\paragraph{Metrics.} For each cell (scenario $\times$ estimator $\times$ sample size), we report: Bias (mean estimate minus true ATT); Coverage (fraction of 95% CIs containing the true ATT); DEFF (ratio of design-based to HC1 variance, $\hat{V}_\text{design} / \hat{V}_\text{HC1}$); and RMSE (root mean squared error). Each cell consists of 2,000 independent replications at sample sizes $n \in \{500, 2{,}000, 8{,}000\}$, yielding 66,000 total replications across the study. The four estimators span both variance paths from Section (ref); simulation validation of the remaining estimators in Table (ref) is left for future work. With 2,000 replications, the Monte Carlo standard error for a coverage estimate near the nominal 95% level is approximately $\sqrt{0.95 \times 0.05 / 2000} \approx 0.5$ percentage points, so reported coverage values within roughly $\pm$1 percentage point of nominal should be interpreted as indistinguishable from $0.95$ at conventional simulation precision.
The baseline scenario generates data under standard parallel trends with moderate intra-class correlation ($\text{ICC} = 0.1$) and weight variation ($\text{CV} = 0.5$), producing a realistic survey design where design effects are present but not extreme. Table (ref) reports the results.
The core finding is stark: HC1 standard errors produce dramatically under-covered confidence intervals under complex survey designs, and the problem worsens with larger samples. At $n = 500$, HC1 coverage for CS-reg is 83.3% --- already below the nominal 95% --- but at $n = 8{,}000$ it collapses to 34.2%. This is not a small-sample artifact; it reflects the fundamental mismatch between the iid assumption embedded in HC1 SEs and the stratified-cluster structure of the data. The estimated design effect grows from 2.1$\times$ at $n = 500$ to 17.3$\times$ at $n = 8{,}000$, as larger samples make the cluster structure increasingly salient relative to the within-cluster variation.
Design-based SEs from Proposition (ref) correct this: coverage is 92.8--95.0% across all modern estimators and sample sizes, closely approximating the nominal 95% level. The correction is consistent across CS-reg, CS-DR, and Sun--Abraham, confirming that the theoretical framework applies uniformly to both IF-based and regression-based variance paths.
\paragraph{Does PSU clustering alone suffice?} A reader familiar with cluster-robust inference might ask whether simply clustering at the PSU level --- without strata or finite population corrections --- is enough once the point estimate is properly survey-weighted. The “Cluster” coverage column in Table (ref) reports this directly: the survey-weighted point estimate (which kills any informative-sampling bias) combined with PSU-level clustering only (no strata, no FPC). In Scenario 1, PSU clustering matches full design-based coverage at every sample size --- 93.8%, 95.0%, and 94.2% for CS-reg, versus 93.5%, 94.9%, and 93.7% under the full design. Strata and FPC contribute little to coverage once PSU clustering and survey weights are in place. As Scenarios 2--4 will show, the same pattern holds under informative sampling and conditional parallel trends.
TWFE shows 0% coverage at every sample size --- not because of the survey design, but because it is biased under staggered treatment adoption dechaisemartin2020. The ATT estimate of 0.31 versus the true value of 2.0 reflects the well-documented negative weighting problem. Design-based SEs correctly estimate the variance of this biased estimand but cannot rescue the point estimate.
This scenario adds two complications: (i) sampling weights are correlated with potential outcomes (informative sampling), so high-outcome units receive heavier weights; and (ii) treatment effects vary across strata (ranging from 0.6 to 3.4, with a population-weighted ATT of approximately 2.6). Under informative sampling, unweighted estimators target the sample ATT rather than the population ATT. Table (ref) reports the results.
The unweighted (no-survey-design) estimates are biased by approximately $-0.19$ across all modern estimators --- they systematically underestimate the population ATT because the sample over-represents strata with smaller treatment effects. Design-based (survey-weighted) estimates eliminate this bias, with mean bias near zero for CS-reg, CS-DR, and Sun--Abraham.
HC1 coverage is severely depressed (10--44%), while design-based coverage remains near 95%. The design effects are large (8--136$\times$), driven by the high weight variation in this scenario. The corrected practitioner heuristic --- weighted point estimate plus PSU-level clustering --- also maintains near-nominal coverage (93--97% across estimators), and is slightly conservative for CS-reg and CS-DR (96--97% versus the full design's 94%). The over-coverage of the cluster specification reflects the absence of stratification and finite-population reduction, both of which shrink the variance toward its design-consistent target. The substantive takeaway is the same as in Scenario 1: survey weighting is doing the estimand-correction work and PSU clustering is doing the variance work; strata and FPC tighten the variance further but are not required for nominal coverage.
Many federal surveys (BRFSS, NHANES, MEPS) are repeated cross-sections rather than panels. This scenario verifies that the design-based framework extends to this setting, using CS-reg with $\texttt{panel} = \texttt{False}$. Table (ref) reports the results.
Design-based coverage is 93.5--94.7% across all sample sizes, confirming that Proposition (ref) applies to repeated cross-sectional survey data. HC1 coverage degrades from 94.5% at $n = 500$ to 69.5% at $n = 8{,}000$, following the same pattern as the panel case in Scenario 1 (though with smaller design effects because the cross-sectional DGP has less persistent clustering).
This scenario examines whether survey-weighted doubly robust estimation --- where both the propensity score and outcome regression are estimated with survey weights --- produces valid inference under conditional parallel trends.
The DGP adds an X-dependent time trend: the outcome includes a term $\delta \cdot x_{1i} \cdot (t / T)$ where $\delta = 1.5$, and treated units' $x_1$ is drawn from $N(1, 1)$ rather than $N(0, 1)$. Because $\mathbb{E}[x_1 \mid \text{treated}] \neq \mathbb{E}[x_1 \mid \text{control}]$, the average time trend differs by group, violating unconditional parallel trends. Conditional on $x_1$, trends are identical --- conditional parallel trends holds. Table (ref) reports the results.
CS-reg, which does not condition on covariates, exhibits large and persistent bias ($+0.57$, roughly 29% of the true ATT) with 0% coverage at every sample size. The bias does not diminish with $n$ because it reflects a misspecification of the identifying assumption, not sampling variability.
CS-DR, which conditions on $x_1$ and $x_2$ via doubly robust estimation with survey-weighted nuisance models, eliminates this bias: mean bias is near zero at $n = 2{,}000$ with 94.7% coverage. The RMSE ratio between the two approaches is 8.2$\times$.
These simulation results provide strong evidence that the design-based variance framework extends to doubly robust estimation with survey-weighted nuisance models. A complete formal derivation would require showing that the santannazhao2020 influence function expansion carries through when the propensity score and outcome regression are estimated with survey-weighted estimating equations --- a nontrivial extension of their iid theory that we leave for future work. The simulation evidence here, with coverage at 92.7--94.7% across 2,000 replications and three sample sizes, suggests that the composition of survey-weighted nuisance estimation with design-based variance is well calibrated in this DGP. The practical implication is that practitioners working with survey data where parallel trends hold only after covariate adjustment can use design-based DR estimation with well-calibrated inference.
We illustrate the practical importance of design-based inference with an analysis of the Affordable Care Act's (ACA) dependent coverage provision using data from the National Health and Nutrition Examination Survey (NHANES).
The ACA's dependent coverage mandate, effective September 2010, allowed young adults to remain on their parents' health insurance until age 26. This created a natural experiment: adults aged 19--25 gained access to dependent coverage (treatment group), while adults aged 27--34, who were demographically similar but ineligible, serve as a comparison group. This age-cutoff DiD design was established by antwi2013 and has been widely used in evaluations of ACA coverage effects sommers2012.
Our outcome is binary health insurance coverage, estimated via a linear probability model. We use two NHANES cycles: 2007--2008 (pre-ACA) and 2015--2016 (post-ACA), yielding $n = 2{,}946$ observations (1,372 in the treatment group, 1,574 in the comparison group).
NHANES employs a stratified multi-stage cluster design to produce nationally representative health and nutrition data nchs2018. The public-use files provide masked survey design variables: 31 pseudo-strata (SDMVSTRA), 2 PSUs per stratum (SDMVPSU), and 2-year examination weights (WTMEC2YR). PSU identifiers are reused across strata, requiring the nested specification when constructing the survey design object. We follow CDC's analytic guidelines for proper variance estimation with these design variables.\footnote{The analytic sample is restricted to MEC-examined respondents (those with $\text{WTMEC2YR} > 0$), and we use WTMEC2YR as the appropriate weight for that domain. CDC's Weighting Module convention selects the sample weight by the data source of the smallest-respondent-set variable used; the MEC weight applies when the analysis is restricted to (or includes any variables from) the MEC examination component. An alternative analysis using WTINT2YR over the full interview sample would target a broader population --- adults who completed only the in-home interview as well as MEC-examined adults. We expect the qualitative lesson --- that weighting and clustering materially affect inference --- to remain, but this robustness check is left for future replication.}
Table (ref) presents DiD estimates under four specifications that progressively incorporate survey design features.
Three findings emerge. First, survey weighting changes the point estimate substantially: the unweighted ATT is 6.5 percentage points, while the survey-weighted ATT is 9.6 percentage points --- a 48% increase. This difference arises because NHANES's sampling design over-represents certain demographic groups. Without weighting, the estimated treatment effect reflects the sample composition rather than the U.S.\ population of young adults.
Second, ignoring PSU clustering understates uncertainty. Comparing the weights-only specification (SE $= 0.037$) with the full design specification (SE $= 0.044$), clustering inflates the standard error by 19% (a design effect of approximately 1.4$\times$ on the variance scale). This is modest compared to the simulation study in Section (ref) (where design effects reached 17$\times$ in the baseline scenario and exceeded 100$\times$ under informative sampling), reflecting the particular variance structure of this analytic sample, which has only 2 PSUs per stratum. Even so, the 95% confidence interval shifts from $[0.023, 0.170]$ to $[0.006, 0.187]$ --- the lower bound moves substantially closer to zero.
Third, the HC1 specification produces a qualitatively different conclusion: the unweighted, no-design ATT of 0.065 (SE $= 0.035$) has a 95% CI of $[-0.002, 0.133]$ that includes zero, while the design-based ATT of 0.096 (SE $= 0.044$) has a CI of $[0.006, 0.187]$ that excludes zero. This significance reversal is driven primarily by the change in the point estimate from proper survey weighting (0.065 to 0.096), not by the variance correction --- the design-based SE is actually larger than the HC1 SE, which widens the confidence interval. Both corrections --- weighting the point estimate upward and widening the confidence interval via design-based SEs --- are needed for valid population-level inference, though they have opposing effects on the test statistic. In this case, the point estimate correction dominates and flips the significance conclusion.
Covariate adjustment for gender and family income-to-poverty ratio reduces the ATT to 6.0 percentage points with a CI that includes zero, suggesting that some of the measured coverage gain is confounded with demographic shifts between the pre- and post-ACA periods. This underscores the importance of conditional parallel trends analysis in survey-weighted DiD --- exactly the setting addressed by Scenario 4 of the simulation study.
A caveat on the empirical design: our post-period (2015--2016) follows not only the 2010 dependent coverage mandate but also the major 2014 ACA provisions (Medicaid expansion and health insurance exchanges), which likely had differential effects on the treatment and control age groups. Standard evaluations of the dependent coverage mandate antwi2013 use post-periods prior to 2014 to avoid this confound. Our purpose here is to illustrate the methodological consequences of ignoring the survey design, not to make a causal claim about the dependent coverage mandate itself.
The theoretical framework and empirical results in this paper are implemented in diff-diff (v3.3.2; diffdiff2026), an open-source Python package for difference-in-differences causal inference. The package provides design-based variance estimation for 15 modern DiD estimators. To our knowledge, this is the first open-source implementation that jointly supports strata, PSU clustering, finite population corrections, and replicate-weight methods for the heterogeneity-robust estimators discussed in Section (ref).
The package follows an sklearn-like API where estimators are instantiated with configuration parameters and fitted to data via a fit() method. Survey design information is encapsulated in a SurveyDesign object that specifies the design variables and is passed to fit() as an optional argument. This design keeps the survey specification separate from the estimator logic, so the same estimator code works with or without survey adjustments:
Table (ref) summarizes the survey variance support for each estimator class. Regression-based estimators use the TSL sandwich (equation (ref)); IF-based estimators apply the stratified-cluster formula directly to influence function values (equation (ref)). Two estimators (Synthetic DiD, TROP) involve non-smooth optimization and rely on resampling for survey variance: Synthetic DiD supports Rao--Wu rescaled bootstrap, stratified placebo permutation, and PSU-level jackknife; TROP supports Rao--Wu rescaled bootstrap only.
diff-diff is available on PyPI (pip install diff-diff) and GitHub. Documentation, tutorials, and the survey tutorial that accompanies this paper are available at the project website.
This paper bridges the survey statistics and modern DiD literatures by showing that Binder's (1983) design-based variance theory applies directly to the influence functions of heterogeneity-robust DiD estimators. The result is simple in retrospect --- the regularity conditions are met, the formulas are known --- but it has not previously been stated for this class of estimators. Empirically, the framework recovers the practitioner heuristic that already sees wide informal use --- weighted point estimates with PSU-level clustering --- and identifies the additional precision available from the full strata/FPC specification.
The simulation results in Section (ref) resolve cleanly into a short set of recommendations. The first two steps are necessary; the third is a precision improvement on top of an already-valid specification.
\paragraph{Always weight the point estimate.} When the target is a population-level treatment effect, use survey weights for estimation regardless of which variance specification you choose. The NHANES illustration in Section (ref) showed a 48% difference in the point estimate between weighted and unweighted specifications --- a gap that no variance estimator can recover. Weights affect what the estimand is, not just how precisely it is measured. Under informative sampling, the bias from skipping weights is large (see Scenario 2, where unweighted estimators have RMSE roughly $1.1\times$ the weighted RMSE despite identical sample sizes).
\paragraph{Cluster at the PSU level for the variance.} For variance, PSU-level clustering through the influence-function representation --- a SurveyDesign with psu and survey weights but no strata or FPC --- is sufficient for nominal coverage in every scenario we examined, including under informative sampling. This matches the cluster=psu heuristic that many applied practitioners already use, and Section (ref) provides the formal justification that has been missing for the modern heterogeneity-robust estimators.
\paragraph{Add strata and FPC for precision, not coverage.} The full design specification (strata, PSU clustering, FPC, and weights) narrows the variance modestly relative to the cluster-only specification, but the gain is incremental. In our Scenario 2 (informative sampling, design effects up to $136\times$), the cluster-only specification was 1--3 percentage points more conservative than the full design --- the type of over-coverage that costs precision in these simulations but does not compromise validity. Practitioners who have access to the full design variables should use them. Practitioners working with public-use files where strata or FPC are masked or suppressed can use the cluster-only specification, which appears conservative for coverage in the designs we study; full design variables remain preferred when available, particularly when stratification is fine-grained or finite-population corrections are large.
\paragraph{Replicate weights are an equivalent path.} Some federal surveys (BRFSS, ACS, RECS) ship replicate-weight files designed to capture the full design variance without requiring strata and PSU identifiers, which are sometimes suppressed for confidentiality. Replicate-weight methods (BRR, Fay's BRR, JK1, JKn, SDR) applied to the influence-function values produce variance estimates asymptotically equivalent to the TSL formula. Practitioners with replicate weights should use them directly rather than reconstructing strata.
\paragraph{Use the design effect as a sanity check.} If the estimated DEFF (ratio of design-based to HC1 variance) is close to unity, the survey design has little impact and reporting HC1 alongside design-based SEs is reasonable. If DEFF exceeds two --- common in health and social surveys with intra-cluster correlation --- HC1 inference is materially misleading and design-based inference (cluster or full) should be primary.
Several limitations should be noted. The TSL variance formula requires at least two PSUs per stratum ($n_h \geq 2$); designs with singleton strata require special handling. The $t$-distribution approximation with $\text{df} = \sum_h n_h - H$ degrees of freedom may be anti-conservative when the total number of PSUs is small (yielding few overall degrees of freedom), a common concern in state-based surveys. Our simulation study used a balanced design ($H = 5$ strata, $n_h = 8$ PSUs per stratum) with a fixed number of PSUs as the sample size $n$ increases. This diverges from the standard PSU-based asymptotic sequence in which both $n$ and the number of PSUs grow together; our simulation instead increases the within-PSU sample size at fixed design complexity. The performance of design-based SEs under severely unbalanced designs, alternative asymptotic regimes, or very few PSUs per stratum merits further investigation. When the number of PSUs per stratum is small, finite-sample degrees-of-freedom corrections such as the bias-reduced linearization of bellmccaffrey2002 or the small-sample adjustments of faygraubard2001 can improve coverage; incorporating these into the design-based DiD framework is a natural next step.
A related practical issue is the level of clustering. The analytical framework aggregates to survey PSUs, but policy evaluations often apply treatments at higher geographic units (states, counties) that may cross-cut the survey stratification. When the treatment assignment level differs from the survey clustering level, reconciling the two sources of dependence requires additional care beyond the scope of this paper.
The framework assumes that the parallel trends assumption holds in the population. Survey weighting ensures that the estimand corresponds to the correct population, but it does not validate the identifying assumption itself. The conditional parallel trends result (Scenario 4) requires correct specification of at least one nuisance model, as with any doubly robust estimator.
Our empirical illustration uses a two-period, two-group design. While the theoretical framework covers the full staggered adoption case, the real-data example does not demonstrate this. Applications of design-based variance estimation to staggered treatment designs using federal survey data (e.g., state-level policy variation in the CPS or ACS) would be a valuable complement.
Several extensions merit investigation. First, calibration and post-stratification weights --- increasingly common in modern survey practice --- interact with the variance estimation framework in ways that are not fully addressed here. Second, multi-level treatment designs, where treatment is assigned at a level that cross-cuts the survey design (e.g., county-level policies in a state-stratified survey), raise additional questions about the appropriate variance estimator. Third, extending the analytical framework to non-smooth estimators such as synthetic control methods would broaden the applicability of design-based DiD inference. Finally, the growing availability of replicate weight files with public-use survey data creates an opportunity for practitioners to implement design-based variance without needing access to the underlying strata and PSU identifiers, which are sometimes suppressed for confidentiality.
The simulation pipeline (DGP, estimator wrappers, table and figure generation) and the NHANES illustration are implemented in Python against diff-diff v3.3.2. The replication code is publicly archived at \url{https://github.com/igerber/design-based-did-replication}, and the v1.0 release accompanying this preprint is deposited at Zenodo gerberreplication2026 under DOI \href{https://doi.org/10.5281/zenodo.20097361}{10.5281/zenodo.20097361}; the concept DOI \href{https://doi.org/10.5281/zenodo.20097360}{10.5281/zenodo.20097360} resolves to the latest version. The NHANES public-use files are available from the U.S.\ National Center for Health Statistics at \url{https://wwwn.cdc.gov/nchs/nhanes}; we ship the analysis-ready golden JSON used to produce Table (ref) so that the empirical illustration can be reproduced without re-downloading the raw NHANES cycles.