Robust and Efficient Estimation of Potential Outcome Means Under Random Assignment

In the past several decades, the potential outcomes framework has become a staple of causal inference in statistics, econometrics, and related fields. Envisioning each unit in a population under different states of intervention or treatment allows one to define treatment or causal effects without reference to a model. One merely needs the potential outcome (PO) means for subpopulations.

When interventions are randomized – as in a clinical trial [Hirano and Imbens (2001)], assignment to participate in a job training program [Calónico and Smith (2017)], receiving a private school voucher [Angrist et al. (2006)], or contingent valuation studies where different bid values are randomized among individuals [Carson et al. (2004)] – one can simply use the subsample means for each treatment level in order to obtain unbiased and consistent estimators of the PO means. In some cases, the precision of the subsample means will be sufficient. Nevertheless, with the availability of good predictors of the outcome or response, it is natural to think that the precision can be improved, thereby shrinking confidence intervals and making conclusions about interventions more reliable.

In this paper we build on Negi and Wooldridge (2021), NW (2021) hereafter, who studied the problem of estimating the average treatment effect (ATE) under random assignment with one control and one treatment group. In the context of random sampling, NW (2021) showed that performing separate linear regressions for the two groups while estimating the ATE never does worse, asymptotically, than the simple difference in means estimator or a pooled regression adjustment estimator. These findings are complementary to Lin (2013), who used a finite population framework to study linear regression adjustment and arrived at essentially the same conclusions. NW (2021) also characterized the class of nonlinear regression adjustment methods that produce consistent estimators of the ATE without any additional assumptions (except regularity conditions).

In the current paper, we allow for $G\geq 2$ treatment levels and study the problem of joint estimation of the vector of PO means. We assume that the assignment to treatment is random – independent of both potential outcomes and observed predictors of the POs. In terms of the assignment mechanism, we assume a Bernoulli scheme whereby each unit has an equally likely probability of being assigned to a given treatment level, $g$. Importantly, other than substantive restrictions such as random assignment, $i.i.d$ sampling, and stable unit treatment value assumption, we impose only standard regularity conditions (such as finite second moments of the covariates). In other words, the regression adjustment (RA) estimators are consistent under essentially the same assumptions as the subsample means estimator with, generally, smaller asymptotic variances. Interestingly, even if the predictors are unhelpful, or the slopes in the linear projections are the same across all groups, no asymptotic efficiency is lost from using the most general RA method.

The extension from the binary treatment case to more than two treatment levels is nontrivial, both in terms of derivations and scope of applications. In the binary case, NW (2021) compared variances of estimators of the scalar ATE. In the current paper, we show that differences in asymptotic variance matrices of the PO estimators are positive semi-definite. Naturally, this result implies that any linear functions of the PO means are more efficiently estimated using separate RA. In addition, nonlinear functions of PO means –such as ratios – are also more efficiently estimated. Even in the binary treatment case, our results significantly improve over NW (2021).¹¹1It is not enough, even in the $G=2$ case, to use NW (2021) for each of the pairwise difference in means. Knowing that each pairwise ATE is more efficiently estimated does not imply that all linear combinations are.

