Experimental Evaluation of Individualized Treatment Rules

TL;DR

This paper introduces new metrics, PAPE and AUPEC, for evaluating individualized treatment rules using experimental data, providing exact variance estimates without modeling assumptions, applicable to complex algorithms.

Contribution

It proposes novel evaluation metrics (PAPE and AUPEC) and their variance estimation methods for assessing ITRs, including when using the same data for estimation and evaluation.

Findings

PAPE and AUPEC effectively measure ITR performance.

Variance formulas enable precise evaluation without resampling.

Method applicable to complex machine learning-based ITRs.

Abstract

The increasing availability of individual-level data has led to numerous applications of individualized (or personalized) treatment rules (ITRs). Policy makers often wish to empirically evaluate ITRs and compare their relative performance before implementing them in a target population. We propose a new evaluation metric, the population average prescriptive effect (PAPE). The PAPE compares the performance of ITR with that of non-individualized treatment rule, which randomly treats the same proportion of units. Averaging the PAPE over a range of budget constraints yields our second evaluation metric, the area under the prescriptive effect curve (AUPEC). The AUPEC represents an overall performance measure for evaluation, like the area under the receiver and operating characteristic curve (AUROC) does for classification, and is a generalization of the QINI coefficient utilized in uplift modeling. We use Neyman's repeated sampling framework to estimate the PAPE and AUPEC and derive their exact finite-sample variances based on random sampling of units and random assignment of treatment. We extend our methodology to a common setting, in which the same experimental data is used to both estimate and evaluate ITRs. In this case, our variance calculation incorporates the additional uncertainty due to random splits of data used for cross-validation. The proposed evaluation metrics can be estimated without requiring modeling assumptions, asymptotic approximation, or resampling methods. As a result, it is applicable to any ITR including those based on complex machine learning algorithms. The open-source software package is available for implementing the proposed methodology.

Figures

Tables

Table 1. Table 1 : The Results of the Simulation Study for Fixed Individualized Treatment Rules (ITRs). The table presents the bias and standard deviation of each estimator as well as the coverage of its 95% confidence intervals under the “Low treatment effect” and “High treatment effect” scenarios. The first three estimators shown here are for BART f 𝑓 f : Population Average Prescription effect (PAPE; τ ^ f subscript ^ 𝜏 𝑓 \hat{\tau}_{f} ), PAPE with a budget constraint of 20% treatment proportion ( τ ^ f ​ ( c 0.2 ) subscript ^ 𝜏 𝑓 subscript 𝑐 0.2 \hat{\tau}_{f}(c_{0.2}) ), Area Under the Prescriptive Effect Curve (AUPEC; Γ ^ f subscript ^ Γ 𝑓 \widehat{\Gamma}_{f} ). We also present the results for the difference in the PAPE between BART and Causal Forest g 𝑔 g ( Δ ^ 0.2 ​ ( f , g ) subscript ^ Δ 0.2 𝑓 𝑔 \widehat{\Delta}_{0.2}(f,g) ) and between BART and LASSO h ℎ h ( Δ ^ 0.2 ​ ( f , g ) subscript ^ Δ 0.2 𝑓 𝑔 \widehat{\Delta}_{0.2}(f,g) ) under the budget constraint.

		$𝒏=\mathrm{𝟏𝟎𝟎}$			$𝒏=\mathrm{𝟓𝟎𝟎}$			$𝒏=\mathrm{𝟐𝟎𝟎𝟎}$
Estimator	truth	coverage	bias	s.d.	coverage	bias	s.d.	coverage	bias	s.d.
Low treatment effect
${\hat{\tau }}_f$	$0.066$	$94.3\%$	$0.005$	$0.124$	$96.2\%$	$0.001$	$0.053$	$95.1\%$	$0.001$	$0.026$
${\hat{\tau }}_f\left(c_{0.2}\right)$	$0.051$	$93.2$	$-0.002$	$0.109$	$94.4$	$0.001$	$0.046$	$95.2$	$0.002$	$0.021$
${\hat{\Gamma }}_f$	$0.053$	$95.3$	$0.001$	$0.106$	$95.1$	$0.001$	$0.045$	$94.8$	$-0.001$	$0.024$
${\hat{\Delta }}_{0.2}(f,g)$	$-0.022$	$94.0$	$0.006$	$0.122$	$95.4$	$0.002$	$0.051$	$96.0$	$0.000$	$0.026$
${\hat{\Delta }}_{0.2}(f,h)$	$-0.014$	$93.9$	$-0.001$	$0.131$	$94.9$	$-0.000$	$0.060$	$95.3$	$-0.000$	$0.030$
High treatment effect
${\hat{\tau }}_f$	$0.430$	$94.7\%$	$-0.000$	$0.163$	$95.7\%$	$0.000$	$0.064$	$94.4\%$	$-0.000$	$0.031$
${\hat{\tau }}_f\left(c_{0.2}\right)$	$0.356$	$94.7$	$0.004$	$0.159$	$95.7$	$0.002$	$0.072$	$95.8$	$0.000$	$0.035$
${\hat{\Gamma }}_f$	$0.363$	$94.3$	$-0.005$	$0.130$	$94.9$	$0.003$	$0.058$	$95.7$	$0.000$	$0.029$
${\hat{\Delta }}_{0.2}(f,g)$	$-0.000$	$96.9$	$0.008$	$0.151$	$97.9$	$-0.002$	$0.073$	$98.0$	$-0.000$	$0.026$
${\hat{\Delta }}_{0.2}(f,h)$	$0.000$	$94.7$	$-0.004$	$0.140$	$97.7$	$-0.001$	$0.065$	$96.6$	$0.000$	$0.033$

Table 2. Table 2 : The Results of the Simulation Study for Cross-validated ITR. The table presents the true value (truth) of each quantity along with the bias and standard deviation of each estimator as well as the coverage of its 95% confidence intervals under the “Low treatment effect” and “High treatment effect” scenarios. All of the results shown here are for LASSO.

	$𝒏=\mathrm{𝟏𝟎𝟎}$				$𝒏=\mathrm{𝟓𝟎𝟎}$				$𝒏=\mathrm{𝟐𝟎𝟎𝟎}$
Estimator	truth	coverage	bias	s.d.	truth	coverage	bias	s.d.	truth	coverage	bias	s.d.
Low Effect
${\hat{\lambda }}_F$	$0.073$	$96.4\%$	$0.001$	$0.216$	$0.095$	$96.7\%$	$0.002$	$0.100$	$0.112$	$97.2\%$	$0.002$	$0.046$
${\hat{\tau }}_F$	$0.021$	$94.6$	$-0.002$	$0.130$	$0.030$	$95.5$	$-0.002$	$0.052$	$0.032$	$94.4$	$-0.000$	$0.027$
${\hat{\tau }}_F\left(c_{0.2}\right)$	$0.023$	$95.4$	$-0.003$	$0.120$	$0.034$	$95.4$	$-0.002$	$0.057$	$0.043$	$96.8$	$0.001$	$0.029$
${\hat{\Gamma }}_F$	$0.009$	$98.2$	$0.002$	$0.117$	$0.029$	$96.8$	$-0.001$	$0.048$	$0.039$	$95.9$	$0.001$	$0.001$
High Effect
${\hat{\lambda }}_H$	$0.867$	$96.9\%$	$-0.007$	$0.261$	$0.875$	$96.5\%$	$-0.003$	$0.125$	$0.875$	$97.3\%$	$0.001$	$0.062$
${\hat{\tau }}_F$	$0.338$	$93.6$	$-0.000$	$0.171$	$0.358$	$93.0$	$0.000$	$0.093$	$0.391$	$95.3$	$0.001$	$0.041$
${\hat{\tau }}_F\left(c_{0.2}\right)$	$0.341$	$94.8$	$-0.002$	$0.170$	$0.356$	$96.2$	$-0.005$	$0.075$	$0.356$	$95.8$	$0.001$	$0.037$
${\hat{\Gamma }}_F$	$0.344$	$98.5$	$0.001$	$0.126$	$0.362$	$98.9$	$0.005$	$0.053$	$0.363$	$99.0$	$0.001$	$0.026$

