Treatment Effects¶
Average Treatment Effects¶
In the simplest setting, we are interested in estimating the causal effect of a binary variable \(W\) on an outcome variable \(Y\). Specifically, we are interested in the average treatment effect
The mathematical operator \(\mathrm{do}(W = w)\) denotes an intervention that holds the value of \(W\) constant at level \(w\) (see Pearl 2009 for more details). The average treatment effect can also be defined using the Neyman-Rubin potential outcomes framework (Imbens and Rubin 2015). WhyNot is agnostic to which framework the practitioner applies.
WhyNot allows the user to generate observational datasets \((X_i, W_i, Y_i)_{i=1}^n\) consisting of \(n\) samples of covariates \(X_i\), treatment assignment \(W_i\), and observed outcome \(Y_i\). Importantly, WhyNot also also returns the ground truth (sample) average treatment effect.
Case-study: Opioid Epidemic Simulation
The Opioid Epidemic Simulator is a system dynamics model of the US opioid epidemic. In the model, we ask:
What is the effect of lowering non-medical prescription opioid use by 10% in 2015 on the number of opioid overdose deaths in the United States in 2025?
Each unit is a run of the simulator from a different initial state. Treatment assignment \(W\) corresponds to whether the run receives the policy intervention in 2015, and the outcome \(Y\) is the number of overdose deaths in 2025.
To generate confounding, we imagine governments are more likely to intervene to reduce opioid abuse if the number of opioid overdose deaths is high. Therefore, runs with high levels of opioid overdose deaths in 2015 are more likely to receive treatment. A set of covariates \(X\) sufficient to estimate the treatment effect is the entire system state in 2015.
The OverdoseConfounding experiment implements this logic. The code below
calls the experiment class and generates the causal inference
Dataset for \(n=500\) samples.
>>> import whynot as wn
>>> import numpy as np
>>> overdose_experiment = wn.opioid.Confounding
>>> dataset = overdose_experiment.run(num_samples=500)
# True average treatment over the sample
>>> sample_ate = dataset.sate
>>> print(f"{sample_ate:.2f")
-16601.53
>>> X, W, Y = dataset.covariates, dataset.treatments, dataset.outcomes
# Replace with your favorite causal estimator
# Confounding significantly biases unadjusted estimates!
>>> estimate = np.mean(Y[W == 1.0]) - np.mean(Y[W == 0.0])
>>> print(f"{estimate:.2f}")
101444.41
One key parameter in the Confounding experiment is the strength of
the confounding. If runs with high levels of treatment get treated with
probability \(p\) and otherwise get treated with probability \(1-p\),
then how does the performance of finite-sample estimators change as \(p \to
1\)? We can easily generate a sequence of datasets to as \(p\) varies to
check this since \(p\) is a Parameter of the experiment.
>>> import whynot as wn
>>> overdose_experiment = wn.opioid.Confounding
>>> overdose_experiment.get_parameters()
Params:
Name: nonmedical_incidence_delta
Description: Percent decrease in new nonmedical users of prescription opioids.
Default: -0.1
Sample Values: [-0.073, -0.1]
Name: propensity
Description: Probability of treatment assignment in high overdose group.
Default: 0.9
Sample Values: [0.5, 0.6, 0.7, 0.9, 0.99]
>>> datasets = []
>>> for p in np.arange(0.5, 0.99, 0.05):
... dataset = overdose_experiment.run(num_samples=500, propensity=p)
... datasets.append(dataset)
Estimating Average Treatment Effects¶
With a causal dataset in hand, WhyNot provides a collection of estimators in the
causal_suite() to estimate average treatment effects. For a
detailed list of estimators, see Average Treatment Effect Estimators. Each estimator returns
an InferenceResult that includes the estimate of the
ATE, as well as a confidence interval (if provided by the estimator).
>>> import whynot as wn
>>> import numpy as np
>>> overdose_experiment = wn.opioid.Confounding
>>> data = overdose_experiment.run(num_samples=100)
# Replace with your favorite causal estimator
# Estimate ATE using a linear model
>>> estimate = wn.algorithms.ols.estimate_treatment_effect(data.covariates, data.treatments, data.outcomes)
# Compare estimate with ground truth
>>> relative_error = np.abs((estimate.ate - data.sate) / data.sate)
>>> print(f"{relative_error:.2f}")
0.01
Heterogeneous Treatment Effects¶
While average treatment effects are concerned with the causal effect over an entire population, heterogeneous treatment effects are concerned with the treatment effect for each individual or for each group defined by covariates \(X=x\). In particular, the Conditional Average Treatment Effect (CATE) for covariates \(x\) is defined as
Given an observational dataset \((X_i, A_i, Y_i)_{i=1}^n\), it is a challenging problem to estimate heterogeneous effects. WhyNot allows benchmarking of individual treatment effect estimations by returning individual level counterfactuals, i.e. both \(Y_{i, \mathrm{do}(W=0)}\) and \(Y_{i, \mathrm{do}(W=1)}\) for each sample \(i\).
Case-study: Opioid Epidemic Simulator To illustrate this, we consider the same study using the opioid epidemic simulator presented in the previous section.
>>> import whynot as wn
>>> overdose_experiment = wn.opioid.Confounding
>>> dataset = overdose_experiment.run(num_samples=200)
# True effects is a n x 1 vector of individual
# level contrasts Y_i(1) - Y_i(0)
>>> dataset.true_effects
# array([-16436.31184686, -16448.84326423, -16063.21459659, -16671.28477321,
# -16393.9382686 , -16533.25323364, ...])
Estimating Heterogeneous Treatment Effects¶
WhyNot provides a collection of estimators in the causal_suite() to
estimate heterogeneous treatment effects. See Heterogeneous Treatment Effect Estimators for a
detailed list. Each estimator returns an InferenceResult with
the property individual_effects. The code below shows how to use the
causal forest estimator to estimate
individual treatment effects for the OverdoseConfounding experiment in the
previous section.
>>> import whynot_estimators
>>> experiment = wn.opioid.Confounding
>>> dataset = experiment.run(num_samples=100)
# Estimate CATE using a causal forest
>>> estimate = whynot_estimators.causal_forest(
dataset.covariates, dataset.treatment, dataset.outcome)
# Compute MSE for HTE estimates
>>> mse = np.mean((estimate.individual_effects - dataset.true_effects) ** 2)