We also extend the nonlinear RA results in NW (2021) to this general PO framework with $G$ treatment levels. Gail et al. (1984), who study nonlinear regression for the binary treatment case, have cautioned against using certain nonlinear regression models as being biased for the ATE. Imbens and Rubin (2015) also leave the impression that nonlinear RA will be more efficient but at the cost of consistency, and therefore should be avoided. We show that for particular kinds of responses– the leading cases being binary, fractional, and nonnegative – it is possible to consistently estimate PO means using pooled and separate RA by combining results from Quasi Maximum Likelihood Estimation (QMLE) in the linear exponential family (LEF) with features of the canonical link function.²²2In the generalized linear model literature, a link function $m^{-1}(\cdot)$ relates the conditional mean of the outcome to a linear predictor such that $\mathbb{E}[Y(g)|\mathbf{X}]=m(\mathbf{X}\bm{\beta})$. A canonical link is a special link function which guarantees that the model $m(\mathbf{X}\bm{\beta})$ fits the overall mean of the outcome without necessarily being correctly specified. While Gail et al. (1984) discussed maximum likelihood estimation in the exponential family, they do not distinguish between canonical and non-canonical link functions. As shown in NW (2021), and extended to the multiple treatment case here, one needs to use the mean function associated with the canonical link to ensure consistency without imposing additional restrictions. Guo and Basse (2023) reach a similar conclusion in the binary treatment case studied in NW (2021), showing that generalized linear models with canonical links are examples of imputation models that are prediction unbiased for the sample average treatment effect under the design-based framework. Unlike the linear RA case, we do not have a general result that shows separate nonlinear RA is asymptotically more efficient than the subsample means estimator. However, if we add the assumption that the conditional mean functions are correctly specified, we show that doing separate RA is (weakly) more efficient than the subsample means estimator.³³3Cohen and Fogarty (2024) propose a two-step calibrated Generalized Oaxaca Blinder estimator that is efficient even under misspecification. We further establish that separate RA attains the semiparametric efficiency bound under correct conditional mean specification.

A Monte Carlo exercise substantiates our theoretical results.⁴⁴4We consider three different population models. The outcomes are generated to be either fractional, non-negative, or continuous with an unrestricted support. The results for fractional and non-negative outcomes can be found in section B of the online appendix. We find that RA estimators, both linear and nonlinear, improve over subsample means in terms of standard deviations and have small biases, even with fairly small sample sizes. Not surprisingly, the magnitude of the precision gains with RA methods generally depends on how well the covariates predict the potential outcomes. Finally, we conclude with an empirical application that uses the RA methods to estimate the lower bound mean willingness to pay for an oil spill prevention program along California’s coast.

In terms of contribution, this paper provides a unified framework for studying regression adjustment in experiments, allowing for multiple treatments under the infinite population (or superpopulation) setting. Besides studying linear and nonlinear RA methods and making efficiency comparisons, we also characterize the semiparametric efficiency bound for all consistent and asymptotically linear estimators of the PO means and show that separate RA attains this bound under correct conditional mean specification. Because our results compare asymptotic variance matrices for different PO mean estimators, the efficiency results derived here extend to linear and (smooth) nonlinear functions of the PO means. We also fill an important gap in the regression adjustment literature by studying pooled nonlinear methods which have not been discussed elsewhere. We show how nonlinear pooled RA is consistent whenever nonlinear separate RA is, making it an important alternative when the researcher lacks sufficient degrees of freedom for estimating separate slopes. On the practical side, the proposed RA estimators are easy to implement with common statistical software like Stata (see online appendix D for reference).

Related Literature: This paper contributes to a vast and rich literature on regression adjustment in experiments. In the design-based framework, Freedman (2008a, b) highlighted that adjusting for covariates additively in a regression is not guaranteed to produce efficiency gains in experiments compared to the simple difference in means estimator. For a binary treatment, this debate was followed up in Lin (2013) which showed that separate RA is weakly more efficient than both the pooled RA and simple difference in means estimators of the ATE.⁵⁵5More recently, Chang, Middleton, and Aronow (2024) propose exact bias-corrected estimators to alleviate the finite sample concerns of RA. This is different from the focus of our paper which is interested in improving efficiency with RA and assumes that sample sizes are large enough to ensure consistent estimation of the PO means. A similar argument for bias correction in high dimensional settings can be found in Lei and Ding (2021) and Chiang, Matsushita, and Otsu (2023).

