Chapter 23 Simulation under the Potential Outcomes Framework
If we are in the business of evaluating how various methods such as matching or propensity score weighting work in practice, we would probably turn to the potential outcomes framework for our simulations. The potential outcomes framework is a framework typically used in the causal inference literature to make very explicit statements regarding the mechanics of causality and the associated estimands one might target when estimating causal effects. While we recommend reading, for a more thourough overview, either [CITE Raudenbush or Field Experiments textbook], we briefly outline this framework here to set out our notation.
Take a sample of experimental units, indexed by \(i\). For each unit, we can treat it or not. Denote treatment as \(Z_i = 1\) for treated or \(Z_i = 0\) for not treated. Now we imagine each unit has two potential outcomes being the outcome we would see if we treated it (\(Y_i(1)\)) or if we did not (\(Y_i(0)\)). Finally, our observed outcome is then \[ Y_i^{obs} = Z_i Y_i(1) + (1-Z_i)Y_i(0) .\] For a unit, the treatment effect is \(\tau_i = Y_i(1) - Y_i(0)\); it is how much our outcome changes if we treat vs. not treat. Frustratingly, for each unit we can only see one of its two potential outcomes, so we can never get an estimate of these 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 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 proportion increase, by contrast, would be
\[ API_{\mathcal{S}} = \frac{1}{N} \sum_{i} \frac{Y_i(1)}{Y_i(0)} \]
23.1 Finite vs. Superpopulation inference
Consider a sample of \(n\) units, \(\mathcal{S}\), along with their set of potential outcomes. We can talk about the true ATE of the sample, or, if we thought of the sample as being drawn from some larger population, 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, our estimand could be (maybe only slightly) 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 for each sample generated.
The catch is when we calculate our performance metrics, we now have two possible targets to pick from. Furthermore, if we are targeting the superpopulation 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.
23.2 Data generation processes for potential outcomes
If we want to write a simulation using the potential outcomes framework, it is clear and transparent to first generate a complete set of potential outcomes, then generate a random assignment based on some assignment mechanism, and finally generate the observed outcomes as a function of assignment and original potential outcomes.
For example, we might say that our data generation process is as follows: First 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 is saying that 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\), assuming \(\beta_1\) and \(\tau_1\) are non-zero.
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 to illustrate the first part of the data generating process (we leave treatment assignment to later):
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 )
}And now we see our estimand can change:
## [1] 0.6374925## [1] 0.5479788In reviewing our 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.
Also in looking at the above, there are a few details to call out:
- We can store the latent, intermediate 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.
- We are using a trick to index our DGP by an R2 value rather than coefficients 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 heterogeniety of the treatment impacts.
- If we were generating data with a constant treatment impact, then \(ATE_{\mathcal{S}} = ATE\) always; this is typical for many similations in the literature. That being said, treatment variation is what causes a lot of methods to fail and so having simulations that have this variation is usually important.
Once we have our schedule of potential outcomes, we would then 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 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:
## # 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) is randomly assigning these whole units to various conditions and observing the consequences.
In general, we might instead 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).
23.3 Finite sample performance measures
Let’s 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:
## [1] 1.154018Our finite-population simulation would be:
rps <- rerun( 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"]] )
}) %>%
bind_rows()## Warning: `rerun()` was deprecated in purrr 1.0.0.
## ℹ Please use `map()` instead.
## # Previously
## rerun(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"]])
## })
##
## # Now
## map(1: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"]])
## })
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.## # A tibble: 1 × 3
## EATE_hat SE ESE_hat
## <dbl> <dbl> <dbl>
## 1 1.16 0.248 0.307We 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. Now this opens up some room for critique: what if our single dataset is non-standard?
Our super-population simulation would be, by contrast:
rps_sup <- rerun( 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"]] )
}) %>%
bind_rows()## Warning: `rerun()` was deprecated in purrr 1.0.0.
## ℹ Please use `map()` instead.
## # Previously
## rerun(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"]])
## })
##
## # Now
## map(1: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"]])
## })
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.## # A tibble: 1 × 3
## EATE_hat SE ESE_hat
## <dbl> <dbl> <dbl>
## 1 1.00 0.381 0.378First, 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 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.)
23.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 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 phenomonen 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. This is what we did in [CITE estimands paper]
Borrowing from the simulation appendix of [CITE estimands paper], 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. Each dataset includes, for each observation, the potential outcome under treatment or control.
Record the true finite-sample ATE, both person and site weighted.
Then, three times, do a finite simulation as follows:
- Randomize units to treatment and control.
- 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 variability for each dataset, but the average across the \(R\) datasets in a given scenario gives a reasonable estimate of expected variability 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 <- rerun( 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:
## Warning: `rerun()` was deprecated in purrr 1.0.0.
## ℹ Please use `map()` instead.
## # Previously
## rerun(3, {
## 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"]])
## })
##
## # Now
## map(1:3, ~ {
## 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"]])
## })
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.## # 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.
## Warning: `rerun()` was deprecated in purrr 1.0.0.
## ℹ Please use `map()` instead.
## # Previously
## rerun(500, one_finite_run())
##
## # Now
## map(1:500, ~ one_finite_run())
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.We use .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 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 observations with every estimate.
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.
We need to keep the squaring. If we look at the SEs themselves, we have further apparent bias due to our estimated ESE_hat being so unstable due to too few observations:
## [1] 0.2944556We can 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):
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.9