Package {badcontrols}

badcontrols: Difference-in-Differences with Bad Controls

Description

Implements methods for difference-in-differences with bad controls, i.e., time-varying covariates that are affected by the treatment. Provides imputation, doubly robust, and machine learning estimators that are based on Caetano, Callaway, Payne, and Sant'Anna (2026) doi:10.48550/arXiv.2608.03881.

Author(s)

Maintainer: Hugo Sant'Anna hsantanna@uab.edu

Authors:

• Hugo Sant'Anna hsantanna@uab.edu

• Carolina Caetano ccaetano@uga.edu

• Brantly Callaway brantly.callaway@uga.edu

• Stroud Payne stroud.payne@vanderbilt.edu

See Also

Useful links:

• https://github.com/hugosantanna/badcontrols

• Report bugs at https://github.com/hugosantanna/badcontrols/issues

Estimate ATT(g,t) with Bad Controls

Description

Main function for difference-in-differences with bad controls. Wraps the ptetools infrastructure to handle time-varying covariates affected by treatment. Supports both imputation and doubly robust ML estimation methods.

Usage

didbc(
  yname,
  gname,
  tname,
  idname = NULL,
  data,
  bad_control_formula = NULL,
  xformula = ~1,
  d_covs_formula = NULL,
  bad_control_cov_formula = NULL,
  bad_control_d_cov_formula = NULL,
  est_method = c("imputation", "dr_ml"),
  bad_control_identification_strategy = c("unconfoundedness", "did"),
  nuisance_method = c("ml", "parametric"),
  control_group = "notyettreated",
  anticipation = 0,
  base_period = "varying",
  weightsname = NULL,
  cband = TRUE,
  alp = 0.05,
  boot_type = "multiplier",
  bstrap = TRUE,
  cl = 1,
  biters = 100,
  nfolds = 5,
  overlap_threshold = 0.99,
  min_group_size = 5,
  num_threads = NULL,
  ...
)

Arguments

Details

This function implements the methods from Caetano, Callaway, Payne, and Sant'Anna (2026). Two estimation approaches are available:

Imputation (est_method = "imputation"):

1. Among controls, learn X_t ~ f(X(t-1), W, Z)

2. For treated, predict counterfactual X_t(0)

3. Run DiD using imputed X_t(0) instead of observed X_t

Doubly Robust ML (est_method = "dr_ml"):

1. Estimate four nuisance functions via random forests

2. Combine into a doubly robust score with cross-fitting

3. Returns influence function for fast inference

Requires the grf package. Consistent if either the outcome regression or the propensity score is correctly specified.

Value

A pte_results object containing:

overall_att

Overall ATT estimate and SE

att_gt

Group-time specific ATT estimates

event_study

Dynamic (event-time) ATT estimates and SEs

ptep

Internal ptetools parameters object, used by summary()/plot()/autoplot() methods

References

Caetano, C., Callaway, B., Payne, S., and Sant'Anna, H. (2026). "Difference-in-Differences with Bad Controls." arXiv preprint arXiv:2608.03881. doi:10.48550/arXiv.2608.03881

Examples

# Simulate data
sim <- simulate_bad_controls(n = 500, T_max = 4)
head(sim$data)

# Imputation approach

res_imp <- didbc(
  yname = "Y", gname = "G", tname = "period", idname = "id",
  data = sim$data,
  bad_control_formula = ~X,
  xformula = ~Z,
  est_method = "imputation"
)

Doubly Robust Estimator for Bad Controls

Description

Implements the semiparametric doubly robust estimator from Caetano, Callaway, Payne, and Sant'Anna, with nuisance functions estimated either via machine learning or via parametric (linear/logit) working models, and K-fold cross-fitting. Estimates ATT(g,t) in the presence of post-treatment covariates (bad controls).

Usage

dr_ml_attgt(
  gt_data,
  xformula = ~1,
  bad_control_formula = NULL,
  d_covs_formula = NULL,
  bad_control_cov_formula = NULL,
  bad_control_d_cov_formula = NULL,
  nuisance_method = c("ml", "parametric"),
  bad_control_binary = FALSE,
  overlap_threshold = 0.99,
  min_group_size = 5,
  n_folds = 5,
  num_threads = NULL,
  ...
)

Arguments

Details

Implements Algorithm 1 from the paper: four nuisance functions, with K-fold cross-fitting:

• m_0(X_t*, X_t*-1, Z): outcome regression ⁠E[DeltaY | X_t*, X_t*-1, Z, D=0]⁠

