Learning Causal Effects From Many Randomized Experiments Using Regularized Instrumental Variables

Randomized experiments (i.e. A/B tests, randomized controlled trials) are a popular practice in medicine, business, and public policy [6, 25]. When decision-makers employ experimentation they have a far greater chance of learning true causal relationships and making good decisions than via observation alone [21, 26, 29]. However, a single experiment is often insufficient to learn about the causal mechanisms linking multiple variables. Learning such multivariate causal structures is important for both theory building and making decisions [17, 22].

Consider the situation of a internet service for watching videos. The firm is interested in how watching different types of videos (e.g., funny vs. serious, short vs. long) affects user behaviors (e.g., by increasing time spent on the site, inducing subscriptions, etc.). Such knowledge will inform decisions about content recommendation or content acquisition. Even though the firm can measure all relevant variables, training a model on observational data will likely be misleading. Existing content recommendation systems and heterogeneous user dispositions will produce strong correlations between exposure and time spent or subscription, but the magnitude of this correlation will, in general, not match what would occur if the decision-maker intervened and changed the promotion or availability of various video types. Thus, we are interested not just in prediction but prediction under intervention [9, 10, 31].

The standard solution is to run a randomized experiment exposing some users to more of some type of video. However, a single test will likely change many things in the complex system. It is hard to change the number of views of funny videos without affecting the number of views of serious videos or short videos. This is sometimes called the problem of ‘fat hand’ interventions because such interventions touch multiple causal variables at once and so the effect on a single variable is not identified. To solve this issue the company would need to experiment with several factors simultaneously, perhaps conducting new experiments specifically to measure effects via each mechanism [22].

However, because routine product experimentation is common in internet companies [5, 25, 39], this firm has likely already run many A/B tests, including on the video recommendation algorithm. The method proposed in this paper can either be applied to a new set of experiments run explicitly to learn a causal effect vector [as in, e.g., 13], or can be applied to repurpose already run tests by treating them as random perturbations injected into the system and using that randomness in a smart way.

Our contributions arise from adapting the econometric method of instrumental variables [IV; 41, 33, 1] to this setting. It is well known that a standard IV estimator — two-stage least squares (TSLS) — is biased in finite samples [3, 36]. For our case, it also has asymptotic bias. We show that this bias depends on the distribution of the treatment effects in the set of experiments under consideration.

Our main technical contribution is to introduce a multivariate l ₀ regularization into the first stage of the TSLS procedure and show that it can reduce the bias of estimated causal effects. Because in finite samples this regularization procedure reduces bias but adds variance, we introduce a method to trade these off and select a regularization parameter. We call this procedure instrumental variables cross-validation (IVCV). In an empirical evaluation that combines simulation and data from hundreds of real randomized experiments, we show that the l ₀ regularization with IVCV outperforms TSLS and a Bayesian random effects model.

Finally, we show how to perform this estimation in a computationally and practically efficient way. Our regularization and cross-validation procedures only require summary statistics at the level of experimental groups. This is advantageous when using raw data is computationally or practically burdensome, e.g., in the case of internet companies. This means the computational and data storage complexities of the method are actually quite low. In addition, standard A/B testing platforms [5, 42] should already compute and store all the required statistics, so the method here can be thought of as an “upcycling” of existing statistics.

Suppose we have some (potentially vector valued) random variable X and a scalar valued outcome variable Y. We want to ask: what happens to Y if we change some component of X by one unit, holding the rest constant? Formally, we study a linear structural (i.e. data generating) equation pair

In general, we are interested in estimating the causal effect β because we are interested in intervention, e.g., one which will change our data-generating model to

In the presence of unobserved confounders trying to learn causal relationships using predictive models naively can lead us astray [9, 10, 31, 34]. Suppose that we have observational data of the form (X, Y) with U completely unobserved. If we use this data to estimate the causal effect β we can, due to the influence of the unobserved confounder, get an estimate that is (even in infinite samples) larger, smaller or even the opposite sign of the true causal effect β.