Table 3. Table 3 : The Estimated Population Average Prescription Effect (PAPE) for BART, Causal Forest, and LASSO with and without a Budget Constraint. We estimate the PAPE for fixed and estimated individualized treatment rules (ITRs). The fixed ITRs are based on the training (70%) and test data (30%), whereas the estimated ITRs are based on 5 fold cross-validation. In addition, the average treatment effect estimates using the entire dataset. For each of the three outcomes, the point estimate, the standard error, and the average proportion treated are shown. The budget constraint considered here implies that the maximum proportion treated is 20%.

	BART			Causal Forest			LASSO
	est.	s.e.	treated	est.	s.e.	treated	est.	s.e.	treated
Fixed ITR
No budget constraint
Reading	$0$	$0$	100%	$-0.38$	$1.14$	84.3%	$-0.41$	$1.10$	84.4%
Math	$0.52$	$1.09$	86.7	$0.09$	$1.18$	80.3	$1.73$	$1.25$	78.7
Writing	$-0.32$	$0.72$	92.7	$-0.70$	$1.18$	78.0	$-0.30$	$1.26$	80.0
\hdashline     Budget constraint
Reading	$-0.89$	$1.30$	20	$0.66$	$1.23$	20	$-1.17$	$1.18$	20
Math	$0.70$	$1.25$	20	$2.57$	$1.29$	20	$1.25$	$1.32$	20
Writing	$2.60$	$1.17$	20	$2.98$	$1.18$	20	$0.28$	$1.19$	20
Estimated ITR
No budget constraint
Reading	$0.19$	$0.37$	99.3%	$0.31$	$0.77$	86.6%	$0.32$	$0.53$	87.6%
Math	$0.92$	$0.75$	84.7	$2.29$	$0.80$	79.1	$1.52$	$1.60$	75.2
Writing	$1.12$	$0.86$	88.0	$1.43$	$0.71$	67.4	$0.05$	$1.37$	74.8
\hdashline     Budget constraint
Reading	$1.55$	$1.05$	20	$0.40$	$0.69$	20	$-0.15$	$1.41$	20
Math	$2.28$	$1.15$	20	$1.84$	$0.73$	20	$1.50$	$1.48$	20
Writing	$2.31$	$0.66$	20	$1.90$	$0.64$	20	$-0.47$	$1.34$	20

Table 4. Table 4 : The Estimated Population Average Prescriptive Effect Difference (PAPD) between Causal Forest, BART, and LASSO under a Budget Constraint of 20%. The point estimates and 95% confidence intervals are shown.

	Causal Forest				BART
	vs. BART		vs. LASSO		vs. LASSO
	est.	$95\%$ CI	est.	$95\%$ CI	est.	$95\%$ CI
Fixed ITR
Math	$1.55$	[$-$0.35, 3.45]	$1.83$	[$-$0.50, 4.16]	$0.28$	[$-$2.39, 2.95]
Reading	$1.86$	[$-$0.79, 4.51]	$1.31$	[$-$1.49, 4.11]	$-0.55$	[$-$4.02, 2.92]
Writing	$0.38$	[$-$1.66, 2.42]	$2.69$	[$-$0.27, 5.65]	$2.32$	[$-$0.53, 5.15]
Estimated ITR
Reading	$-1.15$	[$-$3.99, 1.69]	$0.55$	[$-$1.05, 2.15]	$1.70$	[$-$0.90, 4.30]
Math	$-0.43$	[$-$2.57, 3.43]	$0.34$	[$-$1.32, 2.00]	$0.77$	[$-$1.99, 3.53]
Writing	$-0.41$	[$-$1.63, 0.80]	$2.37$	[0.76, 3.98]	$2.79$	[1.32, 4.26]

Table 5. Table A1 : A Numerical Example for Binary Treatment Assignment and Outcomes

Individual	$T_i$	$f\left({𝐗}_i\right)$	$Y_i$	$Y_i\left(0\right)$	$Y_i\left(1\right)$
A	1	1	2	0	2
B	1	0	3	1	3
C	0	0	$-$1	$-$1	$-$1
D	0	1	1	1	0
E	1	0	3	0	3

Table 6. Table A2 : The Estimated Population Average Prescription Effect (PAPE) for Bayesian Causal Forest (BCF), Causal Forest, and R-Learner with and without a Budget Constraint. For each of the three outcomes, the point estimate, the standard error, and the proportion treated are shown. The budget constraint considered here implies that the maximum proportion treated is 20%.

	BCF			Causal Forest			R-Learner
	est.	s.e.	treated	est.	s.e.	treated	est.	s.e.	treated
No budget constraint	$-0.104$	$0.128$	48.4%	$-0.349$	$0.137$	47.5%	$0$	$0$	100%
20% Budget constraint	$-0.02$	$0.121$	20%	$-0.120$	$0.107$	20%	$0$	$0.104$	20%

Equations

$$f:\mathcal{X}⟶\{0,1\}.$$

$${\lambda }_f=\mathbb{E}\{Y_i(f({\mathbf{X}}_i))\}.$$

$${\tau }_f=\mathbb{E}\{Y_i(f({\mathbf{X}}_i))-p_f{Y}_i(1)-(1-p_f)Y_i(0)\}.$$

$$f({\mathbf{X}}_i,c)=1\{s({\mathbf{X}}_i)>c\}.$$

$$c_p(f)=\inf \{c\in \mathbb{R}:\mathrm{Pr}(f({\mathbf{X}}_i,c)=1)\le p\}.$$

$$s(\mathbf{x})=\mathbb{E}(Y_i(1)-Y_i(0)\mid {\mathbf{X}}_i=\mathbf{x}).$$

$${\tau }_{fp}=\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c_p(f)))-p{Y}_i(1)-(1-p)Y_i(0)\}.$$

$$\Delta (f,g)={\lambda }_f-{\lambda }_g=\mathbb{E}\{Y_i(f({\mathbf{X}}_i))-Y_i(g({\mathbf{X}}_i))\}.$$

$${\Delta }_p(f,g)={\tau }_{fp}-{\tau }_{gp}=\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c_p(f)))-Y_i(g({\mathbf{X}}_i,c_p(g)))\}.$$

$${\Gamma }_f={\int }_0^{p_f}\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c_p(f)))\}\mathrm{d}p+(1-p_f)\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c^{\ast }))\}-\frac{1}{2}\mathbb{E}(Y_i(0)+Y_i(1)),$$

$$\text{QINI}=n\left({\int }_0^{1}p\mathbb{E}\{Y_i(1)-Y_i(0)\mid f({\mathbf{X}}_i,c_p(f))=1\}\mathrm{d}p-\frac{1}{2}\mathbb{E}(Y_i(1)-Y_i(0))\right).$$

$$n\left({\int }_0^1\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c_p(f)))\}\mathrm{d}p-\frac{1}{2}\mathbb{E}(Y_i(0)+Y_i(1))\right)$$

$${\Gamma }_f^{\ast }=\frac{1}{\tau }\left[{\int }_0^{p_f}\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c_p(f)))\}\mathrm{d}p+(1-p_f)\mathbb{E}\{Y_i(f({\mathbf{X}}_i,c^{\ast }))\}-\mathbb{E}(Y_i(0))\right]-\frac{1}{2}.$$

$$(Y_i(1),Y_i(0),{\mathbf{X}}_i)\stackrel{i.i.d.}{\sim }\mathcal{P}$$

$$\mathrm{Pr}(T_i=1\mid Y_i(1),Y_i(0),{\mathbf{X}}_i)=\frac{n_1}{n}$$

$${\hat{\tau }}_{fp}(\mathcal{Z})=\frac{1}{n_1}\sum _{i=1}^n{Y}_i{T}_{i}f({\mathbf{X}}_i,{\hat{c}}_p(f))+\frac{1}{n_0}\sum _{i=1}^n{Y}_i(1-T_i)(1-f({\mathbf{X}}_i,{\hat{c}}_p(f)))-\frac{p}{n_1}\sum _{i=1}^n{Y}_i{T}_i-\frac{1-p}{n_0}\sum _{i=1}^n{Y}_i(1-T_i)$$

$${\mathbb{P}}_{{\hat{c}}_p(f)}(|\mathbb{E}\{{\hat{\tau }}_{fp}(\mathcal{Z})-{\tau }_{fp}\mid {\hat{c}}_p(f)\}|\ge ϵ)$$

$${\gamma }_{fp}(ϵ)=\frac{ϵ}{{\max }_{c\in [c_p(f)-ϵ,c_p(f)+ϵ]}\mathbb{E}({\tau }_i\mid s({\mathbf{X}}_i)=c)}.$$

