gen_data <- function( n = 100,
R2 = 0.5,
beta_0 = 0, beta_1 = 1,
tau_0 = 1, tau_1 = 1,
alpha = 1, df = 3 ) {
stopifnot( R2 >= 0 && R2 < 1 )
X_i = rexp( n, rate = 1 ) - 1
beta_1 = sqrt( 1 - R2 )
sigma_e = sqrt( R2 )
Y0_i = beta_0 + beta_1 * X_i + rnorm( n, sd=sigma_e )
tau_i = tau_0 + tau_1 * X_i + alpha * rt( n, df = df )
Y1_i = Y0_i + tau_i
tibble( X = X_i, Y0 = Y0_i, Y1 = Y1_i )
}20 Simulation under the Potential Outcomes Framework
In this chapter we give some general approaches for running simulations in the context of causal inference. The broader point of this chapter is to underscore that one can be very creative in how one designs a simulation, insofar as there is a well-defined way to make data with a well-defined “answer” that we are trying to recover by analyzing the synthetic data.
When evaluating causal inference methods, we would probably turn to the potential outcomes framework. The potential outcomes framework is a framework typically used in the causal inference literature to make explicit statements regarding the mechanics of causality and the associated estimands one might target. While we recommend reading, for a more thorough overview, either Rosenbaum (2017) or Gerber and Green (2017), we briefly outline this framework here to set out our notation.
Take a sample of \(n\) experimental units, indexed by \(i\). For each unit, we can treat it or not, with \(Z_i = 1\) if we do and \(Z_i = 0\) if we do not. Now we imagine each unit has two potential outcomes we might see depending on whether we treated it (\(Y_i(1)\)) or not (\(Y_i(0)\)). We observe only one of these outcomes, depending on what we did: \[ Y_i^{obs} = Z_i Y_i(1) + (1-Z_i)Y_i(0) .\] The “Fundamental Problem of Causal Inference” is that we can never see both potential outcomes for a given unit. But with simulation, since we will be making our own data, we can—and that allows us to evaluate causal inference methods in many natural-seeming scenarios.
Each unit’s treatment effect is then \(\tau_i = Y_i(1) - Y_i(0)\); it is how much our outcome changes if we treat vs. not treat. Frustratingly, because we can only see one of a unit’s two potential outcomes, we can never directly estimate an individual \(\tau_i\). Under this view, causality is a missing data problem: if we only were able to impute the missing potential outcomes, we could have a dataset where we could calculate any estimands we wanted.
E.g., the most common estimand of the true average treatment effect for the sample \(\mathcal{S}\) would be:
\[ ATE_{\mathcal{S}} = \frac{1}{N} \sum_{i} Y_i(1) - Y_i( 0 ) . \] The average proportional increase, by contrast, would be
\[ API_{\mathcal{S}} = \frac{1}{N} \sum_{i} \frac{Y_i(1)}{Y_i(0)} \] Importantly, under a potential outcomes framework, our sample is defined by the potential outcomes, and randomness comes from the random assignment of treatment, not the sample itself. We call the process of how units are assigned to treatment the assignment mechanism.
20.1 Finite vs. Superpopulation inference
In the above we consider the sample as fixed. That said, if we thought of the sample as being drawn from some larger population (often called a “superpopulation”), we could talk about the true ATE of that larger population.
This is a tension that often arises in potential outcomes based simulations: if we are focused on \(ATE_{\mathcal{S}}\) then for each sample we generate, and different units have different \(\tau_i\), then our estimand could be different, depending on whether our sample has more or fewer units with high \(\tau_i\). If, on the other hand, we are focused on where the units came from (which is our data generating model), our estimand is a property of the DGP, and would be the same across our samples.
Now if we were generating data with a constant treatment impact, then all the \(\tau_i\) are the same and \(ATE_{\mathcal{S}} = ATE\) always; this is typical for many simulations in the literature. That being said, treatment variation is what causes a lot of methods to fail, and so having simulations with treatment variation is usually important.
The catch is when we calculate our performance metrics, we now have two possible targets to pick from: do we compare our estimate to the sample’s true treatment effect, or the population’s? Furthermore, if we are targeting the population ATE, then our error in estimation may be due in part to the representativeness of the sample, not the estimation or uncertainty due to the random assignment.
We will follow this theme throughout this chapter.
20.2 Generating data under the potential outcome framework
When writing a simulation using the potential outcomes framework, we recommend first generating a complete set of potential outcomes, then a random assignment based on some assignment mechanism, and finally the observed outcomes as a function of assignment and original potential outcomes.
For the first step, we generate each unit \(i = 1, \ldots, n\), as \[ \begin{aligned} X_i &\sim exp( 1 ) - 1 \\ Y_i(0) &= \beta_0 + \beta_1 X_i + \epsilon_i \mbox{ with } \epsilon_i \sim N( 0, \sigma^2 ) \\ \tau_i &= \tau_0 + \tau_1 X_i + \alpha u_i \mbox{ with } u_i \sim t_{df} \\ Y_i(1) &= Y_i(0) + \tau_i , \end{aligned} \] with \(exp(1)\) being the standard exponential distribution and \(t_{df}\) being a \(t\) distribution with \(df\) degrees of freedom. We subtract 1 from \(X_i\) to zero-center it (it is often convenient to have zero-centered covariates so we can then, e.g., interpret \(\tau_0\) as the true superpopulation ATE of our experiment).
The above model says we first, for each unit, generate a covariate. We then generate our two potential outcomes. I.e., we are generating what the outcome would be for each unit if it were treated and if it were not treated. We are driving both the level and the treatment effect with \(X_i\).
The result of the above is a schedule of potential outcomes for our \(n\) units. The schedule of potential outcomes encodes, for each unit, what would happen to it given the possible treatment assignements it might receive.
One advantage of generating all the potential outcomes is we can then calculate the finite-sample estimands such as the true average treatment effect for the generated sample: we just take the average of \(Y_i(1) - Y_i(0)\) for our sample.
Here is some code implementing the above model (we leave treatment assignment to later):
We have reparameterized our model so we can index our DGP by an R2 value rather than a coefficient on X so we can have a standardized control-side outcome (the expected variation of \(Y_i(0)\) will be 1). The treatment outcomes will have more variation due to the heterogeneity of the treatment impacts.
Here is some sample data:
set.seed( 40454 )
d1 <- gen_data( 50 )
head( d1 )# A tibble: 6 × 3
X Y0 Y1
<dbl> <dbl> <dbl>
1 -0.769 -0.572 2.50
2 0.136 -0.549 1.95
3 -0.164 -1.12 2.21
4 0.417 0.572 0.623
5 -0.623 0.493 -2.77
6 -0.128 0.0174 0.277We store latent quantities (both potential outcomes, in particular) so we can calculate the estimands of interest or learn about our data generating process. When we hand the data to an estimator, we would not provide this “secret” information.
Note how our estimand can change with each sample generated:
mean( d1$Y1 - d1$Y0 )[1] 0.6374925d2 <- gen_data( 50 )
mean( d2$Y1 - d2$Y0 )[1] 0.5479788In reviewing our model and code, we know our superpopulation ATE should be tau, or 1 exactly. If our estimate for d1 is 0.6 do we say that is close or far from the target? From a finite sample performance approach, we nailed it. From superpopulation, less so.
Once we have our schedule of potential outcomes, we generate the observed outcomes by assigning our (synthetic, randomly generated) \(n\) units to treatment or control. For example, say we wanted to simulate an observational context where chance of treatment was a function of our covariate. We could model each unit as flipping a weighted coin with some probability that was a function of \(X_i\) as so:
\[ \begin{aligned} p_i &= logit^{-1}( \xi_0 + \xi_1 X_i ) \\ Z_i &= Bern( p_i ) \\ Y_i &= Z_i Y_i(1) + (1-Z_i) Y_i(0) \end{aligned} \]
Here is code for assigning our data to treatment and control:
assign_data <- function( dat,
xi_0 = -1, xi_1 = 1 ) {
n = nrow(dat)
dat = mutate( dat,
p = arm::invlogit( xi_0 + xi_1 * X ),
Z = rbinom( n, 1, prob=p ),
Yobs = ifelse( Z == 1, Y1, Y0 ) )
dat
}We can then add our assignment variable to our given data as so:
assign_data( d2 )# A tibble: 50 × 6
X Y0 Y1 p Z Yobs
<dbl> <dbl> <dbl> <dbl> <int> <dbl>
1 0.670 0.667 2.58 0.418 1 2.58
2 0.371 0.314 4.57 0.348 1 4.57
3 1.94 1.29 3.03 0.719 0 1.29
4 -0.244 0.119 -10.0 0.224 1 -10.0
5 0.00850 1.44 2.88 0.271 0 1.44
6 1.41 1.14 5.02 0.600 1 5.02
7 -0.864 0.461 0.802 0.134 1 0.802
8 -0.00533 -0.914 -1.17 0.268 0 -0.914
9 -0.907 -0.202 0.555 0.129 1 0.555
10 -0.363 -0.141 1.16 0.204 1 1.16
# ℹ 40 more rowsNote how Yobs is, depending on Z, either Y0 or Y1. Separating our our DGP and our random assignment underscores the potential outcomes framework adage of the data are what they are, and we the experimenters (or nature) are randomly assigning these whole units to various conditions and observing the consequences.
We might put the p_i part of the model in our code generating the outcomes, if we wanted to view the chance of treatment assignment as inherent to the unit (which is what we usually expect in an observational context).
20.3 Finite sample performance measures
For a finite sample evaluation, we can generate a single dataset with our DGP from above, and run a small experiment where we actually randomize units to treatment and control:
n = 100
set.seed(442423)
dat = gen_data(n, tau_1 = -1)
dat = mutate( dat,
Z = 0 + (sample( n ) <= n/2),
Yobs = ifelse( Z == 1, Y1, Y0 ) )
mod = lm( Yobs ~ Z, data=dat )
coef(mod)[["Z"]][1] 0.8914992We can compare this to the true finite-sample ATE:
tauS = mean( dat$Y1 - dat$Y0 )
tauS[1] 1.154018Our finite-population simulation would be:
rps <- simhelpers::repeat_and_stack( 1000, {
dat = mutate( dat,
Z = 0 + (sample( n ) <= n/2),
Yobs = ifelse( Z == 1, Y1, Y0 ) )
mod = lm( Yobs ~ Z, data=dat )
tibble( ATE_hat = coef(mod)[["Z"]],
SE_hat = arm::se.coef(mod)[["Z"]] )
})
rps %>%
summarise( EATE_hat = mean( ATE_hat ),
bias = EATE_hat - tauS,
SE = sd( ATE_hat ),
ESE_hat = sqrt( mean( SE_hat^2 ) ) )# A tibble: 1 × 4
EATE_hat bias SE ESE_hat
<dbl> <dbl> <dbl> <dbl>
1 1.16 0.00482 0.248 0.309We are simulating on a single dataset. In particular, our set of potential outcomes is entirely fixed; the only source of randomness (and thus the randomness behind our SE) is the random assignment. The finite sample simulation opens up some room for critique: what if our single dataset is non-standard? How can we claim our findings are general? For example, perhaps the difference in the true SE and average estimated one is a quirk of our particular dataset, and not a general phenomenon.
Our super-population simulation would be, by contrast:
rps_sup <- simhelpers::repeat_and_stack( 1000, {
dat = gen_data(n)
dat = mutate( dat,
Z = 0 + (sample( n ) <= n/2),
Yobs = ifelse( Z == 1, Y1, Y0 ) )
mod = lm( Yobs ~ Z, data=dat )
tibble( ATE_hat = coef(mod)[["Z"]],
SE_hat = arm::se.coef(mod)[["Z"]] )
})
rps_sup %>% summarise( EATE_hat = mean( ATE_hat ),
bias = EATE_hat - 1,
SE = sd( ATE_hat ),
ESE_hat = sqrt( mean( SE_hat^2 ) ) )# A tibble: 1 × 4
EATE_hat bias SE ESE_hat
<dbl> <dbl> <dbl> <dbl>
1 1.000 -0.00000591 0.381 0.389First, note our superpopulation simulation is not biased for the superpopulation ATE. Also note the true SE is larger than our finite-sample simulation; this is because part of the uncertainty in our estimator is the uncertainty of whether our sample is representative of the superpopulation.
Finally, this clarifies that our linear regression estimator is estimating standard errors assuming a superpopulation model. The true finite sample standard error is less than the expected estimated standard error: from a finite sample perspective, our estimator is giving overly conservative uncertainty estimates. (This discrepancy is often called the correlation of potential outcomes problem.)
20.4 Nested finite simulation procedure
We just saw a difference between a specific, single, finite-sample dataset and a superpopulation. What if we wanted to know if this phenomenon was more general across a set of finite datasets? This question can be levied more broadly: if we run a simulation on a single dataset, this is even more narrow than running on a single scenario: if we compare methods and find one is superior to another for our single dataset, how do we know this is not an artifact of some specific characteristic of that data and not a general phenomenon at all?
One way forward is to run a nested simulation, where we generate a series of finite sample datasets, and then for each dataset run a small simulation. We then calculate the expected finite sample performance across the datasets. One could almost think of the datasets themselves as a “factor” in our multifactor experiment. See Miratrix et al. (2021) for an example of this kind of simulation framework.
Borrowing from the simulation appendix of Miratrix et al. (2021), for a nested simulation repeat \(R\) times:
Generate a dataset using a particular DGP. This data generation is the “sampling step” for a superpopulation (SP) framework. The DGP represents an infinite superpopulation. We store both potential outcomes for all generated units.
Record the true finite-sample ATE for our generated data.
Then, three times, do a finite simulation as follows:
- Randomize units to treatment and control using the assignment mechanism.
- Calculate the corresponding observed outcomes.
- Analyze the results using the methods of interest, recording both the point estimate and estimated standard error for each.
Having only three trials will give a poor estimate of within-dataset (finite sample) performance, but the average performance across the \(R\) datasets in a given scenario will give a reasonable estimate of the expected performance across datasets of the type we would see given the scenario parameters.
To demonstrate we first make a mini-finite sample driver:
one_finite_run <- function( R0 = 3, n = 100, ... ) {
dat = gen_data( n = n, ... )
rps <- simhelpers::repeat_and_stack( R0, {
dat = mutate( dat,
Z = 0 + (sample( n ) <= n/2),
Yobs = ifelse( Z == 1, Y1, Y0 ) )
mod = lm( Yobs ~ Z, data=dat )
tibble( ATE_hat = coef(mod)[["Z"]],
SE_hat = arm::se.coef(mod)[["Z"]] )
}) %>%
bind_rows()
rps$ATE = mean( dat$Y1 - dat$Y0 )
rps
}This driver also stores the finite sample ATE for future reference:
one_finite_run()# A tibble: 3 × 3
ATE_hat SE_hat ATE
<dbl> <dbl> <dbl>
1 0.348 0.421 0.768
2 1.32 0.472 0.768
3 1.17 0.549 0.768We then run a bunch of finite runs.
runs <- simhelpers::repeat_and_stack( 500, one_finite_run(),
id = "runID" )
head( runs ) runID ATE_hat SE_hat ATE
1 1 0.61007830 0.4908961 0.7583748
2 1 0.04792879 0.3968365 0.7583748
3 1 0.76594448 0.4617217 0.7583748
4 2 0.62535954 0.3123307 0.8821758
5 2 0.59101011 0.3247807 0.8821758
6 2 1.04007072 0.3666993 0.8821758We set an id because we will need to separate out each finite run and analyze separately, and then aggregate.
Each finite run is a very noisy simulation for a fixed dataset. This means when we calculate performance measures we have to be careful to avoid bias in the finite-level calculations; in particular, we need to focus on estimating \(SE^2\) across the finite runs, not \(SE\), to avoid the bias caused by having a few replicates within each finite sample run:
fruns <- runs %>%
group_by( runID ) %>%
summarise( EATE_hat = mean( ATE_hat ),
SE2 = var( ATE_hat ),
ESE_hat = mean( SE_hat ),
.groups = "drop" )And then we aggregate our finite sample runs:
res <- fruns %>%
summarise( EEATE_hat = mean( EATE_hat ),
EESE_hat = sqrt( mean( ESE_hat^2 ) ),
ESE = sqrt( mean( SE2 ) ) ) %>%
mutate( calib = 100 * EESE_hat / ESE )
res# A tibble: 1 × 4
EEATE_hat EESE_hat ESE calib
<dbl> <dbl> <dbl> <dbl>
1 0.996 0.380 0.331 115.We see our expected standard error estimate is, across the collection of finite sample scenarios all sharing a similar parent superpopulation DGP, 15% too large for the true expected finite-sample standard error.
20.5 Superpopulation performance measures
For a superpopulation simulation, the easy road is to simply generate a new dataset with each simulation trial. We can also use our collection of mini-finite-sample runs to estimate superpopulation quantities as well. Given that the simulation datasets are i.i.d. draws, we can simply take expectations across all our simulations. The only concern is our estimates of MCSE will be off due to the clustering in our simulation runs.
Here we calculate superpopulation performance measures (both with the squared SE and without; we prefer the squared version for assessing calibration):
runs %>%
summarise( EATE_hat = mean( ATE_hat ),
SE_true = sd( ATE_hat ),
SE_hat = mean( SE_hat ),
SE2_true = var( ATE_hat ),
SE2_hat = mean( SE_hat^2 ) ) %>%
pivot_longer( cols = c(SE_true:SE2_hat ),
names_to = c( "estimand", ".value" ),
names_sep ="_" ) %>%
mutate( inflate = 100 * hat / true )# A tibble: 2 × 5
EATE_hat estimand true hat inflate
<dbl> <chr> <dbl> <dbl> <dbl>
1 0.996 SE 0.389 0.377 96.9
2 0.996 SE2 0.151 0.142 93.9We see our estimated SEs may be a bit too small, on average, as compared to the true SE.
20.6 Calibrating simulations within the potential outcomes framework
The potential outcomes framework provides a natural path for generating calibrated simulations for causal inference problems (See Section 15.1.6 for more discussion of calibration.)
Under the calibration framework, we would take an existing randomized experiment or observational study and then impute all the missing potential outcomes under some specific scheme. This fully defines the sample of interest and thus any target parameters, such as a measure of heterogeneity, are then fully known.
Data in hand, we then synthetically, and repeatedly, randomize and “observe” outcomes to be analyzed with the methods we are testing. We could also resample from our dataset to generate datasets of different size, or to have a superpopulation target as our estimand.
The key feature here is the imputation step: how do we build the full set of covariates and potential outcomes? There are a variety of options that have different pros and cons.
Regardless of how we generate the missing potential outcomes, once we have a “fully observed” sample with the full set of treatment and control potential outcomes for all of our units, we can calculate any target estimands we like on our population, and then compare our estimators to these ground truths (even if they have no parametric analog) as desired, just as we saw above.
20.6.1 Matching-based imputation
One approach for generating pairs of potential outcomes is to generate a matched-pairs dataset by, for each unit, finding a close match in terms of the covariates. The opposite pair of each matched unit then gives the imputed potential outcome.
By doing this (with replacement) for all units we can generate a fully imputed dataset which we can then use as our population, with all outcomes being “real,” as they are taken from actual data. Such matching can mostly preserve complex relationships between covariates and outcomes in a manner that is not model dependent. For example, if an outcome is coarsely defined (e.g., on an integer scale) or has specific clumps (such as from zero-inflation or rounding), this structure will be preserved.
One concern with the matching approach is the inexactness of the matching could dilute the structure of the treatment effect as the control- and treatment-side potential outcomes may be very unrelated, creating a lot of so-called idiosyncratic treatment variation (unit-to-unit variation in the treatment effects that is not explained by the covariates). This is akin to measurement error diluting found relationships in linear models.
We could try to reduce such variation by first imputing missing outcomes using some model (e.g., a random forest) fit to the original data, and then matching units including the imputed potential outcomes as additional “covariates.” This is not a data analysis strategy, but instead a method of generating synthetic data that both has a given structure of interest and also remains faithful to the idiosyncrasies of an actual dataset.
20.6.2 Model-based imputation
A second approach for imputing the missing potential outcomes is to specify a treatment effect model, predict treatment effects for all units and use the predicted effects to impute the treatment potential outcome for all control units and control potential outcomes for all treated units. For example, if we fit \(\hat{\tau}(X)\) to the original data, where \(\hat{tau}(X)\) estimates the treatment effect for a unit with covariate value \(X\), we can then impute \(Y_i(1) = Y_i(0) + \hat{\tau}(X_i)\) for all control units and \(Y_i(0) = Y_i(1) - \hat{f}(X_i)\) for all treated units. This imputation perfectly preserves the complex structure between the covariates and the \(Y_i(0)\)s for the original controls and covariates and \(Y_i(1)\) for the original treated.
One drawback is the imputed outcomes may not be values we actually see on the counterfactual side. For example, if the original outcomes are binary, but if \(\hat{\tau}(X)\) estimates fractional effects, then the imputed outcomes will not be binary, and could even be negative or above 1.
A second (probably minor) drawback is this form of imputation does not allow for any idiosyncratic treatment variation. To add in idiosyncratic variation we would need to generate a distribution of perturbations and add these to the imputed outcomes just as we would an error term in a regression model.
One advantage of the modeling approach is it allows for varying the amount of treatment effect. We can, for example, scale \(\hat{\tau}(X)\) by some factor to increase or decrease the size of the treatment effects in our imputed dataset.
20.7 Concluding Thoughts
Potential outcomes are a way of thinking about causal inference where we fully specify characteristics of the data, and then define estimands as summary measures across the full data (including aspects we would not normally see in the real world). This idea is a general one: in any simulation context, it is important to remember that the model for generating the data need not have any relation to how we estimate or analyze the data. By keeping these pieces separate in a simulation, it is much easier to generate data that is relevant to a specific context (e.g., by calibrating to an existing dataset) than when trying to use the same model used for estimation.
Keeping these parts separate has another important quality of helpling keep the researcher honest: it is usually the case that an estimator will work best on data generated from the estimator’s own working model. By generating using a different process, then the data does not inherently align with the estimator, meaning the estimator performance is less likely to be entirely a special case.