To see this consider the linear structure equation model above and suppose that we only observe (X, Y) where both are scalar. Since the underlying model is linear, we can try to estimate it using a linear regression. However, not including the confounder U in the regression yields the estimator:

When all variables are scalar algebra yields

Thus, the best linear predictor of Y given X ($\hat{\beta }_\text{obs}$ ) may not be lead to a good estimate of what would happen to Y if we intervened (β).

Figure 1: DAG representing our structural equations, in which the relationship between X and Y is confounded by U, and including the instrumental variable Z. Crosses represent causal relationships that are ruled out by the IV assumptions.

We now discuss instrumental variable (IV) estimator as a method for learning the causal effects. Suppose that we have some variable Z that has two properties (see Figure 1 for the directed acyclic graph which represents these assumptions):

1. Z is not caused by anything in the (X, U, Y) system; that is, Z is as good as randomly assigned.
   
2. Z affects Y only via X. This is often called the exclusion restriction or complete mediation assumption [3]
   

In terms of the structural equations, this modifies the equation for X to be

The standard IV estimator for β is two-stage least squares (TSLS) and works off the principle that the variance in X can be broken down into two components. The first component is confounded with the true causal effect (i.e. comes from U). The second component, on the other hand, is independent of U. Thus, if we could regress Y only on the random component, we could recover the causal effect β. Knowing Z allows us to do nearly this (i.e. by using only the variation in X caused by Z not U).

TSLS can be thought of as follows: in the first stage we regress X on Z. We then replace X by the predicted values from the regression. In the second stage, we regress Y on these fitted values. It is straightforward to show that as n approaches infinity this estimator converges to the true causal effect β [40, Theorem 5.1].

All IV methods make a full rank assumption. In order to estimate the effect of each variable X_j on Y with the other X’s held constant it must be the case that Z is such that it causes independent variation in all dimensions of X. This implies that we must, at least, have as many instruments as the dimension of β for TSLS to work. An interesting and fruitful direction for future work is what to do when some subspace of X is well spanned by our instruments but some some subspace is not.

In our setting of interest, randomly assigned groups from a large collection of experiments are the instruments.

Formally, here the IV is a categorical variable indicating which of K test groups a unit (e.g., user) was assigned to in one of many experiments. For simplicity of notation, we assume that each treatment group g ∈ {1, ..., K} has exactly n_g = n _per units assigned to it at random.

We now introduce a regularization procedure that can decrease bias in the TSLS estimator. We show that, in this setting a l ₀-regularized first stage is computationally feasible and can help reduce this bias under some conditions on the distribution of the latent treatment effects.

4.1 Intuition via a Mixture Model

There are many types of A/B tests conducted — some are micro-optimizations at the margin and some are larger explorations of the action space. Consider the stylized case with two types of tests calling the smaller variance type ‘weak’ tests while the larger variance ones are ‘strong.’ Here we can model the first stage $\bar{Z}$ as being drawn from a two-component mixture distribution:

\[ \bar{Z} = {E}[X | Z] \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\mathcal {N} (0, \sigma ^2_{weak}) & \text{with probability } p\\ \mathcal {N} (0, \sigma ^2_{strong}) & \text{with probability } (1-p)\\ \end{array}\right. \]

If we knew which group was drawn from which component and ran two separate TSLS procedures using only groups whose $\bar{Z}$ is drawn from the same component, we would asymptotically get two estimators:

\[ \mathop{plim}_{K \rightarrow \infty } \hat{\beta }_{\text{TSLS}, j} = \beta + \gamma \frac{ \psi \frac{\sigma ^2_{U}}{{n_{\text{per}}}}}{\psi ^{2} \frac{\sigma ^2_U}{{n_{\text{per}}}} + \frac{\sigma ^2_{\epsilon _{X}}}{{n_{\text{per}}}} + \sigma ^2_{j}} \]