$$\mathbb{V}({\hat{\tau }}_{fp}(\mathcal{Z}))=\frac{\mathbb{E}(S_{fp1}^2)}{n_1}+\frac{\mathbb{E}(S_{fp0}^2)}{n_0}+\frac{\lfloor np\rfloor (n-\lfloor np\rfloor )}{n^2(n-1)}\left\{(2p-1){\kappa }_{f1}(p{)}^2-2p{\kappa }_{f1}(p){\kappa }_{f0}(p)\right\},$$

$$\widehat{\mathbb{E}(S_{fpt}^2)}$$

$$\widehat{{\kappa }_{ft}(p)}$$

$${\hat{\Delta }}_p(f,g,\mathcal{Z})=\frac{1}{n_1}\sum _{i=1}^n{Y}_i{T}_i\{f({\mathbf{X}}_i,{\hat{c}}_p(f))-g({\mathbf{X}}_i,{\hat{c}}_p(g))\}+\frac{1}{n_0}\sum _{i=1}^n{Y}_i(1-T_i)\{g({\mathbf{X}}_i,{\hat{c}}_p(g))-f({\mathbf{X}}_i,{\hat{c}}_p(f))\}.$$

$${\hat{\Gamma }}_f(\mathcal{Z})$$

$${\mathbb{P}}_{{\hat{p}}_f}(|\mathbb{E}({\hat{\Gamma }}_f(\mathcal{Z})-{\Gamma }_f\mid {\hat{p}}_f)|\ge ϵ)$$

$${\gamma }_{p_f}(ϵ)=\frac{ϵ}{2{\max }_{c\in [c^{\ast }-ϵ,c^{\ast }+ϵ]}\mathbb{E}({\tau }_i\mid s({\mathbf{X}}_i)=c)}.$$

$$\mathbb{V}({\hat{\Gamma }}_f(\mathcal{Z}))$$

$$\widehat{\mathbb{E}(S_{ft}^{\ast 2})}$$

$${\hat{\kappa }}_{ft}(z/n)$$

$${\hat{\Gamma }}_f^{\ast }(\mathcal{Z})$$

$$F:\mathcal{Z}\to \mathcal{S}.$$

$${\hat{f}}_{{\mathcal{Z}}^{tr}}({\mathbf{X}}_i,c_p({\hat{f}}_{{\mathcal{Z}}^{tr}}))=1\{{\hat{s}}_{{\mathcal{Z}}^{tr}}({\mathbf{X}}_i)>c_p({\hat{f}}_{{\mathcal{Z}}^{tr}})\},$$

Full text

Key Words: causal inference, heterogenous treatment effects, machine learning, precision medicine, uplift modeling

1 Introduction

In today’s data-rich society, the individualized (or personalized) treatment rules (ITRs), which assign different treatments to individuals based on their observed characteristics, play an essential role. Examples include personalized medicine and micro-targeting in business and political campaigns (e.g., Hamburg and Collins, 2010; Imai and Strauss, 2011). In the causal inference literature, a number of researchers have developed methods to estimate optimal ITRs using a variety of machine learning algorithms (see e.g., Qian and Murphy, 2011; Zhang et al., 2012; Fu et al., 2016; Luedtke and van der Laan, 2016a, b; Zhou et al., 2017; Athey and Wager, 2018; Kitagawa and Tetenov, 2018). In addition, applied researchers often use machine learning algorithms to estimate heterogeneous treatment effects and then construct ITRs based on the resulting estimates.

In this paper, we consider a common setting, in which a policy-maker wishes to experimentally evaluate the empirical performance of an ITR before implementing it in a target population. Such evaluation is also essential for comparing the efficacy of alternative ITRs. Specifically, we show how to use a randomized experiment for evaluating ITRs. We propose two new evaluation metrics. The first is the population average prescriptive effect (PAPE), which compares an ITR with a non-individualized treatment rule that randomly assigns the same proportion of units to the treatment condition. The PAPE represents the difference between the average outcome under the ITR and that under the random treatment rule. The key idea is that a well-performing ITR should outperform the random treatment rule, which does not utilize any individual-level information.

Averaging the PAPE over a range of budget constraints yields our second evaluation metric, the area under the prescriptive effect curve (AUPEC). Like the area under the receiver and operating characteristic curve (AUROC) for classification, the AUPEC represents an overall summary measure of how well an ITR performs over the random treatment rule that treats the same proportion of units.

We estimate these evaluation metrics using Neyman (1923)’s repeated sampling framework (see Imbens and Rubin, 2015, Chapter 6). An advantage of this approach is that it does not require any modeling assumption or asymptotic approximation. As a result, we can evaluate a broad class of ITRs including those based on complex machine learning algorithms. We show how to estimate the PAPE and AUPEC with a minimal amount of finite sample bias and derive the exact variance solely based on random sampling of units and random assignment of treatment.

We further extend this methodology to a common evaluation setting, in which the same experimental data is used to both estimate and evaluate ITRs. In this case, our finite-sample variance calculation is exact and directly incorporates the additional uncertainty due to random splits of data used for cross-validation. We implement the proposed methodology through an open-source R package.

Our simulation study demonstrates the accurate coverage of the proposed confidence intervals in small samples (Section 5). We also apply our methods to the Project STAR (Student-Teacher Achievement Ratio) experiment and compare the empirical performance of ITRs based on several popular methods (Section 6). Our evaluation approach addresses theoretical and practical difficulties of conducting reliable statistical inference for ITRs.

Relevant Literature.

A large number of existing studies have focused on the derivation of optimal ITRs that maximize the population average value. For example, Qian and Murphy (2011) use penalized least squares whereas Zhao et al. (2012) show how a support vector machine can be used to derive an optimal ITR. Another popular approach is based on doubly-robust estimation (e.g., Dudík et al., 2011; Zhang et al., 2012; Chakraborty et al., 2014; Jiang and Li, 2016; Athey and Wager, 2018; Kallus, 2018).

We propose a general methodology for empirically evaluating and comparing the performance of various ITRs including the ones proposed by these and other authors. While many of these methods come with uncertainty measures, even those that produce standard errors rely on asymptotic approximation, modeling assumptions, or resampling methods. In contrast, our methodology utilizes Neyman’s repeated sampling framework and does not require any of these assumptions or approximations.

There also exists a related literature on policy evaluation. Starting with Manski (2004), many studies focus on the derivation of regret bounds given a class of ITRs. For example, Kitagawa and Tetenov (2018) show that an ITR, which maximizes the empirical average value, is minimax optimal without a strong restriction on the class of ITRs, whereas Athey and Wager (2018) establish a regret bound for an ITR based on doubly-robust estimation in observational studies (see also Zhou et al., 2018). In addition, Luedtke and van der Laan (2016a, b) propose consistent estimators of the optimal average value even when an optimal ITR is not unique (see also Rai, 2018).

Our goal is different from these studies. We focus on statistical inference using the Neyman’s repeated sampling framework for the experimental evaluation of arbitrary ITRs including optimal or non-optimal and simple or complex ones. Our evaluation metric is also different from the existing metrics. In particular, to the best of our knowledge, we are the first to formally study the AUPEC as an AUROC-like summary measure for evaluation.

In contrast, much of the policy evaluation literature focus on the optimal average value, which is required to compute the regret of an ITR. Athey and Wager (2018) briefly discusses a quantity related to the PAPE in their empirical application, but this quantity evaluates an ITR against the treatment rule that randomly treats exactly one half of units rather than the same proportion as the one treated under the ITR. Although empirical studies in the campaign and marketing literatures have used “uplift modeling,” which is based on the PAPE (e.g., Imai and Strauss, 2011; Rzepakowski and Jaroszewicz, 2012; Gutierrez and Gérardy, 2016; Ascarza, 2018; Fifield, 2018), none develops formal estimation and inferential methods. We show that the AUPEC is a generalization of the QINI coefficient, which is a widely utilized statistic in uplift modeling (Radcliffe, 2007; Diemert et al., 2018). Thus, our theoretical results for the AUPEC apply directly to the QINI coefficient as well.

Another related literature is concerned with the estimation of heterogeneous effects. Researchers have explored the use of tree-based methods (e.g., Imai and Strauss, 2011; Athey and Imbens, 2016; Wager and Athey, 2018; Hahn et al., 2020), regularized regressions (e.g., Imai and Ratkovic, 2013; Künzel et al., 2018), and ensamble methods (e.g., van der Laan and Rose, 2011; Grimmer et al., 2017). In practice, the estimated heterogeneous treatment effects based on these machine learning algorithms are used to construct ITRs.