• nu_0(X_t*-1, W, Z): ⁠E[m_0(X_t*, X_t*-1, Z) | X_t*-1, W, Z, D=0]⁠, a second-stage regression of m_0's fitted values on (X_t*-1, W, Z) among untreated units

• p_2(X_t*-1, W, Z): the propensity score P(D=1 | X_t*-1, W, Z)

• omega_0(X_t*, X_t*-1, Z): ⁠E[p_2/(1-p_2) | X_t*, X_t*-1, Z, D=0]⁠, a regression of the fitted propensity odds ratio on (X_t*, X_t*-1, Z) among untreated units

nuisance_method = "ml" (the default) estimates all four via cross-fitted random forests (grf package). nuisance_method = "parametric" estimates m_0/nu_0/omega_0 via OLS and p_2 via logistic regression; this assumes all four working models are correctly specified, in which case the doubly robust score's Neyman orthogonality means no further correction to the influence function is needed beyond what is already here for either choice of nuisance_method – parametric estimators converge faster than the ML rate conditions this orthogonality argument requires. Cross-fitting is retained under "parametric" for a single, shared code path, even though it is not required for validity there.

Unlike the imputation estimator, no counterfactual imputation of the bad control is needed: X_t* is observed directly for untreated units (X_t* = X_t*(0)), and nu_0/omega_0 target m_0's/p_2's fitted values through a second regression rather than plugging in a predicted X_t*(0).

nu_0 and omega_0 are themselves nested conditional expectations of m_0/p_2 (e.g. nu_0 = E[m_0 | X_t*-1, W, Z, D=0]), so regressing m_0's/ p_2's fitted values on the coarser feature set using the same training-fold observations that fit m_0/p_2 would be a generated- regressor problem: the pseudo-outcome's estimation error is correlated with the very sample the nested regression is fit on, not just small in L^2. Under nuisance_method = "ml", this is avoided using each forest's out-of-bag (OOB) predictions – grf's predict(fit) with no newdata, which averages only over trees that did not include a given row in their subsample – as the pseudo-outcome for nu_0/omega_0's training data, instead of the in-sample predict(fit, newdata = <training data>). This costs nothing extra to fit (no additional forests, unlike an explicit training-fold split) and keeps the full training fold available to nu_0/omega_0, at the cost of a less clean-cut independence argument than literal disjoint subsamples would give. The main m_0/p_2 fits used directly in the doubly robust score are unaffected – they are already evaluated out-of-sample via the outer K-fold. Under nuisance_method = "parametric", no OOB equivalent applies (lm/glm have no such notion); that path is unchanged, per the paragraph above.

The doubly robust score is consistent if either (m_0, nu_0) or (p_2, omega_0) are correctly specified.