Here j ∈ {weak, strong} representing the component from which a particular group's $\bar{Z}$ is drawn. Because $\sigma ^2_\text{strong} {\gt} \sigma ^2_\text{weak}$ we will have that $\hat{\beta }_{\text{TSLS, strong}}$ is a less asymptotically biased (and thus asymptotically better) estimator than $\hat{\beta }_{\text{TSLS, weak}}.$ Thus, if we could choose, we would choose to only use strong tests for our estimation of the causal effect.

In reality, we would likely not know which component each group is drawn from and if simply ran a TSLS on the full data set, this estimator will be a weighted combination of the two estimators.

Figure 2: Comparison of stagewise vs. IVCV method. X-axis is the strength of regularization (lower p-value implies stronger regularization). Optimizing for stagewise loss would imply using almost no regularization whereas optimizing for IVCV loss implies strong regularization. Causal loss coincides much more with IVCV loss than stagewise loss.

Within this discrete mixture model, we are limited to how much we can reduce bias (since $\mathop{plim}_{K \rightarrow \infty } \hat{\beta }_\text{TSLS, strong} \ne \beta$ ). However if the treatment effects are drawn from a distribution which is an infinite mixture of normals that has full support on normals of all variances (for example a t distribution) then we can asymptotically reduce the bias below any ϵ by using only observations which come from components with arbitrarily large variances.

We now introduce a regularization procedure that tries to perform this selection. Because using this regularization effectively decreases our dataset size, decreasing the bias increases the variance. Thus, afterwards we will turn to a procedure to set the regularization parameter to make good finite-sample bias–variance tradeoffs.

We now turn to an important practical question: because there is a bias–variance tradeoff how should one set the regularization parameter when K is finite to optimize for prediction under intervention?

First, let us suppose that we have access to the raw data where a row is a (X_i , Z_i , Y_i ) which is a unit i’s, X, Y and treatment assignment Z. We propose a procedure to set our hyperparameter q. We describe 2-fold version as it conveys the full intuition, but extension to k-folds is straightforward.

Instrumental variables cross-validation algorithm (IVCV):

1. Split each treatment in the data set into 2 folds, call these new data sets $\lbrace (X^1_i, Y^1_i, Z^1_i) \rbrace$ and $\lbrace (X^2_i, Y^2_i, Z^2_i) \rbrace$ .
   
2. Compute treatment level averages $\lbrace (\bar{X}^1_g, \bar{Y}^1_g) \rbrace$ and $\lbrace (\bar{X}^2_g, \bar{Y}^2_g) \rbrace$ as described above where g now indexes experimental groups.
   
3. Compute $\hat{\beta }_{q}$ for a variety of thresholds q using $\lbrace (\bar{X}^1_g, \bar{Y}^1_g) \rbrace$ .
   
4. Compute treatment level predictions of Y using fold 1 for each level of q: $\hat{Y}^q_g = \bar{X}^1_g \hat{\beta }_q$ .
   
5. Choose q which minimizes $\text{IVCV}(q) = \sum _{j} (\bar{Y}^2_g - \hat{Y}^q_g)^2.$
   

The intuition behind IVCV is similar to the main idea behind IV in general. Recall that our objective is to use variation in X that is not caused by U. The IVCV algorithm uses the X value from fold 1 and compares the prediction to the Y value in fold 2 because fold 1 and fold 2 share a Z but differ in U (since U is independent across units but Z is the same within group). This intuition has been exploited in split-sample based estimators [2, 19, 23].

We can demonstrate the importance of using the full causal loss by comparing the IVCV procedure to other two candidates. The first is simply applying naive CV in the second stage (i.e., splitting each group into 2, training a model on fold 1 and computing the CV loss naively as $\Vert Y_2 - X_2 \hat{\beta }_q \Vert ^2$ ). The second is stagewise, in which the regularization parameter is chosen to minimize MSE in the first stage, and then the second stage is fit conditional on the selected model [as in 8, 20]. We compare these approaches in a simple linear model with scalar X, such that