However, as Chernozhukov et al. (2019) point out, most machine learning algorithms, which require data-driven tuning parameters, cannot be regarded as consistent estimators of the conditional average treatment effect (CATE) unless strong assumptions are imposed. They propose a methodology to estimate heterogeneous treatment effects without such assumptions. Similar to theirs, our methodology does not depend on any modeling assumption and accounts for the uncertainty due to splitting of data. The key difference is that we focus on the evaluation of ITRs. In addition, our variance calculation is based on randomization and does not rely on asymptotic approximation.

Finally, Andrews et al. (2020) develops a conditional inference procedure, based on normal approximation, for the average value of the best-performing policy based on experimental or observational data. In contrast, we develop an unconditional exact inference for the difference in the average value between any pair of policies under a budget constraint. We also consider the evaluation of estimated policies based on the same experimental data using cross-validation whereas Andrews et al. focus on the evaluation of fixed policies.

2 Evaluation Metrics

In this section, we introduce our evaluation metrics. We first propose the population average prescriptive effect (PAPE) that, unlike the population average value, adjusts for the proportion of units treated by an ITR. The idea is that an efficacious ITR should outperform a non-individualized treatment rule, which randomly assigns the same proportion of units to the treatment condition. We extend the PAPE to the settings with a binding budget constraint. Finally, we propose the area under the prescriptive effect curve (AUPEC) as a univariate summary performance measure of an ITR under a range of budget constraint.

2.1 The Setup

Following the literature, we define an ITR as a deterministic map from the covariate space $\mathcal{X}$ to the binary treatment assignment (e.g., Qian and Murphy, 2011; Zhao et al., 2012),

[TABLE]

Let $T_{i}$ denote the treatment assignment indicator variable, which is equal to $1$ if unit $i$ is assigned to the treatment condition, i.e., $T_{i}\in\mathcal{T}=\{0,1\}$. For each unit, we observe the outcome variable $Y_{i}\in\mathcal{Y}$ as well as the vector of pre-treatment covariates, $\mathbf{X}_{i}\in\mathcal{X}$, where $\mathcal{Y}$ is the support of the outcome variable. We assume no interference between units and denote the potential outcome for unit $i$ under the treatment condition $T_{i}=t$ as $Y_{i}(t)$ for $t=0,1$. Then, the observed outcome is given by $Y_{i}=Y_{i}(T_{i})$.

Assumption 1 (No Interference between Units)

The potential outcomes for unit $i$ do not depend on the treatment status of other units. That is, for all $t_{1},t_{2},\ldots,t_{n}\in\{0,1\}$, we have, $Y_{i}(T_{1}=t_{1},T_{2}=t_{2},\ldots,T_{n}=t_{n})=Y_{i}(T_{i}=t_{i}).$

Under this assumption, the existing literature almost exclusively focuses on the derivation of an optimal ITR that maximizes the following population average value (e.g., Qian and Murphy, 2011; Zhao et al., 2012; Zhou et al., 2017),

[TABLE]

Next, we show that $\lambda_{f}$ may not be the best evaluation metric in some cases.

2.2 The Population Average Prescriptive Effect

We now introduce our main evaluation metric, the Population Average Prescriptive Effect (PAPE). The PAPE is based on two ideas. First, it is reasonable to expect a good ITR to outperform a non-individualized treatment rule, which does not use any information about individual units when deciding who should receive the treatment. Second, a budget constraint should be considered since the treatment is often costly. This means that a good ITR should identify units who benefit from the treatment most. These two considerations lead to the random treatment rule, which assigns, with equal probability, the same proportion of units to the treatment condition, as a natural baseline for comparison.

Figure 1 illustrates the importance of accounting for the proportion of units treated by an ITR. In this figure, the horizontal axis represents the proportion treated and the vertical axis represents the average outcome under an ITR. The example shows that an ITR $g$ (red), which treats 80% of units, has a greater average value than another ITR $f$ (blue), which treats 20% of units, i.e., $\mathbb{E}[Y_{i}(f)]<\mathbb{E}[Y_{i}(g)]$. Despite this fact, $g$ is outperformed by the random treatment rule (black), which treats the same proportion of units, whereas $f$ does a better job than the random treatment rule. This is indicated by the fact that the black solid line is placed above $\mathbb{E}[Y_{i}(g)]$ and below $\mathbb{E}[Y_{i}(f)]$.

To overcome this undesirable property of the average value, we propose an alternative evaluation metric that compares the performance of an ITR with that of the random treatment rule. The random treatment rule serves as a natural baseline because it treats the same proportion of units without any individual information. This is analogous to the predictive setting, in which a classification algorithm is often compared to random classification.

Formally, let $p_{f}=\Pr(f(\mathbf{X}_{i})=1)$ denote the population proportion of units assigned to the treatment condition under ITR $f$. Without loss of generality, we assume a positive average treatment effect $\tau=\mathbb{E}\{Y_{i}(1)-Y_{i}(0)\}>0$ so that the random treatment rule assigns the exactly proportion $p_{f}>0$ of the units to the treatment group. If the treatment is on average harmful (a testable condition using the experimental data), the best random treatment rule is to treat no one. In that case, the estimation of the average value is sufficient for the evaluation. We define the population average prescription effect (PAPE) of ITR $f$ as the following difference in the average value between the ITR and random treatment rule,

[TABLE]

One motivation for the PAPE is that administering a treatment is often expensive. Consider a costly treatment that does not harm anyone but only benefits a relatively small fraction of people. If we do not impose a budget constraint, treating everyone is the best ITR but such a policy does not use any individual-level information. Thus, to further evaluate the efficacy of an ITR, we extend the PAPE to the settings with a budget constraint.

2.3 Incorporating a Budget Constraint

With a budget constraint, we cannot simply treat all units who are predicted to benefit from the treatment. Instead, an ITR must be based on a scoring rule that sorts units according to their treatment priority: a unit with a greater score has a higher priority to receive the treatment. Let $s:\mathcal{X}\longrightarrow\mathcal{S}$ be such a scoring rule where $\mathcal{S}\subset\mathbb{R}$. For simplicity, we assume that the scoring rule is bijective, i.e., $s(\mathbf{x})\neq s(\mathbf{x}^{\prime})$ for any $\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$ with $\mathbf{x}\neq\mathbf{x}^{\prime}$. This assumption is not restrictive as we can always redefine $\mathcal{X}$ such that the assumption holds.

We define an ITR based on a scoring rule by assigning a unit to the treatment group if and only if its score is higher than a threshold, $c$,

[TABLE]

Under a binding budget constraint $p$, we define the threshold that corresponds to the maximal proportion of treated units under the budget constraint, i.e.,

[TABLE]

Our framework allows for any arbitrary scoring rule. A popular scoring rule is the conditional average treatment effect (CATE),

[TABLE]

Researchers have studied the estimation of the CATE using various machine learning algorithms such as tree-based methods and regularized regressions.

We emphasize that the scoring rule need not be based on the CATE. In fact, policy makers rely on various indexes. Such examples include the MELD (Kamath et al., 2001) that determines liver transplant priority, and the IDA Resource Allocation Index that informs the World Bank about the provision of economic aid.

Given this setup, we generalize the PAPE to the setting with a budget constraint. As before, without loss of generality, we assume the treatment is on average beneficial, i.e., $\tau=\mathbb{E}\{Y_{i}(1)-Y_{i}(0)\}>0$, so that the constraint is binding for the random treatment rule treating at most $100\times p\%$ of units. The PAPE with a budget constraint $p$ is defined as,

[TABLE]

A budget constraint facilitates the comparison of multiple ITRs on the same footing. Suppose that we compare two ITRs, $f$ and $g$, using the difference in their average values,

[TABLE]

While this quantity is useful, like the average value, it also fails to take into account the proportion of units assigned to the treatment condition under each ITR.

We can address this issue by comparing the efficacy of two ITRs under the same budget constraint. Formally, we define the Population Average Prescriptive Effect Difference (PAPD) under budget $p$ as,

[TABLE]

2.4 The Area under the Prescriptive Effect Curve

Since the PAPE (Eqn (3)) varies as a function of budget constraint $p$, it would be useful to develop a summary performance metric of an ITR over a range of $p$. We propose the area under the prescriptive effect curve (AUPEC) as a metric analogous to the area under the receiver operating characteristic curve (AUROC) for classification performance.