With no bad control at all (bad_control_formula = NULL), nu_0 and omega_0 are not estimated at all: nu_0 = m_0 and omega_0 = p_2/(1-p_2) exactly (m_0 no longer depends on X_t*, so there is nothing left for nu_0 to marginalize over; p_2 is already a function of Z alone, so omega_0's further conditioning on Z is a no-op). The doubly robust score then reduces exactly to the classical Sant'Anna and Zhao (2020) AIPW-DiD estimator.

There are two triggers that fall back to imputation_attgt() for a whole (g,t) cell, rather than dropping p_2/omega_0 alone – doing that would leave a moment that is not Neyman orthogonal, which is exactly why imputation_attgt() needs its own first-stage correction terms that this function's nuisances don't have:

1. Small treated group. If the cell has fewer treated units than the number of p_2 covariates plus min_group_size, a warning names the group, time period, and treated count. Fitting p_2 needs a nontrivial treated sample to identify at all, relative to how many covariates it has to fit (a logit can hit perfect separation, a probability_forest is essentially memorizing a handful of points), which can happen even when the population-level covariate distributions overlap fine – a distinct problem from overlap below.

2. Overlap. Before cross-fitting, a preliminary propensity model (not used in the final estimation) is fit once on the whole cell: if any unit's fitted propensity exceeds overlap_threshold, a warning names the group, time period, and offending unit IDs.

Value

attgt_if object with ATT estimate and influence function

References

Caetano, C., Callaway, B., Payne, S., and Sant'Anna, H. (2026). "Difference-in-Differences with Bad Controls." arXiv preprint arXiv:2608.03881. doi:10.48550/arXiv.2608.03881

Sant'Anna, P.H.C. and Zhao, J. (2020). "Doubly Robust Difference-in-Differences Estimators." Journal of Econometrics.

Wager, S. and Athey, S. (2018). "Estimation and Inference of Heterogeneous Treatment Effects using Random Forests." Journal of the American Statistical Association.

Examples

## Not run: 
# dr_ml_attgt() is normally called internally as didbc()'s attgt_fun
# (via ptetools::pte()), not invoked directly. This is the standard way
# to reach it end-to-end; slow due to cross-fitted random forests plus
# the multiplier bootstrap, hence \dontrun.
sim <- simulate_bad_controls(n = 500)
res <- didbc(
  yname = "Y", gname = "G", tname = "period", idname = "id",
  data = sim$data,
  bad_control_formula = ~X,
  xformula = ~Z,
  est_method = "dr_ml",
  nuisance_method = "parametric"
)
summary(res)

## End(Not run)

Imputation ATT(g,t) estimator for bad controls

Description

Imputation ATT(g,t) estimator for bad controls

Usage

imputation_attgt(
  gt_data,
  xformula = ~1,
  bad_control_formula = NULL,
  d_covs_formula = NULL,
  bad_control_cov_formula = NULL,
  bad_control_d_cov_formula = NULL,
  bad_control_binary = FALSE,
  bad_control_identification_strategy = c("unconfoundedness", "did"),
  ...
)

Arguments

Details

This function only handles the shared leadin (pivoting gt_data to one row per unit, and constructing the covariate columns/names), then dispatches to imputation_unconfoundedness or imputation_did depending on bad_control_identification_strategy.

Imputation estimator under parallel trends for the bad control

Description

Implements the two-step imputation estimator and its influence function under parallel trends for the bad control itself, rather than Covariate Unconfoundedness: ass:bad-control-parallel-trends and cor:att-under-bad-control-parallel-trends-and-linearity in app:bad-control-parallel-trends of the supplementary appendix, and dev/bad_control_parallel_trends_influence_function.md for the influence function derivation. Step 1 regresses the bad control's own change, rather than its post-period level, on (W, Z) among the comparison group; Step 2 (the outcome regression) is unchanged from imputation_unconfoundedness. Called by imputation_attgt when bad_control_identification_strategy = "did".

Usage

imputation_did(
  wide_data,
  pre_data,
  post_data,
  D,
  n,
  x_names,
  dx_names,
  bc_cov_names,
  bc_dcov_names,
  comparison_idx,
  bad_control_formula
)

Arguments

Imputation estimator under Covariate Unconfoundedness for the bad control

Description

Implements the two-step imputation estimator and its influence function (eq:psi-ra-final in the supplementary appendix) under Covariate Unconfoundedness: the bad control's untreated potential evolution is modeled from (bc_pre, W, Z) among the comparison group. Called by imputation_attgt when bad_control_identification_strategy = "unconfoundedness".

Usage

imputation_unconfoundedness(
  wide_data,
  pre_data,
  post_data,
  D,
  n,
  x_names,
  dx_names,
  bc_cov_names,
  bc_dcov_names,
  comparison_idx,
  bad_control_formula,
  bad_control_binary
)

Arguments

NLSY79 Job Displacement Application Data

Description

A balanced panel of 3,231 NLSY79 respondents observed biennially from 1992 through 2002. The sample applies the paper's positive-earnings restriction and excludes respondents first treated in 1992, who have no pre-treatment period within the application window. The outcome is log earnings, the bad control is an occupation-based wage score, and treatment is job displacement.

Usage

data(nlsy_job_displacement)

Format

A data.frame with 19,386 rows and 8 columns: id, year, log_earnings, occ_score, group, race, female, and educ_max_grade.

Details

Job displacement is defined as involuntarily leaving a job because of a layoff/job elimination or plant/company/workplace closure. The occupation score is time-varying at the individual level because respondents can change occupations, even though the score is fixed for a given occupation. Respondents first displaced in 1992 are excluded because the application requires a pre-treatment period within the 1992–2002 window.

This is a researcher-created subset of public-use NLSY79 data and an occupation-score merge based on the public-use IPUMS USA 1990 5\ It is not an official NLSY79 data release. See the source files in ⁠data-raw/⁠ for the construction script and provenance.

Source

National Longitudinal Survey of Youth 1979 (NLSY79), U.S. Bureau of Labor Statistics, https://www.nlsinfo.org/investigator/; IPUMS USA, https://usa.ipums.org/usa/.

References

Bureau of Labor Statistics, U.S. Department of Labor. National Longitudinal Survey of Youth 1979 cohort, 1979-2022 (rounds 1-30). Produced and distributed by the Center for Human Resource Research (CHRR), The Ohio State University. Columbus, OH.

Steven Ruggles, Sarah Flood, Matthew Sobek, Daniel Backman, Grace Cooper, Julia A. Rivera Drew, Stephanie Richards, Renae Rogers, Jonathan Schroeder, and Kari C.W. Williams. IPUMS USA: Version 16.0 [dataset]. Minneapolis, MN: IPUMS, 2025. https://doi.org/10.18128/D010.V16.0

Examples

data(nlsy_job_displacement)
head(nlsy_job_displacement)
table(nlsy_job_displacement$group[!duplicated(nlsy_job_displacement$id)])

Simulate Panel Data with Bad Controls

Description

Generates a staggered difference-in-differences panel dataset matching the Monte Carlo designs in Caetano, Callaway, Payne, and Sant'Anna (2026), where a time-varying covariate X is affected by treatment (a bad control). Returns the known group-time, event-study, and overall ATT alongside the data for testing estimators against.

Usage

simulate_bad_controls(
  n = 2000,
  T_max = 4,
  groups = 2:T_max,
  dgp = c("dgp1", "dgp2", "dgp3", "dgp4", "dgp5"),
  lambda = 0.5,
  delta = 0.5,
  kappa = 0.5,
  beta_drift = 0.2,
  binary_bad_control = FALSE
)

Arguments

Details

Common structure (shared across all dgp choices): Z_i,\eta_i \sim N(0,1) and W_i = 0.8\eta_i + 0.3Z_i + 0.2\varepsilon^W_i are time-invariant unit characteristics; \eta_i is unobserved. Treatment group is assigned by splitting units into equal-sized bins of the latent index 0.2Z_i + 0.4W_i + 0.3\eta_i + \varepsilon^D_i, with the lowest bin never-treated and successive bins assigned to groups in increasing order.

X_{i1} = 0.5\eta_i + 0.4Z_i + 0.3\varepsilon^{X_1}_i, and for t \ge 2, X_{it}(0) evolves according to the dgp-specific equation (redrawing fresh noise every period):

• dgp1: 0.7 X_{i,t-1}(0) + 0.3Z_i + 0.2W_i + 0.15. Matches the paper's actual Monte Carlo design, including the "Exclude BC" finding that dropping the bad control entirely stays nearly unbiased there specifically.

• dgp2: 0.7 X_{i,t-1}(0) + 0.3Z_i + 0.2W_i + 0.03 W_i^2 + 0.15. The "nonlinear W" design, matching the paper's own DGP2 exactly. An earlier version of the paper's DGP2 used a larger quadratic coefficient (0.15W_i^2 added to dgp1's equation), which put the quadratic term's vertex at W = -0.667, well inside W's support (SD(W) \approx 0.88), making the map from W to X_{it}(0) two-to-one for nearly the whole population – a pathological non-monotonicity, not just "nonlinearity," that specifically broke dr_ml_attgt()'s \omega_0 (which conditions on (X_{it}, X_{i,t-1}, Z), not W, and so must implicitly marginalize over a non-invertible W). This package's dgp2 keeps the same linear coefficient on W as dgp1 (0.2, so it still reads as "dgp1 plus one added wrinkle"), but with a small enough quadratic coefficient (0.03) that the vertex sits at W = -3.33 (about 3.8 SDs out) – safely beyond what any realistic sample reaches, while still contributing a real, detectable nonlinearity (about 13\ W). The paper's own DGP2 was subsequently redesigned to use this same 0.03W_i^2 coefficient, so the two now match exactly.

• dgp3: 0.7 X_{i,t-1}(0) + 0.3Z_i + 0.4 X_{i,t-1}(0) Z_i + 0.2 X_{i,t-1}(0)^2 + 0.15 (no W). Matches the paper's actual Monte Carlo design.

• dgp4: X_{i,t-1}(0) + 0.3Z_i + 0.2W_i + 0.15. The coefficient of 1 on the lag (a random walk with drift in (Z,W)), unlike dgp1's 0.7, is deliberate: it makes \Delta X_{it}(0) = X_{it}(0) - X_{i,t-1}(0) itself a function of (Z,W) alone (plus independent noise), so parallel trends for the bad control given (W,Z) (app:bad-control-parallel-trends in the supplementary appendix) holds exactly under dgp4, alongside Covariate Unconfoundedness given (X_{t-1},W,Z). Does not attempt to match any Monte Carlo design in the paper, and exists only to test the parallel-trends identification strategy.

• dgp5: 0.7 X_{i,t-1}(0) + 0.3Z_i + 0.2W_i + 0.03 W_i^2 + 0.05 X_{i,t-1}(0) W_i + 0.15. dgp2's equation plus an X_{t-1} \times W interaction. Unlike a curved-but-additive term, an interaction is something a linear model (as in imputation_attgt()'s Step 1) cannot represent at all, no matter how strong – a categorical, not just approximate, misspecification. Meant to show a bigger separation between Imputation and dr_ml_attgt() than dgp2 does, since forests capture interactions natively without needing to be told where to look. The interaction coefficient (0.05) was chosen empirically: it's the smallest of several tried that still gives Imputation clearly larger (roughly double or more) bias than dgp2 does, while keeping dr_ml_attgt()'s own analytical-SE calibration reasonable (larger coefficients tried, e.g. 0.10, gave a starker bias gap but degraded nuisance_method = "ml"'s SE calibration, apparently worsening rather than improving with n – a real, persistent problem rather than finite-sample noise, unlike dgp5's own milder behavior, which does improve with n).

each plus 0.3\varepsilon^{X_t}_i.

The untreated outcome is Y_{it}(0) = 0.3t + 0.5\eta_i + 0.3Z_i + \beta_t X_{it}(0) + 0.3\varepsilon^Y_{it}, with a time-varying loading \beta_t = 1 + \code{beta\_drift}(t-2) (equal to 1 at t=2 regardless of beta_drift). Once treated (period t \ge G_i, event time e = t - G_i), the observed covariate and outcome are X_{it} = X_{it}(0) + \lambda_e and Y_{it} = Y_{it}(0) + \beta_t \lambda_e + \delta_e, where \lambda_e = \lambda(1+\kappa e) and \delta_e = \delta(1+\kappa e). The \beta_t \lambda_e term routes the treatment's effect on X through the (time-varying) X-Y relationship, giving true \mathrm{ATT}(g,t) = \beta_t \lambda_e + \delta_e.

The true \mathrm{ATT}(g,t) only ever depends on \beta_t at the post-treatment period t, which equals 1 whenever t=2 regardless of beta_drift – so with the defaults, T_max = 2, groups = 2 already gives true \mathrm{ATT} = \delta+\lambda = 1.00. But beta_drift != 0 still makes \beta_1 \ne \beta_2, so an estimator that assumes a constant X-Y loading across periods (as linear Imputation implicitly does) will show bias in that case even though the true ATT value matches the paper. To reproduce the paper's two-period designs exactly, including the constant-loading assumption, set beta_drift = 0 as well.

When binary_bad_control = TRUE, the same equations for X_i1 and X_it(0) (dropping their additive noise term) are used as a logit index instead of a direct value: X_{i1} \sim \text{Bernoulli}(\Lambda(0.5 \eta_i + 0.4 Z_i)) and X_{it}(0) \sim \text{Bernoulli}(\Lambda(\text{evolve}_x)), where \Lambda is the logistic CDF. Treatment shifts the index rather than the level: X_{it}(1) \sim \text{Bernoulli}(\Lambda(\text{evolve}_x + \lambda_e)). This reuses each dgp's existing evolution equation unchanged, so dgp1 + binary_bad_control = TRUE is the case where the imputation estimator's binary Step 1 (logistic regression) is correctly specified, the same role dgp1 plays for the continuous, OLS Step 1.

Value

A list with:

data

Panel data.frame with columns id, period, G, D, Y, X, Z, W

true_att_gt

data.frame(g, t, att): the true ATT(g,t) for every valid group-time cell. Computed exactly from the DGP parameters when binary_bad_control = FALSE; when TRUE, there is no closed form (the shift in E(X) varies by unit), so this is instead the realized average over this sample's group-g cohort, same as true_att_by_e/true_att_overall below

true_att_by_e

data.frame(e, att): the true event-study ATT at each realized event time, averaged over the (g,t) cells that share that event time in this sample

true_att_overall

Realized sample average of the true individual effect across every treated (i,t) observation in this sample

Examples

sim <- simulate_bad_controls(n = 500, dgp = "dgp1")
head(sim$data)
sim$true_att_gt
sim$true_att_overall

# Collapses exactly to the paper's two-period design
sim2 <- simulate_bad_controls(n = 500, T_max = 2, groups = 2, dgp = "dgp1",
                               beta_drift = 0)
sim2$true_att_overall