Among those adopting a superpopulation approach with a binary treatment, Yang and Tsiatis (2001) and Leon, Tsiatis, and Davidian (2003) focus on a pre-test post-test trial with no baseline covariates. Much like the current paper, specialized to the binary treatment case with ATE as the parameter of interest, Yang and Tsiatis (2001) show that the interacted estimator is the most efficient among simple difference in means and pooled RA estimators. Leon et al. (2003) identify the most efficient ATE estimator using the semiparametric theory of Robins et al. (1994) but do not provide any efficiency comparisons between commonly employed estimators of the ATE. Tsiatis et al. (2008) characterize the class of all consistent and asymptotically normal (CAN) estimators of the ATE. They further show that such estimators have a common form which can be used to rank the different linear estimators in terms of efficiency. The current paper (with $G=2$) and NW (2021) arrive at the same efficiency results with respect to linear RA as Tsiatis et al. (2008). In addition, the former explicitly discuss consistency with nonlinear RA (separate and pooled) involving quasi-log-likelihoods and canonical links from the generalized linear model literature. Furthermore, this paper establishes the asymptotic efficiency of separate RA relative to subsample averages when the conditional mean functions are correctly specified. Pitkin et al. (2013) compare the asymptotic variances of simple difference in means and linear RA (pooled and separate) under essentially a similar assumption as this paper wherein the conditional mean functions are not assumed to be linear. However, Pitkin et al. (2013) neither discuss nonlinear RA nor conditions under which it is efficient relative to simple difference in means. Similar to Leon et al. (2003), Zhang, Tsiatis, and Davidian (2008) also study regression adjustment from a semiparametric perspective. However, their focus is on CAN estimators of pairwise treatment comparisons in a multiarmed trial. As mentioned above, the efficiency results derived in this paper extend to linear and nonlinear functions of the PO means, including pairwise ATEs. The current paper also explicitly ranks the different linear and nonlinear RA estimators, which is not explored in Zhang et al. (2008). We further expand our set of results to include the semiparametric efficiency bound for all CAN estimators of the PO means and establish that SRA of the QML variety attains the efficiency bound when the conditional means are correctly specified. Another paper that studies nonlinear RA for estimating the ATE is Rosenblum and Van Der Laan (2010). The efficiency results discussed there are nested within this paper, which additionally considers pooled nonlinear methods that have not been discussed elsewhere. An added advantage of the nonlinear estimators proposed here is that they are computationally simple and easy to implement.

In the design-based setting, Guo and Basse (2023) and Cohen and Fogarty (2024) also discuss nonlinear RA in binary experiments. The former proposes ATE estimators which impute the missing potential outcomes using nonlinear models that are prediction unbiased. This is similar to the “mean-fitting” property discussed here and in NW (2021) which is crucial for consistent estimation of the PO means. This property is guaranteed by choosing appropriate combinations of canonical link and quasi-log-likelihood functions in the LEF estimated using QMLE. Cohen and Fogarty (2024) extend Guo and Basse (2023)’s efficiency result to the case of mean misspecification.⁶⁶6They argue that a second-step calibration of the Generalized Oaxaca Blinder estimator produces an estimator of the ATE that is efficient even under misspecification of the nonlinear model.

A recent strand also looks at multiarmed treatments in the context of factorial experiments under the design-based framework [Zhao and Ding (2022a, b), Zhao and Ding (2023), Pashley and Bind (2023)]. Such experiments are designed to accommodate multiple factors of interest in a single experiment and define treatment levels as all possible combinations of the factors involved. The typical concern in this literature is degree of freedom conservation while estimating all main and interaction effects of these factors simultaneously. An exception is Zhao and Ding (2023) which discusses linear RA for multiarmed experiments (factorial experiments being a special case). The efficiency results concerning linear RA in Zhao and Ding (2023) have a direct analogue in the current paper under the superpoulation framework. For more complex experimental designs, see Cytrynbaum (2023), Jiang et al. (2023), Roth and Sant’Anna (2023), and Chang (2023).

The rest of the paper is organized as follows. Section 2 describes the potential outcomes framework with $G$ treatment levels along with a discussion of the main assumptions. Section 3 presents the asymptotic variances of the linear RA estimators, whereas section 4 ranks them in terms of asymptotic efficiency. Section 5 considers a class of nonlinear RA estimators that ensure consistency for the PO means without imposing additional assumptions. Section 6 characterizes the semiparametric efficiency bound for all regular estimators of the PO means and establishes local efficiency of separate RA. Section 7 constructs a simulation exercise for studying the finite sample behavior of the RA estimators. Section 8 applies the RA methods to study willingness to pay for an oil spill prevention program, and section 9 concludes.