Figure 2 graphically illustrates the AUPEC. Similar to Figure 1, the vertical and horizontal axes represent the average outcome and the budget, respectively. The budget is operationalized as the maximal proportion treated. The red solid curve corresponds to the average value of an ITR $f$ as a function of budget constraint $p$, i.e., $\mathbb{E}\{Y_{i}(f(\mathbf{X}_{i},c_{p}(f)))\}$, whereas the black solid line represents the average value of the random treatment rule. The AUPEC corresponds to the area under the red curve minus the area under the black line, which is shown as a red shaded area.

Thus, the AUPEC represents the average performance of an ITR relative to the random treatment rule over the entire range of budget constraint (one could also compute the AUPEC over a specific range of budgets). Unlike the previous work (e.g., Rzepakowski and Jaroszewicz, 2012), we do not require an ITR to assign the maximal proportion of units to the treatment condition though such a constraint can also be imposed if desired. For example, treating more than a certain proportion of units will reduce the average outcome if these additional units are harmed by the treatment. This is indicated by the flatness of the red line after $p_{f}$ in Figure 2.

Formally, for a given ITR $f$, we define the AUPEC as,

[TABLE]

where $c^{\ast}$ is the pre-determined minimum score such that one would be treated in the absence of a budget constraint, and $p_{f}=\Pr(f(\mathbf{X}_{i},c^{*})=1)$ denotes the maximal proportion of units assigned to the treatment condition under the ITR with no budget constraint. The last term represents the area under the random treatment rule.

Different values of the threshold $c^{\ast}$ are possible. If the goal is to treat only those who on average benefit from the treatment, then we could use the CATE as a scoring rule and set $c^{\ast}=0$. Another example is the use of the MELD score as the scoring rule and choose an appropriate value of $c^{\ast}$ so that a sufficiently healthy patient is never considered for a transplant. Finally, setting $c^{\ast}=-\infty$ would represent the setting where the maximum possible units should be treated regardless of the scoring rule.

We further note that the AUPEC is a generalization of the QINI coefficient widely utilized in literature for uplift modeling (Radcliffe, 2007). Formally, the population-level QINI coefficient is commonly defined as,

[TABLE]

After some algebra, we can rewrite this quantity in the following form,

[TABLE]

Thus, the QINI coefficient is (up to a constant factor $n$) a special case of AUPEC when $c^{\ast}=-\infty$ and $p_{f}=1$. The choice of $c^{\ast}=-\infty$ may be reasonable in the applications where the treatment can be assumed to be never harmful, i.e., $Y_{i}(1)\geq Y_{i}(0)$. In such cases, under no budget constraint one would treat the entire population.

To enable a comparison of efficacy across different datasets, the AUPEC can be normalized to be scale-invariant by shifting $\Gamma_{f}$ by $\mathbb{E}(Y_{i}(0))$ and dividing by $\tau=\mathbb{E}(Y_{i}(1)-Y_{i}(0))$,

[TABLE]

This normalized AUPEC is invariant to the affine transformation of the outcome variables, $Y_{i}(1),Y_{i}(0)$, while the standard AUPEC is only invariant to a constant shift. The normalized AUPEC takes a value in $[0,1]$, and has an intuitive interpretation as the average percentage outcome gain using the ITR $f$ compared to the random treatment rule, under the uniform prior distribution over the percentage treated.

3 Estimation and Inference

Having introduced the evaluation metrics, we show how to estimate them and compute standard errors under the repeated sampling framework of Neyman (1923). Here, we consider the setting, in which researchers are interested in evaluating the performance of fixed ITRs. In other words, throughout this section, ITR $f$ is assumed to be known and has no estimation uncertainty. For example, one may first construct an ITR based on an existing data set (experimental or observational), and then conduct a new experiment to evaluate its performance. In Section 4, we extend the methodology developed here to the setting, in which the same experimental data set is used to both construct and evaluate an ITR.

For the rest of this paper, we assume that we have a simple random sample of $n$ units from a super-population, $\mathcal{P}$. We conduct a completely randomized experiment, in which $n_{1}$ units are randomly assigned to the treatment condition with probability $n_{1}/n$ and the rest of the $n_{0}(=n-n_{1})$ units are assigned to the control condition. While it is straightforward to allow for unequal treatment assignment probabilities across units, for the sake of simplicity, we assume complete randomization. We formalize these assumptions below.

Assumption 2 (Random Sampling of Units)

Each of $n$ units, represented by a three-tuple consisting of two potential outcomes and pre-treatment covariates, is assumed to be independently sampled from a super-population $\mathcal{P}$, i.e.,

[TABLE]

Assumption 3 (Complete Randomization)

For any $i=1,2,\ldots,n$, the treatment assignment probability is given by,

[TABLE]

where $\sum_{i=1}^{n}T_{i}=n_{1}$.

We now present the results under a binding budget constraint. The results for the average value and PAPE with no budget constraint appear in Appendix A.1.

3.1 The Population Average Prescription Effect (PAPE)

To estimate the PAPE with a binding budget constraint $p$ (Eqn (3)), we consider the following estimator,

[TABLE]

where $\hat{c}_{p}(f)=\inf\{c\in\mathbb{R}:\sum_{i=1}^{n}f(\mathbf{X}_{i},c)\leq np\}$ represents the estimated threshold given the maximal proportion of treated units $p$. Unlike the case of no budget constraint (see Appendix A.1.2), the bias is not zero because the threshold $c_{p}(f)$ needs to be estimated without an assumption about the distribution of the score produced by the scoring rule. We derive an upper bound of bias and the exact variance.

Theorem 1

(Bias Bound and Exact Variance of the PAPE Estimator with a Budget Constraint)* Under Assumptions 1, 2, and 3, the bias of the proposed estimator of the PAPE with a budget constraint $p$ defined in Eqn (7) can be bounded as follows,*

[TABLE]

where any given constant $\epsilon>0$, $B(\epsilon,\alpha,\beta)$ is the incomplete beta function (if $\alpha\leq 0$ and $\beta>0$, we set $B(\epsilon,\alpha,\beta):=H(\epsilon)$ for all $\epsilon$ where $H(\epsilon)$ is the Heaviside step function), and

[TABLE]

The variance of the estimator is given by,

[TABLE]

where $S_{fpt}^{2}=\sum_{i=1}^{n}(Y_{fpi}(t)-\overline{Y_{fp}(t)})^{2}/(n-1)$ and $\kappa_{ft}(p)=\mathbb{E}(\tau_{i}\mid f(\mathbf{X}_{i},\hat{c}_{p}(f))=t)$ with $Y_{fpi}(t)=\left(f(\mathbf{X}_{i},\hat{c}_{p}(f))-p\right)Y_{i}(t)$, and $\overline{Y_{fp}(t)}=\sum_{i=1}^{n}Y_{fpi}(t)/n$, for $t=0,1$.

Proof is given in Appendix A.2. The last term in the variance accounts for the variance due to estimating $c_{p}(f)$. The variance can be consistently estimated by replacing each unknown parameter with its sample analogue, i.e., for $t=0,1$,

[TABLE]

where $Y_{fpi}=(f(\mathbf{X}_{i},\hat{c}_{p}(f))-p)Y_{i}$ and $\overline{Y}_{fpt}=\sum_{i=1}^{n}\mathbf{1}\{T_{i}=t\}Y_{fpi}/n_{t}$. To estimate the term that appears in the denominator of $\gamma_{fp}(\epsilon)$ as part of the upper bound of bias, we may assume that the CATE, $\mathbb{E}(\tau_{i}\mid s(\mathbf{X}_{i})=c)$, is continuous in $c$. Continuity is often assumed when estimating the CATE (e.g., Künzel et al., 2018; Wager and Athey, 2018). We may also utilize an upper bound for CATE, if known, to estimate the bias conservatively.

Building on the above results, we also consider the comparison of two ITRs under the same budget constraint, using the following estimator of the PAPD (Eqn (5)),

[TABLE]

Theorem A3 in Appendix A.3 derives the bias bound and exact variance of this estimator. In one of their applications, Zhou et al. (2018) apply the $t$-test to the cross-validated test statistic similar to the one introduced here under no budget constraint. However, no formal justification of this procedure is given under the cross-validation setting and it cannot be readily extended to the case with a budget constraint. In contrast, our methodology is applicable under these settings as well (see Section 4).

3.2 The Area under the Prescriptive Effect Curve