Recall that our goal is to find a hyperparameter (regularization threshold) which gives us the best prediction $\hat{\beta }_q$ of the causal parameter. Formally, we write this as

This gives us 3 candidate cross-validation loss functions to compare to the true casual loss in our simulations

1. The second stage CV loss
   \[ \sum _{g} (\bar{Y}^2_g - \bar{X}^2_g \hat{\beta }_q)^2 \]
   
   
   
2. The first stage CV loss
   \[ \sum _{g} (\bar{X}^2_g - \hat{\bar{X}}^{2q}_g)^2 \]
   
   
   
3. The IVCV loss
   \[ \sum _{g} (\bar{Y}^2_g -\bar{X}^1_g \hat{\beta }_q)^2 \]
   
   
   

Figure 2 shows these losses as a function of the parameter q averaged over 500 simulations of the model above. We see that both the first stage loss curve and the second stage loss curve look very different from the causal loss curve. However, the IVCV loss curve matches almost exactly. Thus, either stage error naively yields a very different objective function from minimizing the causal error. In particular, we see that making the bias–variance tradeoffs for the first stage need not coincide with a desirable bias-variance tradeoff for inferring β.

Figure 3: Left: Two dimensions of the multivariate means for sampled test groups ($\bar{X}_g$ ). The joint distribution exhibits many more outliers than a multivariate normal would suggest. Right: QQ-plots for the dimensions of the sampled test groups ($\bar{X}_g$ ). For each variable the marginal distributions are notably non-normal.

The l ₀-regularized IV estimator only requires the kinds of summary statistics per experimental group that are already recorded in the course of running A/B tests, which has practical and computational utility. However, the cross-validation procedure above requires the use of raw data. We now turn to the following question: if the raw data is unavailable, but summary statistics are, can we use these summary statistics to choose a threshold q?

Suppose that we have access to summary means $\lbrace (\bar{X}_g, \bar{Y}_g) \rbrace$ for each treatment j and the covariance matrix of $(\bar{X}, \bar{Y})$ conditional on Z = 0 which we denote by τ. We note that τ can be estimated very precisely from observational data or, in the case of the experimental meta-analysis just looking at covariances among known control groups. We assume that n _per is large enough such that the distributions of U and ϵ in groups of size $\frac{{n_{\text{per}}}}{2}$ are well approximated by the Gaussian $\mathcal {N}(0, \frac{\sigma ^2_i}{\frac{{n_{\text{per}}}}{2}}).$

To perform IVCV under these assumptions, we use a result from the literature on Monte Carlo [30, ch. 8]. If some vector X is distributed

This means if we know the observational covariance matrix τ then for every group g we can take the group level averages $(\bar{X}_g, \bar{Y}_g)$ and sample using the equation above to get $\bar{X}^1_g$ and $\bar{X}^2_g$ such that $\bar{X}^1_g + \bar{X}^2_g = 2\bar{X}_g$ . Since by the central limit theorem the generating Gaussian model is approximately correct, this procedure simulates the split required by IVCV without having access to the raw data. The algorithm is as follows:

Summary statistics instrumental variables cross-validation algorithm (sIVCV):

1. Start with data comprising of treatment group means $\lbrace (\bar{X}_g, \bar{Y}_g) \rbrace$ .
   
2. Use the covariance matrix to perform Monte Carlo sampling to simulate groups
   \[ \lbrace (X^1_i, Y^1_i, Z^1_i) \rbrace \]
   
   \[ \lbrace (X^2_i, Y^2_i, Z^2_i) \rbrace \]
   
   
   
3. Use the IVCV algorithm to set the hyperparameter using the simulated splits.
   
4. Estimate β using the selected hyperparameters on the full data set.
   