We use the standard potential outcomes framework, also known as the Neyman-Rubin causal model. The goal is to estimate the population means of $G$ potential (counterfactual) outcomes, $Y(g)$, $g=1,...,G$. Define $\mu_{g}=\mathbb{E}\left[Y(g)\right]$ for each $g=1,...,G.$

The vector of assignment indicators is $\mathbf{W}=(W_{1},...,W_{G})$, where each $W_{g}$ is binary and $W_{1}+W_{2}+\cdots+W_{G}=1$. In other words, the groups are exhaustive and mutually exclusive. The setup applies to many situations, including the standard treatment-control setup with $G=2$, multiple treatment levels (with $g=1$ as the control group), and in contingent valuation studies where subjects are presented with a set of $G$ prices or bid values.

Next, let $\mathbf{X}=(X_{1},X_{2},...,X_{K})$ be a vector of observed covariates of dimension $K$ which is assumed to be fixed in our asymptotic analysis. With respect to the assignment process, we make the following assumption.

Assumption 1 puts us in the framework of experimental interventions. We also assume that each group, $g$, has a positive probability of being assigned such that $\rho_{1}+\rho_{2}+\cdots+\rho_{G}=1$.

The i.i.d assumption is not the only one we can make. For example, we could allow for a sampling-without-replacement scheme given a fixed sample size $N$. This would complicate the analysis because it generates a slight correlation across draws from the population. As discussed in NW (2021), Assumption 2 is traditional in studying the asymptotic properties of estimators and is realistic as an approximation in cases where a sample is taken from a large population. It also forces us to account for the sampling error in $\mathbf{\bar{X}}$, as an estimator of $\bm{\mu}_{\mathbf{X}}=\mathbb{E}\left(\mathbf{X}\right)$.

A different approach is to assume that there is no sampling uncertainty and instead assume we observe the entire population. In the design–based approach, all uncertainty is in the assignment of the treatment indicators, $\mathbf{W}=(W_{1},...,W_{G})$. Among others, Freedman (2008a) and Lin (2013) take this approach in studying linear regression adjustment with a binary treatment. In many examples in economics (including the empirical example in our supplement), one does not observe the entire population and so it seems natural to study the random sampling framework. Our results for the random sampling case can be used to further study the efficiency issue in a framework that allows both sampling and assignment uncertainty, as in Abadie et al. (2020). Such a framework will affect inference only in cases where the sample is a substantial fraction of the population. If the sample size is small relative to the population size, the finite-population adjustments of the kind derived in Abadie et al. (2020) would have a trivial effect.

Given Assumption 2, for each draw $i$ from the population we only observe

and so the data are $\left\{\left(\mathbf{W}_{i},Y_{i},\mathbf{X}_{i}\right):i=1,2,...,N\right\}$. Definition of population quantities only requires us to use the random vector $\left(\mathbf{W},Y,\mathbf{X}\right)$, which represents the population.

Along with assumptions 1 and 2, we also assume that the stable unit treatment value assumption (SUTVA) holds. This implies that there are no spillovers or hidden variations of the treatments. Random assignment, $i.i.d$ sampling, and SUTVA are the only substantive restrictions used in this paper and will be maintained throughout. Subsequently, we assume that linear projections exist and that the central limit theorem holds for properly standardized sample averages of i.i.d random vectors. Therefore, we are implicitly imposing at least finite second moment assumptions on each $Y(g)$ and $X_{j}$. We do not make this explicit in what follows.

In this section we discuss the asymptotic variances of three estimators: the subsample means, separate regression adjustment, and pooled regression adjustment.⁷⁷7The asymptotic representation proofs of the three estimators can be found in section A of the online appendix.

We now take the representations given above and use them to compare the asymptotic variances of the three estimators. It is helpful to summarize the conclusions reached in Section 3:

where $\mathbf{L}_{i}$, $\mathbf{Q}_{i}$, $\mathbf{K}_{i}$, and $\mathbf{F}_{i}$ are defined in (2), (4) and (5), respectively.

We now discuss a class of nonlinear regression adjustment methods that preserve consistency without adding additional assumptions other than the ones discussed in section 2 (and weak regularity conditions). In particular, we extend the setup in NW (2021) to more than two treatment levels.

Typically, when the outcome is discrete or has limited support, nonlinear functions of $\mathbb{E}[Y(g)|\mathbf{X}]$ may be more appropriate and may provide a better approximation to the true conditional mean. However, consistency with nonlinear models is closely tied to functional form assumptions. In this section, we discuss nonlinear methods that are consistent for $\mathbb{E}[Y(g)]$ despite possible misspecification in $\mathbb{E}[Y(g)|\mathbf{X}]$. Not surprisingly, using a canonical link function in the context of QMLE in the linear exponential family plays a key role. We show that both separate and pooled nonlinear methods are consistent provided we choose the mean functions and objective functions appropriately. Unlike in the linear case, we can only show that SRA improves over the subsample means estimator when the conditional means are correctly specified.

In this section, we generalize the efficiency result with regression adjustment in experiments for multiple treatments by providing an efficiency benchmark for estimating the vector of PO means. We characterize the efficient influence function and the semiparametric efficiency bound (SEB) for all regular and asymptotically linear estimators of $\bm{\mu}$. The SEB is useful as it provides the semiparametric analogue of the Cramer-Rao lower bound for parametric models. It gives us a standard against which to compare the asymptotic variance of any regular estimator of the potential outcome means, including the RA estimators studied in the previous sections.

In our case, the most efficient estimator for $\bm{\mu}$ has the following asymptotically linear form

such that $\mathbb{E}[\bm{\psi}_{i}]=\bm{0}$ and $\mathbb{E}[\bm{\psi}_{i}\bm{\psi}_{i}^{\prime}]$ is finite and non-singular. Then, the SEB is given by $\mathbf{V_{\ast}}=\mathbb{E}[\bm{\psi}_{i}\bm{\psi}_{i}^{\prime}]$ where the influence function is

This result follows from Cattaneo (2010), who studies efficient semiparametric estimation of treatment effects under the assumption of unconfoundedness. Given that we have the stronger assumption of randomized treatment, this essentially leads us to the same moment conditions as Cattaneo (2010) with the exception that the propensity score here is a constant.

The proof can be found in Appendix A. The above result establishes that the separate regression adjustment estimator achieves the semiparametric efficiency bound if the true conditional mean of the POs is correctly specified. This implies that running separate linear regressions for different treatment level, as advocated in section 4, is optimally efficient if the mean function, $m(\cdot)$, is truly linear. However, if $m(\cdot)$ is truly logistic or exponential, then nonlinear SRA discussed in section 5, is optimally efficient for estimating the PO means.

This section studies the finite sample properties of the different estimators of the PO means, namely, SM, PRA, SRA, and their nonlinear counterparts. For the simulations, we generate a population of one million observations and mimic the asymptotic setting of random sampling from an “infinite” population. The empirical distributions of the RA estimators are simulated for sample sizes $N\in\{500,1000,5000\}$ by randomly drawing the data vector $\left\{\left(\mathbf{W}_{i},Y_{i},\mathbf{X}_{i}\right):i=1,2,...,N\right\}$, without replacement, ten thousand times from the population. We consider three different population models, each of which uses a unique data generating process for the potential outcomes, but we report results only for the linear model. Results for fractional and non-negative outcomes can be found in online appendix B.⁹⁹9See Tables B.1-B.4 for results on fractional and non-negative outcomes.

To simulate multiple treatments, we consider potential outcomes, $Y(g)$, corresponding to three treatment states, $g=1,2,3$. Hence, $G=3$ for all population models. In each of the populations, the treatment vector $\mathbf{W}=\left(W_{1},W_{2},W_{3}\right)$ is generated with probability mass function defined by $\rho_{g}$ such that there is an equally likely probability of being assigned to a particular treatment group.