Next, we consider the estimation and inference about the AUPEC (Eqn (6)). Let $n_{f}$ represent the maximum number of units in the sample that the ITR $f$ would assign under no budget constraint, i.e., $\hat{p}_{f}=n_{f}/n=\sum_{i=1}^{n}f(\mathbf{X}_{i},c^{*})/n$. We propose the following estimator of the AUPEC,

[TABLE]

The following theorem shows a bias bound and the exact variance of this estimator.

Theorem 2 (Bias and Variance of the AUPEC

Estimator)

Under Assumptions 1, 2, and 3, the bias of the AUPEC estimator defined in Eqn (9) can be bounded as follows,

[TABLE]

where any given constant $\epsilon>0$, $B(\epsilon,\alpha,\beta)$ is the incomplete beta function (if $\alpha=0$ and $\beta>0$, we set $B(\epsilon,\alpha,\beta):=H(\epsilon)$ for all $\epsilon$ where $H(\epsilon)$ is the Heaviside step function), and

[TABLE]

The variance is given by,

[TABLE]

where $Z$ is a Binomial random variable with size $n$ and success probability $p_{f}$, and $S^{\ast 2}_{ft}=\sum_{i=1}^{n}(Y^{\ast}_{i}(t)-\overline{Y^{\ast}(t)})^{2}/(n-1)$, $\kappa_{ft}(k/n)\ =\ \mathbb{E}(Y_{i}(1)-Y_{i}(0)\mid f(\mathbf{X}_{i},\hat{c}_{k/n}(f))=t)$, with $Y_{i}^{\ast}(t)=\left[\left\{\sum_{z=1}^{n_{f}}f(\mathbf{X}_{i},\hat{c}_{z/n}(f))+(n-n_{f})f(\mathbf{X}_{i},\hat{c}_{\hat{p}_{f}}(f))\right\}/n-\frac{1}{2}\right]Y_{i}(t)$ and $\overline{Y^{\ast}(t)}=\sum_{i=1}^{n}Y_{i}^{\ast}(t)/n$, for $t=0,1$.

Proof is given in in Appendix A.4. When $c^{\ast}=-\infty$ (i.e., the AUPEC equals the QINI coefficient), the estimator is unbiased, implying that the bias comes from estimating the proportion treated $p_{f}$ under no budget constraints and $c^{\ast}>-\infty$. As before, $\mathbb{E}({S_{ft}^{\ast}}^{2})$ does not equal $\mathbb{V}(Y_{i}^{\ast}(t))$ due to the need to estimate the terms $c_{z/n}$ for all $z$, and the additional terms account for the variance of estimation. We can consistently estimate the upper bound of bias, for example, under by assuming that the CATE is continuous. We may also utilize an upper bound for CATE, if known, to estimate the bias conservatively.

To estimate the variance, we replace each unknown parameter with its sample analogue,

[TABLE]

for $t=0,1$. In the extreme cases with $z\to 1$ for $t=1$ and $z\to n$ for $t=0$, each denominator in Eqn (LABEL:eq:kappa2est) is likely to be close to zero. In such cases, we instead use the estimator $\widehat{\kappa_{f1}(z_{\min}/n)}$ for all $z<z_{\min}$ where $z_{\min}$ is the smallest $z$ such that Eqn (LABEL:eq:kappa2est) for $\kappa_{f1}(z/n)$ does not lead to division by zero. Similarly, for $t=0$, we use $\widehat{\kappa_{f0}(z_{\max}/n)}$ for all $z>z_{\max}$ where $z_{\max}$ is the largest $z$.

For the terms involving the binomial random variable $Z$, we first note that, when fully expanded out, they are the polynomials of $p_{f}=\mathbb{E}(f(\mathbf{X}_{i}))$. To estimate the polynomials, we can utilize their unbiased estimators as discussed in Stuard and Ord (1994), i.e., $\hat{p}_{f}^{z}=s(s-1)\cdots(s-z+1)/\{n(n-1)\cdots(n-z+1)\}\}$ where $s=\sum_{i=1}^{n}f(\mathbf{X}_{i})$ is unbiased for $p_{f}^{z}$ for all $z\leq n$. When the sample size is large, this estimation method is computationally inefficient and unstable due to the presence of high powers. Hence, we may use the Monte Carlo sampling of $Z$ from a Binomial distribution with size $n$ and success probability $\hat{p}_{f}$. In our simulation study, we show that this Monte Carlo approach is effective even when the sample size is small (see Section 5).

Finally, a consistent estimator for the normalized AUPEC is given by,

[TABLE]

The variance of $\widehat{\Gamma}^{\ast}_{f}(\mathcal{Z})$ can be estimated using the Taylor expansion of quotients of random variables to an appropriate order as detailed in Stuard and Ord (1994).

4 Estimating and Evaluating ITRs Using the Same Experimental Data

We next consider a common situation, in which researchers use the same experimental data to both estimate and evaluate an ITR via cross-validation. This differs from the setting we have analyzed so far, in which a fixed ITR is given for evaluation. We first extend the evaluation metrics introduced in Section 2 to the current setting with estimated ITRs. We then develop inferential methods under the Neyman’s repeated sampling framework by accounting for both estimation and evaluation uncertainties. Below, we consider the scenario, in which researchers face a binding budget constraint. Appendix A.5 presents the results for the case with no budget constraint.

4.1 Evaluation Metrics

Suppose that we have the data from a completely randomized experiment as described in Section 3. We first estimate an ITR $f$ by applying a machine learning algorithm $F$ to training data $\mathcal{Z}^{tr}$. Then, under a budget constraint of the maximal proportion of treated units $p$, we use test data to evaluate the resulting estimated ITR $\hat{f}_{\mathcal{Z}^{tr}}$. As before, we assume that this constraint is binding, i.e., $p<p_{F}$ where $p_{F}=\Pr(\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i})=1)$ represents the proportion of treated units under the ITR without a budget constraint.

Formally, an machine learning algorithm $F$ is a deterministic map from the space of training data $\mathcal{Z}$ to that of scoring rules $\mathcal{S}$,

[TABLE]

Then, for a given training data set $\mathcal{Z}^{tr}$, the estimated ITR is given by,

[TABLE]

where $\hat{s}_{\mathcal{Z}^{tr}}=F(\mathcal{Z}^{tr})$ is the estimated scoring rule and $c_{p}(\hat{f}_{\mathcal{Z}^{tr}})\ =\ \inf\{c\in\mathbb{R}:\Pr(\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},c)=1\mid\mathcal{Z}^{tr})\leq p\}$ is the threshold based on the maximal proportion of treated units $p$. The CATE is a natural choice for the scoring rule, i.e., $s_{\mathcal{Z}^{tr}}(\mathbf{X}_{i})=\mathbb{E}(\tau_{i}\mid\mathbf{X}_{i})$. We need not assume that $\hat{s}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i})$ is consistent for the CATE.

To extend the PAPE (Eqn (3)), we first define the population proportion of units with $\mathbf{X}_{i}=\mathbf{x}$ who are assigned to the treatment condition under the estimated ITR as,

[TABLE]

While $c_{p}(\hat{f}_{\mathcal{Z}^{tr}})$ depends on the specific ITR generated from the training data $\mathcal{Z}^{tr}$, the population proportion of treated units averaged over the sampling of training data, $\bar{f}_{Fp}(\mathbf{X}_{i})$, only depends on $p$.

Lastly, the PAPE of the estimated ITR under budget constraint $p$ is defined as,

[TABLE]

This evaluation metric corresponds to neither that of a specific ITR estimated from the whole experimental data set nor its expectation. Rather, we are evaluating the efficacy of a learning algorithm that is used to estimate an ITR using the same experimental data.

We can also compare estimated ITRs by further generalizing the definition of the PAPD (Eqn (5)) to the current setting. Specifically, we define the PAPD between two machine learning algorithms, $F$ and $G$, under budget constraint $p$ as,

[TABLE]

Finally, we consider the AUPEC of an estimated ITR. Specifically, the AUPEC of an machine learning algorithm $F$ is defined as,

[TABLE]

where $p_{\hat{f}}=\Pr(\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i})=1)$ is the maximal population proportion of units treated by the estimated ITR $\hat{f}_{\mathcal{Z}^{tr}}$.

4.2 Estimation and Inference