Figure 4: Left: Causal error (relative to a naive observational estimator) for the full l ₀-regularization path (solid black), TSLS (solid red), IVCV selected parameters (dashed purple) and Bayesian random effects model (dashed teal). IVCV outperforms all other estimation techniques. Right: Error in estimating causal effects for varying numbers of test groups K. IVCV is useful even with a relatively small meta-analysis, while TSLS exhibits asymptotic bias. With a very small number of test groups, the Oracle can actually underperform TSLS because of near collinearity.

We now evaluate the IVCV procedure empirically. True causal effects in real IV data are generally unobservable, so comparisons of methods usually lack a single number to validate against. Examples of the kinds of evaluations usually done in work on IV include comparing different procedures and showing that one yields estimates which are more ‘reasonable.’

Simulations allow us to know the true causal effects, but can lack realism. We strike a middle ground by using simulations where we set the causal effects ourselves but use real data to generate distributions for other variables. In our simulations we use a model given by

6.2 Results

In addition to the l ₀-regularized IV method and TSLS, we examine a Bayesian random effects model, as in Chamberlain and Imbens [12] but with a t, rather than Gaussian, distribution for the instruments. Formally we let

\[ \bar{Z} \sim t(d) \]

with a standard choice of prior

\[ d \sim \text{Gamma}(2, .2). \]

To tilt the evaluation against our procedure we also give the Bayesian model the true covariance matrix for $\bar{U}$ . To fit the Bayesian model we use Stan [ 11]. We compare the Bayesian random effects model and our regularized IV model to the infeasible Oracle estimator where the estimate of the first stage $\mathbb {E}[\bar{X} \mid \bar{Z}]$ is known with certainty.

Figure 4(A) shows the results for various dimensions of X for 1,000 simulations. Because of the high level of confounding in the observational data, the observational (OLS) estimates of the causal effect are highly biased, such that even the standard TSLS decreases our causal MSE by over $70\%.$

We see that the l ₀-regularization path (black line) reduces error compared with TSLS and, with high regularization, approaches the Oracle estimator. Furthermore, feasible selection of this hyperparameter using IVCV leads to near optimal performance (purple line). The Bayesian random effects model can reduce bias, but substantially increases variance and thus MSE.

We also look at how large the collection of experimental groups needs to be to see advantages of a regularized estimator relative to a TSLS procedure.

We repeat the TSLS, Oracle, and l ₀-regularization with IVCV analyses in 100 simulations with smaller K (Figure 4(B)) for the case of the 7 dimensional X. Intuitively, what is important is the relative size of the tails of the distribution of the latent treatment effects $\bar{Z}$ . As the tails get fatter, fewer experiments are required to get draws from the more extreme components of the mixture. We see that in this realistic case where $\bar{Z}$ is determined using a sampled set of Facebook A/B tests, feasible selection of the l ₀-regularization hyperparameter using IVCV outperforms TSLS substantially for many values of K. Thus, meta-analyses of even relatively small collections of experiments can be improved by the first-stage l ₀ regularization.

Figure 5: Performance of various IV estimation techniques under various first stage data generating assumptions (top = independent t, middle = Wishart t, bottom = correlated variances). We see that when the Z induced components of X are independent even for moderate dimensionality that the l ₀ regularization performs less well. However, as soon as there is any correlation the IVCV procedure performs much better than TSLS and can both under or over-perform the Bayesian random effects model.

In addition to the evaluation using a real data set. We also consider the IVCV procedure in several completely synthetic data sets. The completely synthetic experiments allows us to elucidate the important assumptions for our procedure to work while the real data-based experiment shows that these assumptions are indeed satisfied in real world conditions.

We consider the same exact model as in the previous section except that we generate the first stage effects $\bar{Z}$ from a known parametric distribution and let U be normal. First, we consider

1. $\bar{Z}$ drawn from independent t with 3 degrees of freedom
   
