Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
46,468 characters · 10 sections · 36 citation commands
Covariate Balancing Sensitivity Analysis for Extrapolating Randomized Trials across Locations
Randomized control trials (RCTs) have long been the gold standard in evaluating public policy. In order to leverage insights learned from RCTs conducted in one location to drive policy change in other locations, there remains the concern about external validity, namely that the target population affected by the considered policy change may respond to the intervention differently than the population sample included in the RCT study cook2002experimental, olsen2013external. For example, perhaps the randomized study was run in Texas, but we want to use results from the study to inform policy across all 50 states. Or perhaps participation in the randomized study was voluntary. As a concrete example, attanasio2011subsidizing examine the effect of vocational training on labor market outcomes using a randomized study—but participants needed to apply to be part of the study, and so we may be concerned that the effect of vocational training for study participants may differ from the effect of training on the general population.
In cases where any cross-location divergences in populations can be explained away using observed covariates, the generalization problem is addressed by several authors, including dehejia2019local, hirshberg2019minimax, hotz2005predicting, stuart2011use, using ideas that generalize propensity score methods going back to rosenbaum1983central. In this paper, however, we are most interested in the case where locations differ along unobserved features, and so controlling for observed covariates doesn’t solve the cross-location generalization problem.
As discussed further below, our proposal builds upon the literature on sensitivity analysis in observational studies rosenbaum2002design, and in particular on recently proposed methods for sensitivity analysis via linear programming aronow2012interval, miratrix2017shape, zhao2019sensitivity. This literature is focused on cases where a randomized or observational study is confounded by unobserved features, and so the resulting treatment effect estimates do not have internal validity. Here, in contrast, we are not worried about the internal validity of our study; rather, we are concerned about how unobserved features may affect external validity when we aim to extrapolate results from RCT studies to other locations.
Our main finding is that in the cross-study generalization problem we have access to richer information than in the single-study sensitivity analysis problem, and that we can use this information to substantially improve power. At a high level, our approach is driven by the insight that we can collect data on observable features in both the location where the RCT was conducted and the target location we wish to extrapolate to, and that any selection weights that transport treatment effect estimates from one location to the other must respect moments of these observed covariates. As such, our motivation is qualitatively related to recent work on covariate-balancing estimation in causal inference athey2018approximate, graham2016efficient, hainmueller2012entropy, imai2014covariate, zubizarreta2015stable.
Generalizing results from randomized trials run in one location to other locations has been a central question in economic and healthcare policymaking cole2010generalizing, gechter2015generalizing, hotz2005predicting. Several authors have discussed generalizability bias (or transportability bias) that arises when there is either selection bias in a RCT or we try to generalize the results of one RCT to a different location that has a different covariate distribution bareinboim2016causal, olsen2013external, pressler2013use. To address covariate shift in the case of no measured confounding, many existing works use the generalized propensity scores for either matching or weighting cole2010generalizing, dehejia2019local, hirshberg2017augmented,stuart2011use, westreich2017transportability, or covariate matching gretton2009covariate. We also note the recent advances in the literature on transportability in causal graphical models that focuses on designing graphical algorithms for identifying whether transportability is feasible from graphs bareinboim2013meta, bareinboim2014transportability, bareinboim2016causal, pearl2014external. Our work is complementary to theirs in the sense that we don't assume transportability is necessarily feasible due to unmeasured effect modifiers, and instead we compute identification intervals whose length varies with the distributional imbalance of these unobserved factors across locations.
We build upon the literature on sensitivity analysis in observational studies cornfield1959smoking, fogarty2016sensitivity, imbens2003sensitivity, rosenbaum1987sensitivity, rosenbaum2002attributing, rosenbaum2002design, rosenbaum2010design, rosenbaum2011new, rosenbaum2014weighted, shen2011sensitivity, vanderweele2017sensitivity and in particular, the recently developed works on using linear programming aronow2012interval, miratrix2017shape,yadlowsky2018bounds, zhao2019sensitivity. This line of work focuses on using sensitivity analysis to address internal validty of an observational study. On the other hand, our work is focused on addressing external validty of an RCT. The sensitivity model we have employed is closest in connection with the marginal sensitivity model advocated in tan2006distributional and zhao2019sensitivity.
Given that violations of transport assumptions is untestable for generalizing results from RCT to other locations, a few authors have proposed methods that perform a sensitivity analysis on how much violation on this assumption can lead to generalizability bias. andrews2017weighting considers the limit of the role of private information, and their approach only applies to the special case where the randomized study population is a subset of the target population we wish to apply a policy intervention to. gechter2015generalizing considers the range of distributions of the treated potential outcomes conditioning on the controls potential outcomes. nguyen2017sensitivity considers the ranges on the mean of the unobservabed covariates. We view these approach as complementary to our proposal, as we consider the limit on the ratio of the probability density function on the unobserved covariates among locations, and further exploit covariate balance for improved power.
The core of our proposal is that by leveraging covariate information in the locations of interest, the set of weights for weighting outcomes can be used to balance moments on the covariates between locations. This would shorten identification intervals for the target population compared to prior work. In particular, the idea of covariate balancing has become popular in other contexts such as estimating average treatment effect athey2018approximate, bennett2018building, graham2016efficient,hainmueller2012entropy, hirshberg2017augmented, kallus2018balanced, wang2017minimal, zhao2019covariate,zubizarreta2015stable.
We note that there is a growing interest in combining observational data with RCTs athey2016estimating, kaizar2011estimating, kallus2018removing, rosenman2018propensity to improve power. In our work, we focus on only using the RCT data, and leave it to future work on how to incorporate observational data in our framework. We also note the interesting and promising direction of using Bayesian hierarchical models for combining RCT results from multiple studies meager2016aggregating,vivalt2016much. They work in the setting of using aggregated experimental data on the study level, whereas our work leverages individual-level data. There is also a growing literature on transfer learning for domain shifts (see, e.g., bastani2020predicting and references therein for a review); in our case, we focus on transfer results from a randomized trial from one location to another.
We assume there are two locations (e.g., two cities). In the first location, an RCT has been conducted to evaluate the impact of an intervention on some outcome of interest, e.g., whether a job training program improved participants' earnings. In the second location, without running an additional RCT tailored for this new population, we want to ask the question of to what extend the previous RCT result is applicable to the new location. In both locations, we observe a set of covariates such as each citizen's age and marital status, but we might not have access to other important covariates such as their education level.
Formally, we denote the two locations by Location $0$ and Location $1$. Location $0$ is where the RCT is conducted, and we wish to generalize the RCT results to Location $1$. Throughout the paper, we use the subscripts of $0$ and $1$ to denote respective quantities for the two locations. We assume the data are i.i.d generated in each location: $X_i, U_i, W_i, Y_i(0), Y_i(1) \overset{\text{i.i.d}}{\sim} P_k$ for $i=1, \cdots, n_k$ where the subscript $k\in \{0,1\}$ corresponds to the two locations respectively, $X \in \mathcal{X}$ is the observed covariates, $U \in \mathcal{U}$ is the unmeasured covariates, $W$ is the randomized treatment assignment indicator where $W=1$ indicates that the treatment is assigned, and $W=0$ indicates otherwise, and $Y(w)$ for $w=0,1$ is the potential outcome corresponding to having received treatment or lack thereof. In Location $0$, we observe data $X, Y, W$ with $Y=Y(W)$, and the propensity score $e:=P_0(W=1)$ is a known constant. On the other hand, in Location $1$, we assume we only observe the covariates $X$. The causal estimand is $\tau = \EE[1]{Y(1)-Y(0)}$, denoting the average treatment effects in the target Location $1$. For convenience, we use $L_i \in \{0,1\}$ to denote the location indicator for each unit.
We start with a few assumptions.
This assumption implies that the distribution of the potential outcomes conditioning on the observed variables $X$ and any unmeasured effect modifiers $U$ is the same across the two locations. Most existing works in the literature assume transportability only conditioning on the observed covariates hotz2005predicting. Given that the unmeasured effect modifier $U$ can have any association with the potential outcomes $Y(0)$ and $Y(1)$ conditionally on $X$, we note that Assumption (ref) does not impose any meaningful restrictions on its own.
The above proposition shows that if the unmeasured effect modifiers were known and measured in both locations, then the causal quantity of interest $\tau$ is identified. In order to turn the above into an estimator, we could use a standard inverse weighted estimator:
such that $\EE[0]{\hat{\tau}^*} = \tau$. Since $U$ is unobserved, the above estimator is not feasible. Instead, we further allow sensitivity models detailed in the next subsections which assume bounds on the ratio of the conditional probability density of unmeansured effect modifiers $U$. This would in turn allow us to derive bounds for $\hat{\tau}$ via a linear programming optimization.
To relate $r(x,u)$ and $r(x)$, we define the unobserved distribution shift ratio as
and use the shorthand $z_i^*:= z^*(X_i, U_i)$. These $z_i^*$ capture the amount by which the unobserved effect modifiers $U_i$ affect the oracle estimator (ref), which we can re-write as
Given we don't know the density ratio terms $r(X, U)$ in (ref) due to the unmeasured effect modifier $U$, we instead aim to get bounds on $\tau$ by estimating $r(X)$ which does not depend on $U$ and by assuming a bound on $z^*(x,u)$. In particular, we assume the following sensitivity model that directly implies a bound on $z(x,u)$.
By Bayes rule, we immediately have the following:
Analogous sensitivity models are common in the sensitivty analysis literature for observational studies rosenbaum2002design, zhao2019sensitivity to assess robustness of findings to unmeasured confounding.
With a bound on $z^*(x,u)$ and an estimated density ratio $\hat{r}(\cdot)$, we can derive bounds on $\hat{\tau}$ from (ref). Since we don't know the ground truth value of $\Gamma^*$, we supply a $\Gamma$ value and solve the following optimization problem to get the upper bound $\tilde{\tau}^+$.
To get the lower bound $\tilde{\tau}^-$, we instead take the infimium in the optimization above. This key idea of employing a linear program to bound the target estimand when the density ratio is unknown but can vary within some range is also seen in previous literature. For example, aronow2012interval and miratrix2017shape employ similar ideas for identification of a population mean when the sampling selection weight is unknown.
The roadmap for the rest of this section is as follows: We first show how to estimate $\hat{r}$ in a way that allows us to achieve sharper bounds. Next, we use the estimate $\hat{r}$ to get the upper bound by (ref) (and respectively, the lower bound by taking the infimum of the same optization) for $\tilde{\tau}$. Finally, we conclude with showing the bounds $[\tilde{\tau}^-, \tilde{\tau}^+]$ gives consistent coverage for the ground truth treatment effect $\tau$ in Location 1.
First, we estimate $\hat{r}(\cdot)$. While (ref) takes in $\hat{r}(\cdot)$ from any density ratio estimators, we suggest estimating $r(x)$ directly by moments matching bradic2019sparsity, imai2014covariate, ning2017high, qin1998inferences, sugiyama2012density. Doing so would allow us to effectively shorten the estimation bounds as detailed later in Section (ref). To fix ideas, we denote $g(x) = \log r(x)$ and then assume the following:
When the basis function $\phi(\cdot)$ is the identity, this simply implies that the density ratio follows a logistic form. On the other hand, we can take $\phi(\cdot)$ to be any bounded basis expansion, which makes the above assumption not as restrictive. In Section (ref), we further relax this assumption and incorporate the model misspecification error. For the rest of the paper, we assume $\phi(X_i)$ contains an intercept term.
By Assumption (ref) and the fact that an RCT is conducted in Location $0$, we have the following covariates balancing moment condition via a change of measure: for $w=0,1$,
We exploit the empirical counterpart of the above moment equation to estimate the $l(x)$. Given Assumption (ref), we only need to estimate $\beta^{(w)}$. Concretely, we solve the following optimization problem:
and then set
The reason we fit \smash{$\hat{\beta}_\phi^{(w)}$} separately on the treated and control units is that it enables us to exactly match the empirical version of the moment condition in (ref) for both treated and control groups (this follows immediately form the first-order condition in (ref)):
We note that, if \smash{$\beta_\phi$} is unique, then both \smash{$\hat{\beta}_\phi^{(w)}$} will eventually converge to \smash{$\beta_\phi$} in large samples.
By Assumption (ref), the density ratio is identified, and the estimated density ratio is consistent, i.e. $\Norm{\hat{g}_\phi(x) - g(x)}_\infty \to_p 0$ qin1998inferences, sugiyama2012density. There is an exact mapping between estimating density ratios and estimating location propensity by treatment locations as random variables. We use a covariate balancing estimator for the density ratio, which corresponds to the covariate balancing estimator for propensity scores if we treated the locations as random variables. In the line of work for propensity score estimation, our covariate balancing approach is closely related to imai2014covariate, tan2017regularized, tan2018model,zhao2019covariate. We can then substitute in $\exp\p{\hat{g}_\phi(\cdot)}$ for $\hat{r}(\cdot)$ in (ref) to compute the corresponding upper and lower bounds $\tilde{\tau}^+$ and $\tilde{\tau}^-$.
Next, we show the resulting bounds give consistent coverage:
Although the optimization formulation in (ref) gives consistent coverage on $\tau$, the resulting bounds can often be too wide in practice to provide meaningful conclusions about the treatment effect $\tau$. As a concrete example, assume the age distribution is the same for the location that we have conducted the RCT in and the location we wish to extrapolate findings to. We may find that from the RCT, the treatment effect is largest among the young population. The optimization procedure detailed in the previous section does not leverage information on age distributions in these two locations, and could construct weights $z_i$ that overweight or underweight the young population leading to wide bounds.
The key contribution of this work is to sharpen the estimation bounds by further leveraging covariate information available in both locations. In particular, similar to (ref), we can also balance the following moments using a change of measure that conditions on both the observed $X$ and the unmeasured effect modifiers $U$. Let $g(x,u) = \log r(x,u)$. Then the following moment equation holds: for $w=0,1$,
We then add the empirical counterpart of the moment equality constraints in (ref) to (ref) as the last two equalities (a) and (b) below.
Similarly, we can estimate $\hat{\tau}^-$ by taking infimum of the above optimization with the same set of constraints. Compared to (ref), (ref) makes a few modifications. First, it substitutes in the estimate $\exp\p{\hat{g}_{\phi}(\cdot)}$ for $\hat{r}(\cdot)$ using (ref). Second, it adds two equality constraints (a) and (b). By setting the derivative of (a) and (b) to 0 in expectation, we arrive at the moment condition in (ref) for $w=0,1$ respectively.
Covariate balancing via (ref) shortens estimation intervals compared to solving (ref) due to the added equality constraints which limit the plausible range of $z_i$'s. To provide more intuition, we include in the Appendix a more technical comparison in the simple context where the covariates $X$ are discrete and the link function $\phi$ is the identity function.
The constructed intervals from the above optimization (ref) gives consistent coverage of the underlying treatment effect parameter of interest.
This implies that as the sample size increases, if we employ a $\Gamma$ value in the optimization that is no less than what is needed for Assumption (ref) to hold, then the constructed interval from the balancing estimator eventually covers the true treatment effect in the target location that we wish to apply policy interventions to.
So far, we have assumed that the density ratio follows a logistic model with a basis expansion $\phi(\cdot)$ by Assumption (ref). In this section, we relax this assumption and develop bounds that take model misspecification error into account.
Given any basis function $\phi(\cdot)$, let $\beta_\phi^{(w)}$ be the population minimizer, i.e. for $w=0,1$,
and let
be the logistic approximation to the ground truth density ratio logit function $g(\cdot)$. For the rest of this section, instead of Assumption (ref), we build upon the following sensitivity model instead.
Then immediately by Bayes rule, we have
In practice, we don't know the magnitude of $e^{g(x)} / e^{g_\phi(x)}$. We supply an additional sensitivity parameter $M$ such that $M \geq e^{g(\cdot)} / e^{g_\phi(\cdot)}$ based on how much we believe the density ratio is misspecified with a chosen logistic model. We proceed to get bounds on $\tau$ by estimating $\hat{g}_\phi(\cdot)$ via (ref). We solve for $\hat{\tau}_\phi^+$ by (ref) but we substitute in $\Gamma M$ for $\Gamma$, and similarly we take the infimum to derive $\hat{\tau}_\phi^-$. Compared to the previous sections, we relax the constraint that $g(x)$ is well specified with a logistic form, and we conclude with the same consistency result if we assume Assumption (ref) instead of Assumption (ref).
We apply our proposed estimator to the California Greater Avenue for Independence (GAIN) dataset hotz2006evaluating. A policy analyst may run an RCT in one location and wishes to know the treatment effect in another location without running additional RCTs. By using our proposal, they could get bounds on the estimated treatment effects in the second location directly. We validate this proposal on the GAIN dataset which includes data from independent RCTs conducted in several selected counties in California to evaluate the impact of welfare-to-work programs on an individual's future income. We focus on two counties: Los Angeles and Riverside. Suppose we had only run the RCT in Los Angeles and would like to extrapolate the results to Riverside. If there are no unmeasured effect modifiers that would affect both the outcome and location likelihood, we can simply use a weighted Hajek-style estimator. The goal is to get an uncertainty quantification of the extrapolated results in Riverside with varying degrees of how strong the unmeasured covariate shift is. Since the GAIN dataset includes the RCT data from Riverside, it allows us to validate the estimated bounds from our proposal against the ground truth treatment effect estimates from the RCT that had been conducted in Riverside.
We run the proposed estimator with covariate balancing to generalize treatment effect bounds from Los Angeles to Riverside and vice versa (referred to as “covariate balance" in Figure (ref)). For comparison, we also run the proposed estimator without leveraging covariate information (referred to as “no balance" in Figure (ref)). We use the mean quarter income over a follow-up period of 9 years post experiment as the outcome. For both counties, we use a simple difference-in-means between the the treated and the control groups to compute their ground truth treatment estimates (referred to as “ground truth" in the legend of Figure (ref)). Let $\gamma := \log \Gamma$. We quantify the strength of the unmeasured covariate shift by assuming different $\Gamma$ values in (ref) and (ref) for the “no balance" and "covariate balance" approaches respectives.
The left plot in Figure (ref) shows generalizating the RCT estimate from Los Angeles to Riverside and the right plot shows the generalization results the other way. By varying $\gamma$ along the x-axis, we vary the assumed bound on the distributional shift of the unmeasured effect modifiers. We see that the “covariate balance" approach meaningfully shortens the estimated bounds compared to the “no balance" approach, and the bounds give coverage to the ground truth estimates as $\gamma$ increases. To estimate the variance, we generate 1000 bootstrap samples in both locations simultaneously to account for stochastic fluctuations in the data. The dashed lines show the 95% confidence intervals using the percentile bootstrap, following zhao2019sensitivity.
We consider the following setups adapted from the simulation study in yadlowsky2018bounds. For some covariate distribution $P$,
where $L=0,1$ denotes the location indicator. We consider the following two setups:
For Setup A, we let $X_k \sim Beta(0.5, 0.5)$ for $k=1,2,3,4$ and $\mu = [2,2,-2,-2]$ such that $P(X\cond L)$ highly depends on the location $L$. We vary $\gamma \in [0,\,0.1,\,0.2, \,0,3,\,0.4,\,0.5]$ and let $\sigma = 3$ and $\alpha_0=0$. We let $\Gamma^* = exp(0.2)$ be the true sensitivity parameter, but we assume it's unknown to us. For Setup B, we let $X \sim Uniform[0,1]^4$, $\Gamma^* = \exp(0.5)$, $\alpha_0 = -2, \sigma = 0.5$, and vary the sensitivity parameter $\gamma = \log{\Gamma} = 0, \,0.1, \,\ldots,\,0.7$ in the optimization procedure, We let $\mu = [0.709, 0.438, 0.2, 0.767]$ in Setup B and let $\beta = [0.513, 0.045, 0.7, 0.646]$ in both setups.\footnote{Both parameters are taken from the simulation in yadlowsky2018bounds}
The form of $U$ is chosen such that Assumption (ref) holds with $\Gamma^*$, with $l(x) = \exp(\alpha_0 + x^\top \mu )/{\{1+\exp(\alpha_0 + x^\top \mu )\}}$, and $P(U=u\cond X=x, L=1) / P(U=u\cond X=x, L=0) = {\Gamma^*}^{\mathbbm{1}_{u>0} - \mathbbm{1}_{u < 0}}$. and we let the size of the total combined population across the two locations to be 1000.\footnote{We note that for the purpose of this simulation setup, it is natural to define location as an additional random variable in the data generating process to ensure the ground truth sensitivity bound to fall within $[1/\Gamma^*, \Gamma^*]$.} We use the mosek package for optimization, and we compare the percentile bootstrap confidence interval obtained through 1000 bootstrap samples among the difference-in-means estimator in location $L=1$ and the proposed estimator with covariate balancing. We see that the covariate balancing approach signficantly shortens the estimation interval, while Setup A in Figure (ref) also shows that our balancing estimator (in red) is not conservative as its confidence interval just covers the ground truth (in green) once the $\gamma$ parameter is increased to $\gamma^*=0.2$ in this case.