Rather than simply splitting the data into training and test sets (in such a case, the inferential procedure for fixed ITRs is applicable), we maximize efficiency by using cross-validation to estimate the evaluation metrics introduced above. First, we randomly split the data into $K$ subsamples of equal size $m=n/K$ by assuming, for the sake of notational simplicity, that $n$ is a multiple of $K$. Then, for each $k=1,2,\ldots,K$, we use the $k$th subsample as a test set $\mathcal{Z}_{k}=\{\mathbf{X}_{i}^{(k)},T_{i}^{(k)},Y_{i}^{(k)}\}_{i=1}^{m}$ with the data from all $(K-1)$ remaining subsamples as the training set $\mathcal{Z}_{-k}=\{\mathbf{X}_{i}^{(-k)},T_{i}^{(-k)},Y_{i}^{(-k)}\}_{i=1}^{n-m}$.

To simplify notation, we assume that the number of treated (control) units is identical across different folds and denote it as $m_{1}$ ($m_{0}=m-m_{1}$). For each split $k$, we estimate an ITR by applying a learning algorithm $F$ to the training data $\mathcal{Z}_{-k}$,

[TABLE]

We then evaluate the performance of the learning algorithm $F$ by computing an evaluation metric of interest $\tau$ based on the test data $\mathcal{Z}_{k}$. Repeating this $K$ times for each $k$ and averaging the results gives a cross-validation estimator of the evaluation metric. Algorithm 1 formally presents this estimation procedure. The variance formula for the estimated evaluation metric is omitted as it is generally a complex function $v(\cdot)$ (see Theorem 3 below).

We develop the inferential methodology for the evaluation based on the cross-validation procedure described above under the Neyman’s repeated sampling framework. We focus on the case with a binding budget constraint. The results with no budget constraint appear in Appendix A.5. We begin by introducing the cross-validation estimator of the PAPE with a binding budget constraint $p$,

[TABLE]

where $\hat{\tau}_{fp}$ is defined in Eqn (7).

Like the fixed ITR case, the bias of the proposed estimator is not exactly zero. However, we are able to show that the bias can be upper bounded by a small quantity while the exact randomization variance can still be derived.

Theorem 3

(Bias Bound and Exact Variance of the Cross-Validation PAPE Estimator with a Budget Constraint)* Under Assumptions 1, 2, and 3, the bias of the cross-validation PAPE estimator with a budget constraint $p$ defined in Eqn (15) can be bounded as follows,*

[TABLE]

where any given constant $\epsilon>0$, $B(\epsilon,\alpha,\beta)$ is the incomplete beta function (if $\alpha=0$ and $\beta>0$, we set $B(\epsilon,\alpha,\beta):=H(\epsilon)$ for all $\epsilon$ where $H(\epsilon)$ is the Heaviside step function), and

[TABLE]

The variance of the estimator is given by,

[TABLE]

where $S_{\hat{f}pt}^{2}=\sum_{i=1}^{m}(Y_{\hat{f}pi}(t)-\overline{Y_{\hat{f}p}(t)})^{2}/(m-1)$, and $S_{Fp}^{2}=\sum_{k=1}^{K}(\hat{\tau}_{\hat{f}_{-k},p}(\mathcal{Z}_{k})-\overline{\hat{\tau}_{\hat{f}_{-k},p}(\mathcal{Z}_{k})})^{2}/(K-1)$, and $\kappa_{Ft}(p)=\mathbb{E}(\tau_{i}\mid\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},\hat{c}_{p}(\hat{f}_{\mathcal{Z}^{tr}}))=t)$, with $Y_{\hat{f}pi}(t)=\{\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},\hat{c}_{p}(\hat{f}_{\mathcal{Z}^{tr}}))-p\}Y_{i}(t)$, $\overline{Y_{\hat{f}p}(t)}=\sum_{i=1}^{n}Y_{\hat{f}pi}(t)/n$, and $\overline{\hat{\tau}_{\hat{f}_{-k},p}(\mathcal{Z}_{k})}=\sum_{k=1}^{K}\hat{\tau}_{\hat{f}_{-k},p}(\mathcal{Z}_{k})/K$, for $t=0,1$.

Proof is given in Appendix A.6. The estimation of the term $\mathbb{E}(\widetilde{S}_{\hat{f}t}^{2})$ is done similarly as before. For $\kappa_{Ft}(p)$, we replace it with its sample analogue:

[TABLE]

To estimate the term that appears in the denominator of $\gamma_{p}(\epsilon)$ as part of the upper bound of bias, we assume that the CATE, i.e., $\mathbb{E}(Y_{i}(1)-Y_{i}(0)\mid\hat{s}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i})=c)$, is continuous in $c$, and replace the maximum with a point estimate. We may also utilize an upper bound for CATE, if known, to estimate the bias conservatively. Building on this result, Appendix A.7 shows how to compare two estimated ITRs by estimating the PAPD (Eqn (12)).

Finally, we consider the following cross-validation estimator of the AUPEC for an estimated ITR (Eqn (13)),

[TABLE]

where $\widehat{\Gamma}_{f}$ is defined in Eqn (9). This can be seen as an overall statistic that measures the prescriptive performance of the machine learning algorithm $F$ on the dataset under cross-validation. The following theorem derives a bias bound and the exact variance of this cross-validation estimator.

Theorem 4

(Bias Bound and Exact Variance of the Cross-validation AUPEC Estimator)* Under Assumptions 1, 2, and 3, the bias of the AUPEC estimator defined in Eqn (17) can be bounded as follows,*

[TABLE]

where any given constant $\epsilon>0$, $B(\epsilon,\alpha,\beta)$ is the incomplete beta function (if $\alpha=0$ and $\beta>0$, we set $B(\epsilon,\alpha,\beta):=H(\epsilon)$ for all $\epsilon$ where $H(\epsilon)$ is the Heaviside step function), and

[TABLE]

The variance is given by,

[TABLE]

where $Z$ is a Binomial random variable with size $m$ and success probability $p_{\hat{f}}$, $S^{\ast 2}_{\hat{f}t}=\sum_{i=1}^{m}(Y^{\ast}_{\hat{f}i}(t)-\overline{Y_{\hat{f}}^{\ast}(t)})^{2}/(m-1)$, $S_{F}^{2}=\sum_{k=1}^{K}(\widehat{\Gamma}_{\hat{f}_{-k}}(\mathcal{Z}_{k})-\overline{\widehat{\Gamma}_{\hat{f}_{-k}}(\mathcal{Z}_{k})})^{2}/(K-1)$, and $\kappa_{Ft}(z/m)=\mathbb{E}(\tau_{i}\mid\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},\hat{c}_{z/m}(\hat{f}_{\mathcal{Z}^{tr}})=t)$ with $\overline{Y_{\hat{f}}^{\ast}(t)}=\sum_{i=1}^{m}Y_{\hat{f}i}^{\ast}(t)/m$, $\overline{\widehat{\Gamma}_{\hat{f}_{-k}}(\mathcal{Z}_{k})}=\sum_{k=1}^{K}\widehat{\Gamma}_{\hat{f}_{-k}}(\mathcal{Z}_{k})/K$, and $Y_{\hat{f}i}^{\ast}(t)=\left[\left\{\sum_{z=1}^{m_{f}}\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},\hat{c}_{z/m}(\hat{f}_{\mathcal{Z}^{tr}}))+(m-m_{f})\hat{f}_{\mathcal{Z}^{tr}}(\mathbf{X}_{i},\hat{c}_{z/m}(\hat{f}_{\mathcal{Z}^{tr}}))\right\}/m-\frac{1}{2}\right]Y_{i}(t)$ for $t=0,1$.

Proof is similar to that of Theorem 2. The estimation of $\mathbb{E}(S_{\hat{f}1}^{\ast 2})$, $\mathbb{E}(S_{\hat{f}0}^{\ast 2})$, and $\mathbb{E}(S_{F}^{\ast 2})$ is the same as before, and the $\kappa_{Ft}(p)$ term can be estimated using Eqn (16).

5 A Simulation Study

We conduct a simulation study to examine the finite sample performance of our methodology, for both fixed and estimated ITRs. We find that the empirical coverage probability of the confidence interval, based on the proposed variance, approximates its nominal rate even in a small sample. We also find that the bias is minimal even when the proposed estimator is not unbiased and that our variance bounds are tight.

5.1 Data Generation Process

Our data generating process (DGP) is based on the one used in the 2017 Atlantic Causal Inference Conference (ACIC) Data Analysis Challenge (Hahn et al., 2018). A total of 8 covariates $\mathbf{X}$ are taken from the Infant Health and Development Program, which originally had 58 covariates and $4,302$ observations. In our simulation, the population distribution of covariates is assumed to equal the empirical distribution of this data set. Therefore, we obtain each simulation sample via boostrap. We vary the sample size: $n=100,500,2000$.

