Skip to contents

Simulate Panel Data with Bad Controls

Source: R/simulate_data.R
simulate_bad_controls.Rd

Generates a staggered difference-in-differences panel dataset matching the Monte Carlo designs in Caetano, Callaway, Payne, and Sant'Anna (2024), 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

n

Number of units (default 2000)

T_max

Number of time periods (default 4)

groups

Integer vector of possible treatment-adoption periods, besides never-treated (default 2:T_max). Values must be between 2 and T_max, since every unit needs at least one pre-period.

dgp

Which counterfactual evolution equation for X_t(0) to use: "dgp1" (linear W), "dgp2" (nonlinear W), "dgp3" (nonlinear in (X(t-1), Z), no W – Simple Covariate Unconfoundedness holds), "dgp4" (linear W, coefficient 1 on the lag – parallel trends for the bad control given (W,Z) holds exactly, see Details), or "dgp5" (nonlinear W plus an \(X_{t-1} \times W\) interaction – something a linear model cannot represent at all, unlike a curved but additive term; see Details)

lambda

How much treatment shifts X at event time 0 (default 0.5)

delta

Direct effect of treatment on Y at event time 0, net of the effect transmitted through X (default 0.5)

kappa

Growth rate of the treatment effect with event time e = t - g; effects at event time e are (1 + kappa * e) times their event-time-0 value (default 0.5)

beta_drift

Drift rate of the loading on X(0) in the outcome equation across calendar time (default 0.2); see Details. Set to 0 to hold the loading fixed at 1 in every period.

binary_bad_control

Logical; if TRUE, the bad control X is binary instead of continuous, generated by treating the usual X(0) equations (baseline and dgp-specific evolution) as a logit index for a Bernoulli draw, and treatment as a shift in that index rather than a level shift in X. See Details.

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

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\). This package's own "nonlinear W" design – not (yet) a literal replica of the paper's own DGP2. The paper's original coefficients (\(0.15W_i^2\) added to dgp1's equation) 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 breaks 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 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\ linear term's own size at 1 SD of \(W\), versus the paper's 66\

  • 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 \(\ATT(g,t) = \beta_t \lambda_e + \delta_e\).

The true \(\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 \(\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.

Examples

sim <- simulate_bad_controls(n = 500, dgp = "dgp1")
head(sim$data)
#>   id period G D          Y          X         Z          W
#> 1  1      1 0 0 -1.1010691 -0.7336020 -1.675990 -1.3497984
#> 2  1      2 0 0 -1.2024615 -1.0611615 -1.675990 -1.3497984
#> 3  1      3 0 0 -1.0881218 -1.0855523 -1.675990 -1.3497984
#> 4  1      4 0 0 -1.4067053 -1.1984356 -1.675990 -1.3497984
#> 5  2      1 2 0  0.6013597  0.3189168  0.859543  0.1443711
#> 6  2      2 2 1  2.0100607  1.0625812  0.859543  0.1443711
sim$true_att_gt
#>   g t  att
#> 1 2 2 1.00
#> 2 2 3 1.65
#> 3 2 4 2.40
#> 4 3 3 1.10
#> 5 3 4 1.80
#> 6 4 4 1.20
sim$true_att_overall
#> [1] 1.525

# 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
#> [1] 1