2. $\bar{Z}$ drawn from t with 3 degrees of freedom and covariance matrix drawn from inverse Wishart (a standard prior for covariance matrices) with 10 × dim(X) degrees of freedom
   
3. $\bar{Z}$ generated by first drawing σ ² from an inverse Gamma distribution and then $\bar{Z}$ drawn from multivariate normal with mean 0 and covariance matrix σ ² I
   

Note that in case 1 effects are axis aligned while in the next two larger values of one dimension can predict more extreme values of Z (and X) on another dimension. The final setup corresponds to the case where components are mean-uncorrelated but have correlated variance. This is the multivariate analog of our motivating example where some A/B tests are strong explorations of the parameter spaces and others are micro-optimizations at the margin. Note that the marginal distribution for each dimension is, in all cases, a t distribution with 3 degrees of freedom (since the t can be written as a mixture of normals drawn from the inverse gamma).

Figure 5 shows the results of applying IVCV to the data generating processes above (top = independent t, middle = Wishart t, bottom = correlated variances). We restrict to dim(X) ∈ {2, 4} because it is sufficient to illustrate our main points. We see that in the independent t case the IVCV procedure (and indeed our multivariate l ₀ regularization) can underperform the Bayesian random effects model and fail to substantially improve on TSLS. This happens because in the independent t case there is a high probability that a single dimension is extreme enough to pass the regularization threshold and thus even strong regularization does not necessarily remove bias. On the other hand, when outcomes are correlated (or their variances are) we see that multivariate IVCV performs well because being extreme in one X component predicts having extreme outcomes in other components. Note that the fact that IVCV performs well in the distribution generated by real A/B tests suggests that real world A/B test effect variance is correlated within A/B test — ie. if a test moves one metric by a large amount, it likely moves others by a large amount.

Most analysis of randomized experiments, whether in academia, business, or public policy tends to look at each test in isolation. When meta-analyses of experiments are conducted, the goal is usually either to pool data about multiple instances of the same intervention or to find heterogeneity in the effects of interventions across settings. We instead propose that combining many experiments can help us learn richer causal structures. IV methods give a way of doing this pooling. We have shown that in such situations using easily-implemented l ₀ regularization in the first stage can lead to much better estimates of the causal effects, and thus better predictions about interventions, than using standard TSLS methods.

We expand on the literature which uses multi-condition experiments as instruments [13, 16]. Such analyses usually feature a smaller number of experimental groups and a single causal variable. Our work contributes to the growing literature merging experimental techniques with methods from machine learning to allow deeper analyses than has been traditionally possible [4, 9, 14, 18, 32, 34]. In addition, our procedure is intended to be simple to implement and not require strong assumptions about the data-generating process.

Our work is also related to research on IV estimation with weak instruments [35, 36, 37]. In addition, we also contribute to existing research on regularized IV estimation [8, 12, 19]. In the case of univariate X and disjoint groups as instruments, the post-lasso method in Belloni et al. [8] coincides with the proposed l ₀ regularization, however in the case where X is a vector, it does not.

Recently, active learning in the form of bandit optimization [15, 28] and reinforcement learning [38] have become quite popular in the AI community. Such approaches can be used in many of the same contexts as the IV analysis we have discussed here, and so may appear to be substitutes. However, we note there are many important differences the largest of which is that RL and bandit approaches try to perform policy optimization rather than learning of a causal graph. For this reason, often why RL/bandit estimated policies work can be hard to understand. This is in contrast to explicit causal models (e.g., the linear model described above) which can be stated explicitly and in terms that are more natural to human decision-makers [24]. On the other hand, RL/bandit approaches have advantages in that they are explicitly online whereas many causal inference procedures (including the one we have described here) are ‘batch’ procedures that assume that data collection is a passive enterprise, separate from analysis. There is growing interest in combining these approaches [27] and we think that future work would benefit greatly from this fusion of thought.