We use the same outcome model as the one used in the competition,

[TABLE]

where $\pi(\mathbf{x})=1/[1+\exp\{3(x_{1}+x_{43}+0.3(x_{10}-1))-1\}]$, $\mu(\mathbf{x})=-\sin(\Phi(\pi(\mathbf{x})))+x_{43}$, and $\tau(\mathbf{x})=\xi(x_{3}x_{24}+(x_{14}-1)-(x_{15}-1))$ with $\Phi(\cdot)$ representing the standard Normal CDF and $x_{j}$ indicating a specific covariate in the data set. One important difference is that we assume a complete randomized experiment whereas the original DGP generated the treatment using a function of covariates to emulate an observational study. As in the competition, we focus on two scenarios regarding the treatment effect size by setting $\xi$ equal to $2$ (“high”) and $1/3$ (“low”). Although the original DGP included four different error distributions, we use the i.i.d. error, $\sigma(\mathbf{X}_{i})\epsilon_{i}$ where $\sigma(\mathbf{x})=0.25\sqrt{\mathbb{V}(\mu(\mathbf{x})+\pi(\mathbf{x})\tau(\mathbf{x}))}$ and $\epsilon_{i}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}\mathcal{N}(0,1)$.

For fixed ITRs, we can directly compute the true values of our causal quantities of interest using the outcome model specified in Eqn (18) and evaluate each quantity based on the entire original data set. This computation is valid because we assume the population distribution of covariates is equal to the empirical distribution of the original data set. For the estimated ITR case, however, we do not have an analytical expression for the true value of a causal quantity of interest. Therefore, we obtain an approximate true value via Monte Carlo simulation. We generate 10,000 independent datasets based on the same DGP, and train the specified algorithm $F$ on each of the datasets using 5-fold cross-validation (i.e., $K=5$). Then, we use the sample mean of our estimated causal quantity across 10,000 data sets as our approximate truth.

We evaluate Bayesian Additive Regression Trees (BART) (Chipman et al., 2010; Hahn et al., 2020), which had the best overall performance in the original competition. We compare this model with two other popular methods: Causal Forest (Athey et al., 2019) as well as the LASSO, which includes all main effects and two-way interaction effects between the treatment and all covariates (Tibshirani, 1996). All three models are trained on the original data from the 2017 ACIC Data Challenge. The number of trees was tuned through the 5-fold cross validation for BART and Causal Forest. The regularization parameter was tuned similarly for LASSO. All models were cross-validated on the PAPE. For implementation, we use R 3.4.2 with bartMachine (version 1.4.2) for BART, grf (version 0.10.2) for Causal Forest, and glmnet (version 2.0.13) for LASSO. Once the models are trained, an ITR is derived based on the magnitude of the estimated CATE $\hat{\tau}(\mathbf{x})$, i.e., $f(\mathbf{X}_{i})=\mathbf{1}\{\hat{\tau}(\mathbf{X}_{i})>0\}$.

5.2 Results

We first present the results for fixed ITRs followed by those for estimated ITRs. Table 1 presents the bias and standard deviation of each estimator for fixed ITRs as well as the coverage probability of its 95% confidence intervals based on 1,000 Monte Carlo trials. The results are shown separately for the high and low treatment effects scenarios. We estimate the PAPE $\tau_{f}$ for BART without a budget constraint as well as the PAPE with a budget constraint of 20% as the maximal proportion of treated units, $\tau_{f}(c_{0.2})$. In addition, we estimate the AUPEC $\Gamma_{f}$ and compute the difference in the PAPE or PAPD between BART and Causal Forest ($\Delta(f,g)$), and between BART and LASSO ($\Delta(f,h)$).

Under both scenarios and across sample sizes, the bias of our estimator is small. Moreover, the coverage rate of 95% confidence intervals is close to their nominal rate even when the sample size is small. Although we can only bound the variance when estimating the PAPD between two ITRs (i.e., $\Delta_{0.2}(f,g)$ and $\Delta_{0.2}(f,h)$), the coverage stay close to 95%, implying that the bound for covariance has little effect on the variance estimation.

For estimated ITRs, Table 2 presents the results of LASSO under cross-validation. For BART and Causal Forest, obtaining an accurate Monte Carlo estimate of the true causal parameter values under cross-validation takes a prohibitively large amount of time. While the out-of-bag estimates of such true values can be computed, they have been shown to create bias under certain scenarios (Janitza and Hornung, 2018). These true values are generally greater for a larger sample size because LASSO performs better with more data.

The proposed cross-validated estimators are approximately unbiased even for $n=100$. The coverage is generally around or above the nominal $95\%$ value, reflecting the conservative estimate of the variance. For the PAPE without and with budget constraint, i.e., $\hat{\tau}_{F}$ and $\hat{\tau}_{F}(c_{0.2})$, the coverage is close to the nominal value. This indicates that the bias of the proposed conservative variance estimator is relatively small even though the number of folds for cross-validation is only $K=5$. The performance of the proposed methodology is good even when the sample size is as small as $n=100$. When the sample size is $n=500$, the standard deviation of the cross-validated estimator is roughly half of the corresponding $n=100$ fixed ITR estimator (this is a good comparison because each fold has 100 observations). This confirms the theoretical efficiency gain that results from cross-validation.

6 An Empirical Application

We apply the proposed methodology to the data from the Tennessee’s Student/Teacher Achievement Ratio (STAR) project, which was a longitudinal study experimentally evaluating the impacts of class size in early education on various outcomes (Mosteller, 1995). Another application based on a canvassing experiment is shown in Appendix A.8.

6.1 Data and Setup

The STAR project randomly assigned over 7,000 students across 79 schools to three different groups: small class,regular class,and regular class with a full-time teacher’s aid. The experiment began when students entered kindergarden and continued through third grade. To create a binary treatment, we focus on the first two groups: small class and regular class without an aid. The treatment effect heterogeneity is important because reducing class size is costly, requiring additional teachers and classrooms. Policy makers who face a budget constraint may be interested in finding out which groups of students benefit most from a small class size so that the priority can be given to those students.

We follow the analysis strategies of the previous studies (e.g., Ding and Lehrer, 2011; McKee et al., 2015) that estimated the heterogeneous effects of small classes on educational attainment. These authors adjust for school-level and student-level characteristics, but do not consider within-classroom interactions. Unfortunately, addressing this limitation is beyond the scope of this paper. We use a total of 10 pre-treatment covariates $\mathbf{X}_{i}$ that include four demographic characteristics of students (gender, race, birth month, birth year) and six school characteristics (urban/rural, enrollment size, grade range, number of students on free lunch, number of students on school buses, and percentage of white students). Our treatment variable is the class size to which they were assigned at kindergarten: small class $T_{i}=1$ and regular class without an aid $T_{i}=0$. For the outcome variables $Y_{i}$, we use three standardized test scores measured at third grade: math, reading, and writing SAT scores.

The resulting data set has a total of 1,911 observations. The estimated average treatment effects (ATE) based on the entire data set are $6.78$ (s.e. = $1.71$), $5.78$ (s.e. = $1.80$), and $3.65$ (s.e. = $1.63$), for the reading, math, and writing scores, respectively. For the fixed test data, the estimated ATEs are similar; $5.10$ (s.e. = $3.07$), $2.78$ (s.e.= $3.15$), and $1.48$ (s.e. = $2.96$).

We evaluate the performance of ITRs using two settings considered above. First, we randomly split the data into the training data (70%) and test data (30%). We estimate an ITR from the training data and then evaluate it as a fixed ITR using the test data. This follows the setup considered in Sections 2 and 3. Second, we consider the evaluation of estimated ITRs based on the same experimental data. We utilize Algorithm 1 with 5-fold cross-validation (i.e., $K=5$). For both settings, we use the same three machine learning algorithms. For Causal Forest, we set tune.parameters = TRUE. For BART, tuning was done on the number of trees. For LASSO, we tuned the regularization parameter while including all interaction terms between covariates and the treatment variable. All tuning was done through the 5-fold cross validation procedure on the training set using the PAPE as the evaluation metric. We then create an ITR as $\mathbf{1}\{\hat{\tau}(\mathbf{x})>0\}$ where $\hat{\tau}(\mathbf{x})$ is the estimated CATE obtained from each fitted model. As mentioned in Section 3, we center the outcome variable $Y$ in evaluating the metrics to minimize the variance of the estimators.

6.2 Results