Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
112,242 characters · 21 sections · 59 citation commands
Efficient Discovery of Heterogeneous Quantile Treatment Effects in Randomized Experiments via Anomalous Pattern Detection
The randomized experiment is employed across many empirical disciplines as an important tool for discovery, by estimating the causal impact of a particular stimulus, treatment or intervention. Moreover, the increasing popularity of large-scale experiments kohavi-online_experiments-2013 has resulted in a widespread interest in discovering fine-grained truths about experimental units, most prominently in the form of heterogeneous treatment effects (HTE). Discovering heterogeneity can be challenging because there are exponentially many subpopulations--with respect to the number of observable covariates--to consider, potentially resulting in multiple hypothesis testing issues and raising questions of unprincipled post-hoc investigation: searching for a fortuitously statistically significant result assmann-subgroup-2000, weisberg-subgroup-2015. Nevertheless, uncovering affected subpopulations can lead to important scientific progress. In a “step toward a new frontier of personalized medicine” saul-bidil-2005, the FDA approved the first race-specific drug, whose impact on African-American subjects was first discovered post-hoc from more general experiments cohn-bidil-1986, cohn-bidil-1991. Conversely, the Perry preschool experiment found significant effects of preschool education on educational and life outcomes barnett-perry-1985,schweinhart-perry-1993,angrist-mhe-2008, while a re-analysis focused on heterogeneity and multiple hypothesis testing concluded that only girls experience these benefits anderson-perry-2008. The original Perry preschool results were fundamental to the creation of the Head Start preschool program angrist-mhe-2008 a national social program that provides, among other services, early childhood education to low-income children. If large-scale medical and policy decisions are made as a result of such experiments, then it is clear that identifying whether there is heterogeneity in treatment effects should be an integral component of the analysis.
In this work we propose a novel computationally efficient framework--Treatment Effect Subset Scanning (TESS)--for discovering which subpopulations in a randomized experiment are the most significantly affected by a treatment. The contributions of this work can be summarized as follows:
These contributions are enabled by structuring the question of causal inference as one of anomalous pattern detection and effect maximization, rather than model fitting and risk minimization. In some contexts, the standard approach of learning an overall model of the treatment effect response surface is desirable; however, in many cases, the identification of affected subpopulations is the primary goal and model learning is simply a step toward this goal. For these cases it seems prudent and efficient to circumvent this first step and solve the subpopulation identification problem by framing it as one of pattern or subset discovery. Such a framing has not previously been considered in the literature.
Most contemporary causal methods are estimators for the (conditional) average treatment effects, or CATE, $\tau_{CATE}(x) = \mathbb{E}\left[Y(1)-Y(0) | X = x \right]$, which in turn limits empirical studies of treatment effects from considering effects beyond mean shifts abadie-qte-2002. However, social scientists argue that effects can greatly vary along the outcome distribution, and distributional impacts beyond the average effect are critical for policy-makers, across a wide range of social programs firpo-qte-2007,abadie-qte-2002,chernozhukov-ivqte-2005,schiele-qte-2016. The primary distributional alternative to ATEs has been Quantile Treatment Effects (QTE) and the subsequent conditional QTEs, or CQTE firpo-qte-2007, koenker-QTE-1978, koenker-hqte-2010, chernozhukov-QTE-generalization-2013 at a given quantile $\alpha$:
where $\mathbb{F}_{Y(1)|X=x}^{-1}$ and $\mathbb{F}_{Y(0)|X=x}^{-1}$ denote the inverse cumulative distribution functions of the outcome $Y$, conditional on covariates $X=x$, under the counterfactual assignments to the treatment group ($W=1$) and control group ($W=0$) respectively.
In this work, we consider the challenge of heterogeneous quantile treatment effects, i.e., detecting the existence of a subpopulation $S$ (characterized by a subset of values for each attribute) for which the CQTE is non-zero, at some quantile, even if there is not a significant effect in the overall population. This motivates the need for a measurement of the heterogeneous treatment effect $\tau_{\text{CQTE}_\alpha}(S)$ and a corresponding test statistic $F_\alpha(S)$ that can be optimized over both subpopulations $S$ and quantiles $\alpha$, capturing unknown heterogeneity in both covariates and the treatment effect distribution, respectively. While a simple extension to (ref),
may appear to be an attractive alternative, this formulation is inadequate for detecting subpopulations with distributional effects. Note that in (ref) the effect is represented by a difference in scalar summaries of the potential outcome distributions, instead of capturing full distributional effect chernozhukov-QTE-generalization-2013,van-structural-2014. Moreover, it first aggregates $\mathbb{F}_{Y(W)|X=x}~\forall x\in S$ to construct $\mathbb{F}_{Y(W)|X \in S}$ and then compares these aggregate conditional distributions, instead of first comparing each $\mathbb{F}_{Y(1)|X=x}$ to the corresponding $\mathbb{F}_{Y(0)|X = x} ~\forall x\in S$ and then aggregating.\footnote{When the effects of interest are simply differences in scalar summaries of each potential outcome distribution, the results are equivalent for any order of comparison and aggregation. However, for more general distributional effects, this equivalence does not hold.} The effect of interest can easily be obfuscated when aggregating before comparing: consider that $Y(W)|X=x$ need not be on the same scale for different $x\in S$. Additionally, (ref) tends to be maximized at extreme values $\alpha\approx 0$ and $\alpha\approx 1$ and is thus highly sensitive to outliers, losing power to detect QTEs occurring at non-extreme values of $\alpha$. Finally, (ref) fails to appropriately calibrate the treatment effect across potential subpopulations $S$ of varying sizes, and therefore in (ref) the optimal $S$ equates to $\max_{x,\alpha} \tau_{\text{CQTE}_{\alpha}}(x)$, i.e., a singular covariate profile (see proof in Appendix (ref)). This last issue also arises when maximizing other popular conditional treatment effect estimands in the literature, such as CATE, over subpopulations.
To avoid these limitations, we first recognize that our primary goals are to test the null hypothesis
for $\alpha \in (0,1)$ and $\forall x$, and to detect subpopulations $S$ for which the two counterfactual outcome distributions differ significantly. An equivalent test is for
Moreover, re-defining $\tau_{\text{CQTE}_{\alpha}}(x) = \mathbb{F}_{Y(1)|X=x}(\mathbb{F}_{Y(0)|X=x}^{-1}(\alpha))$ captures the full distributional effect (as demonstrated by chernozhukov-QTE-generalization-2013,van-structural-2014) prior to aggregation, and is constrained to $(0,1)$, allowing coherent aggregation and calibration over $x\in S$. More precisely, we can define $\tau_{\text{CQTE}_\alpha}(S) = \sum_{x \in S} \tau_{\text{CQTE}_\alpha}(x) P(X=x \mid X \in S)$, and then define a test statistic $F_\alpha(S)$ to measure the significance of the divergence between $\tau_{\text{CQTE}_\alpha}(S)$ and $\alpha$.
Therefore, we can make specific, simple, and testable assumptions about the relationships between each $\mathbb{F}_{Y(0)|X=x}$ and $\mathbb{F}_{Y(1)|X=x}$, under the null and alternative hypotheses, and construct a generalized likelihood ratio test that maximizes detection power for distinguishing these hypotheses:
We define $H_1(S)$ as in (ref) (with constant $\beta_x = \beta$) because of our interest in detecting subsets $S$ where $\mathbb{F}_{Y(1)|X=x}$ differs systematically from $\mathbb{F}_{Y(0)|X=x}$ for $x \in S$, thus grouping together covariate profiles that exhibit similar treatment effects, rather than massively overfitting to individual covariate profiles.\footnote{This is analogous to tree-based methods which assume the same conditional average treatment effect (CATE) for all cells assigned to a given leaf, but rather than looking for a mean shift, we measure how much of the probability density of $Y(1)$ has been “shifted” into the $\alpha$-tail of $Y(0)$.} Moreover, if $\mathbb{F}_{Y(1)|X=x} \ne \mathbb{F}_{Y(0)|X=x}$, we know that there exists some $\alpha$ and some subset $S$ such that $\tau_{\text{CQTE}_\alpha}(S) = \beta \ne \alpha$, while if $\mathbb{F}_{Y(1)|X=x} = \mathbb{F}_{Y(0)|X=x}$, then $\beta = \alpha$ everywhere and there is no shift. Therefore, defining the alternative in this way allows us to prove desirable theoretical results both for detection power and for subset correctness, as described in Theorems (ref)-(ref) and (ref)-(ref) respectively.
As we show in Appendix (ref), the log-likelihood ratio statistic for (ref), for a given sample, corresponds to the Berk-Jones nonparametric scan statistic mcfowland-fgss-2013,berk-bj-1979:
where each $\hat{\tau}_{\text{CQTE}_{\alpha}}(x)$ is computed using its potential outcome empirical distribution, $N(S)$ is the number of treatment group units, and $KL(\beta, \alpha)=\beta \log \frac{\beta}{\alpha} + (1-\beta) \log \frac{1-\beta}{1-\alpha}$ is the Kullback-Leibler divergence between Bernoulli distributions with the corresponding parameters. The statistic in (ref) includes a maximization over subpopulations $S$ and thresholds $\alpha$ to identify the most significantly affected (highest scoring) subpopulation, and its significance can then be determined by a randomization test, appropriately controlling for multiple testing.
Treatment Effect Subset Scan (TESS) is a novel framework for identifying subpopulations in a randomized experiment which experience treatment effects, built atop the heterogeneous quantile treatment effect test statistic established in (ref). Unlike previous methods, TESS structures the challenge of treatment effect identification as an anomalous pattern detection problem--where the objective is to identify patterns of systematic deviations away from expectation--which is then solved by scanning over subpopulations. TESS therefore searches for subsets of values of each attribute for which the distributions of outcomes in the treatment groups are systematically anomalous, i.e., significantly different from their expectation as derived from the control group. More precisely, we define a real-valued outcome of interest $Y$ and a set of discrete covariates $X = (X^{1}, \ldots, X^{d}$), where each $X^{j}$ can take on a vector of values $V^{j}=\{v^{j}_{m}\}_{m = 1...|V^j|}$. We note that continuous covariates can be discretized into categories, using the observed covariate distribution or domain knowledge.\footnote{An extension could include considering the intervals of the continuous covariate created by each of its unique split points (realized values) in the data; this is similar to how tree-based methods determine discrete splits on continuous variables.} With continuous and discrete covariates, the distribution and quantile functions are well-defined and unique for all levels $\alpha\in (0, 1)$ when the outcome $Y$ is real-valued. We define the arity of covariate $X^{j}$ as $|V^{j}|$ (i.e., the cardinality of $V^j$) and note that for any covariate profile $x$ (i.e., a realization of $X$) it follows that $x \in V^1 \times \ldots \times V^d$. We then define a dataset as a sample $\mathcal{N}$ composed of $n$ records (units) $\{R_{1}, \ldots, R_{n}\}$, drawn independently and identically distributed from population $\mathcal{P}$. Each 3-tuple $R_{i} = (Y^{\text{obs}}_i, X_i, W_i)$ is described by an observed potential outcome $Y^{\text{obs}}_i = Y_i(W_i)$, covariates $X_i$, and an indicator variable $W_i$, which indicates if the unit was randomly assigned to the treatment condition; see Table (ref) for a demonstrative example. We define the subpopulations $S$ under consideration to be $S = v^1 \times \ldots \times v^d$, where $v^j \subseteq V^j$. Therefore, we consider subsets $S$ representing subspaces of the attribute space, i.e., the Cartesian product of a subset of values for each attribute. This is important because the treatment of interest may affect multiple values, e.g., African-Americans or Hispanics who live in New York or Pennsylvania. Finally, we wish to find the most anomalous subset
where $F(S)$ is commonly referred to in the anomalous pattern detection literature as a score function, to measure the anomalousness of a subset $S$. In the context of TESS, this function is a test statistic of the treatment effect--i.e., the divergence between the treatment and control group--in subpopulation $S$, like the one defined in (ref).
We accomplish this by first partitioning the experimental dataset into control and treatment groups, and passing the groups to the TESS algorithm. For each unique covariate profile $x$ in the treatment group, TESS uses the control group to compute a conditional outcome distribution $\hat{\mathbb{F}}_{Y^C \mid X=x}$, providing an estimate of the conditional outcome distribution under the null hypothesis $H_{0}$ that the treatment has no effect on units with this profile. Then for each record $R_{i}$ in the treatment group, TESS computes an empirical $p$-value $\hat{p}_{i}$, which serves as a measure of how uncommon it is to see an outcome as extreme as $Y^{\text{obs}}_i$ given $X=x_{i}$ under $H_0$. The ultimate goal of TESS is to discover subpopulations $S$ with a large amount of evidence against $H_0$, i.e., the outcomes of units in $S$ are consistently extreme given $H_0$. Thus, TESS searches for subpopulations which contain an unexpectedly large number of low (significant) empirical $p$-values, as such a subpopulation is more likely to have been affected by the treatment.
After partitioning the data into treatment and control groups, the TESS framework obtains an estimate of the reference distribution for each unique covariate profile in the treatment group. To obtain the estimates of $\mathbb{F}_{Y(0) \mid X}~\forall X$, TESS relies on two assumptions: randomization and a sharp null hypothesis of no treatment effect. First, randomization implies that the potential outcomes $ Y_i(0), Y_i(1) \protect\mathpalette{\protect\independenT}{\perp} W_i~\forall R_i$: selection into treatment and control groups is completely random. Secondly, the sharp null hypothesis that no subpopulation is affected by the treatment implies that $\mathbb{F}_{Y(0) \mid X} = \mathbb{F}_{Y(1) \mid X}$. With these two assumptions in hand, the TESS framework includes two options for estimating the necessary reference distributions. The first is more flexible but may encounter estimation challenges in extremely sparse, high-dimensional settings; the second is useful for higher-dimensional settings, but requires additional structural assumptions on the data generating process. Finally, given a chosen procedure for estimating reference distributions, TESS uses the distributions to convert each observed outcome $Y^{\text{obs}}_i$ (for data records $R_i$ in the treatment group) to an empirical $p$-value range, capturing how “anomalous” that outcome is given its reference distribution.
The first option we present for deriving reference distributions involves estimating the empirical conditional probability function as follows:
representing a weighted average across data units, with a weight function defined as
For any record $R_j$ in the treatment group (i.e., with $W_j = 1$) we can use $\hat{\mathbb{F}}_{Y^C|X}(\cdot|x_j)$ as an estimate of its distribution function. When we use (ref) as the weight definition, then (ref) amounts to the empirical density function derived from the control units that share covariate profile $X = x_j$. Moreover, it follows directly from TESS's assumption of randomization and the Glivenko-Cantelli Theorem gaenssler-glivenko_cantelli-2004 that $\hat{\mathbb{F}}_{Y^C \mid X} \xrightarrow{a.s.} \mathbb{F}_{Y(0) \mid X}$. Therefore, TESS can use $\hat{\mathbb{F}}_{Y^C \mid X}$ as an unbiased and strongly-consistent estimator of the unknown $\mathbb{F}_{Y(1) \mid X}$ under $H_{0}$. Intuitively, under this sharp null, the outcomes of the treatment and control groups are drawn from the same distribution, allowing $\hat{\mathbb{F}}_{Y^C \mid X = x}$ to serve as an outcome reference distribution for treatment units with covariate profile $X = x$.
Although we define and estimate (ref) individually for each unique covariate profile $X=x$ using the empirical distribution function, we note that TESS only requires some means of computing the conditional probability of observing each treatment unit outcome. The empirical distribution allows TESS to accommodate arbitrary differences in conditional outcome distributions across covariate profiles, enabling general applicability without a priori contextual knowledge. However, it is also possible to combine data across profiles to estimate the conditional probability distributions. This aggregation of information can help alleviate challenges that arise when there is data sparsity, i.e., when there are covariate profiles present in the treatment group that have few or no corresponding control data records. Intuitively, “neighboring” covariate profiles in the control group can be pooled and leveraged to improve local estimation. However, this improved estimation comes by imposing additional structure or assumptions on the underlying data generating process.
Statistical learning offers many options for distribution (or density) estimation, any of which can be utilized in TESS. We identify the Random Forest estimator that underpins the Quantile Regression Forests algorithm meinshausen-quantileforests-2006 as an attractive alternative to the purely empirical estimator described above. Random Forest can be cast as an adaptive locally weighted estimator, where the forest places more weight on observations with more similar covariates. Therefore, TESS can learn a Random Forest on the control data, still using (ref) as its reference distribution, but redefining its weights as:
where $B$ corresponds to the number of trees in the forest, and $L_{b}(x_j)$ captures the leaf node--i.e., a subset of the covariate space--of tree $b$ that $x_j$ falls into. Therefore, (ref) can be seen as a relaxation of (ref), where a control unit can have non-zero weight even if its profile does not match $x_j$ precisely. Also the weights are adaptive and more smoothly increase with how similar a control unit is to the treatment unit. Intuitively, this adaptive similarity “kernel” is particularly helpful in sparse and/or high-dimensional settings, where the curse of dimensionality makes estimation challenging, because it allows for local estimation within covariate subspaces of similar units. Importantly, this similarity is measured along the subset of dimensions that are discovered as relevant (via the random forest learning procedure), which manifests as how often the two points would appear in the same leaf node of the learned trees. Moreover, it has been shown that such a random forest based weighting scheme for estimation can alleviate the curse of dimensionality athey-generalforest-2019. It has also been shown that with weights as in (ref), (ref) is weakly-consistent, given a set of regularity conditions and the assumption that the true distribution function is Lipschitz continuous meinshausen-quantileforests-2006. Therefore, it follows directly from this property of consistency and TESS's assumption of randomization that $\hat{\mathbb{F}}_{Y^C \mid X} \xrightarrow{p} \mathbb{F}_{Y(0) \mid X}$. Therefore, TESS is able to use a random forest estimator of $\hat{\mathbb{F}}_{Y^C \mid X}$ as a weakly-consistent estimator of the unknown $\mathbb{F}_{Y(1) \mid X}$ under the null hypothesis $H_{0}$.
Given a mechanism for estimating the conditional probabilities of outcomes, TESS calculates an empirical $p$-value range mcfowland-fgss-2013 for each treatment unit to obtain a measure of how “anomalous” or unusual a particular unit's outcome is given its reference distribution. For each unit $R_i$ in the treatment group ($W_i=1$), using (ref) and an appropriate weighting scheme, the standard empirical $p$-value would be
The empirical $p$-value range is an extension of this traditional empirical $p$-value, defined as
where the sums are taken over all control observations. The numerator of $\hat{p}_{\text{max}}$, but not $\hat{p}_{\text{min}}$, includes “tied” observations (i.e., $Y^{\text{obs}}_i = y$). The treatment observation $y$ is also considered part of its own reference distribution, following from the assumption of exchangeability of control and treatment outcomes under $H_0$, and thus adding one to the denominators of $\hat{p}_{\text{min}}$ and $\hat{p}_{\text{max}}$ as well as the numerator of $\hat{p}_{\text{max}}$. Following mcfowland-fgss-2013, we use empirical $p$-value ranges because they improve upon traditional empirical $p$-values. The $p$-value ranges are better equipped for sparsity in high-dimensional data, as the range naturally adapts to the amount of reference data used for estimation: a treatment unit's $p$-value range shrinks as more control units are used to estimate its reference distribution. Additionally, if we represented $R_i$ with an empirical $p$-value $\hat{p}_i$ that is drawn uniformly at random from its empirical $p$-value range $\hat{p}(y_i; x_i)$, then under $H_0$, $\hat{p}_i \sim \text{Uniform}(0,1)$.\footnote{From the exchangeability under $H_0$ of $Y^{\text{obs}}_i \sim \mathbb{F}_{Y(1) \mid X}$ and $Y^{\text{obs}}_j \sim \mathbb{F}_{Y(0) \mid X}$, and the probability integral transform.} Standard empirical $p$-values are only asymptotically distributed as Uniform$[0,1]$ and exhibit finite sample bias, while the ranges are unbiased in finite samples, ensuring that $\mathbb{E}\left[ \hat{\mathbb{F}}_{Y (1)| X=x}(\hat{\mathbb{F}}^{-1}_{Y (0)|X=x} (\alpha)) \right] = \alpha$ under $H_0$.
The left-tailed $p$-value ranges defined in (ref) identify outcomes in the extremes of the lower-tail of the reference distribution. The $p$-value range in relation to only the right-tail of the reference distribution or both tails can be derived from the $p$-value range specified for the left-tail. The right-tail range is
while the two-tailed range is
Finally, the significance of a $p$-value range, for a significance level $\alpha$, is defined as
Intuitively, $n_\alpha(\hat{p}(y; x))$ measures the proportion of the range that is significant at level $\alpha$, or equivalently, the probability that a $p$-value drawn uniformly from $\left[\hat{p}_{\text{min}}(y; x), \hat{p}_{\text{max}}(y; x) \right]$ is less than $\alpha$.
Given $p$-values as a measure of the anomalousness of individual treatment units, we now consider how TESS combines these measures to form subpopulations. For intuition, we propose representing the data as a tensor, where each covariate is represented by a mode of the tensor, $X = (X^{1}, \ldots, X^{d})$, resulting in a $d$-order tensor. $|V^j|$, the arity of the $j^{th}$ covariate, is the size of the $j^{th}$ mode. Therefore, each covariate profile $x$ maps to a unique cell in the tensor, which contains the $p$-values of the treatment units that share $x$ as their covariate profile. As stated above, a subpopulation is $S = v^1 \times \ldots \times v^d$, where $v^j \subseteq V^j$; therefore, an individual cell (i.e., covariate profile $x$) is itself a subpopulation: $S = \{x^1\} \times \ldots \times \{x^d\}$, where $x^j \in V^j$. For a demonstrative example see Table (ref). For a given subpopulation $S$, we define the quantities
where $U_{X}(S)$ is the set of non-empty covariate profiles in $S$, $Y^{Tr}(x)= \{Y_i^{obs} | X_i = x, W_i = 1\}$ is the collection of treatment units' outcomes with covariate profile $x$, $N(S)$ represents the total number of empirical $p$-values contained in $S$, and $N_{\alpha}(S)$ is the number of $p$-values in $S$ that are less than $\alpha$.\footnote{For $p$-value ranges, as in mcfowland-fgss-2013, $N_{\alpha}(S)$ is more precisely the total probability mass less than $\alpha$ over the $p$-value ranges in $S$.} Given that the distribution of each $p$-value is Uniform(0,1) under the null hypothesis that the treatment has no effect, for a subpopulation $S$ consisting of $N(S)$ empirical $p$-values, $\mathbb{E}\left[N_{\alpha}(S)\right] = \alpha N(S)$. Under the alternative hypothesis, we expect the outcomes of the affected units to be more concentrated in the tails of their reference distributions; thus, the $p$-values for these affected units will be lower. Therefore, subpopulations composed of covariate profiles that are systematically affected by the treatment should express higher values of $N_\alpha(S)$ for some $\alpha$. Consequently, a subpopulation $S$ where $N_\alpha(S) > \alpha N(S)$ (i.e., with a higher than expected number of low, significant $p$-values) is potentially affected by the treatment.
TESS utilizes the nonparametric scan statistic mcfowland-fgss-2013,feng-npss_graph-2014 to evaluate the statistical anomalousness of a subpopulation $S$ by comparing the observed and expected number of significantly low $p$-values it contains. The general form of the nonparametric scan statistic is
where $N_{\alpha}(S)$ and $N(S)$ are defined as in (ref), and $\Delta$ is a measure of divergence measuring the anomalousness of the $p$-values in $S$. See Appendix (ref)for a collection of goodness-of-fit scoring functions written in the general form of the nonparametric scan statistic. As described in (ref), in this work we utilize the Berk-Jones scan statistic: $\max_{\alpha} \Delta_{BJ}\left(\alpha,N_{\alpha}(S),N(S)\right) = \max_{\alpha}N(S)KL\left( \frac{N_{\alpha}(S)}{N(S)} , \alpha \right)$, a log-likelihood ratio test statistic of the distributional treatment effect in subpopulation $S$. Maximizing $F(S)$ over a range of $\alpha$, rather than a single arbitrarily-chosen $\alpha$ value, enables TESS to detect a small number of highly anomalous $p$-values, a larger subpopulation with subtly anomalous $p$-values, or anything in between. We consider “significance levels” $\alpha \in [\alpha_{\text{min}},\alpha_{\text{max}}]$, for constants $0 < \alpha_{\text{min}} < \alpha_{\text{max}} < 1$. The range of $\alpha$ to consider can be specified based on the quantile values of interest. The choice of $\alpha_{\text{max}}$ describes how extreme a value must be, as compared to the reference distribution, in order to be considered significant. We often choose $\alpha_{\text{min}} \approx 0$, but larger values can be used to avoid returning subsets with a small number of extremely significant $p$-values.
The next step in the TESS framework is to detect the subpopulation most affected by the treatment, i.e., to identify the most anomalous subset of values for each of the $d$ modes of the tensor, or equivalently for each covariate $X^1 \ldots X^d$. More specifically, the goal is to identify the set of subsets $\{v^{1},\ldots, v^{d}\}$ where each element corresponds to values in a tensor-mode (covariate), such that $F(v^{1} \times \ldots \times v^{d})$ is jointly maximized. The computational complexity of solving this optimization naively is $O^{(2^{\sum_j |V^j|})}$, where $|V^j|$ is the size of mode $j$ (the arity of $X^j$), and is computationally infeasible for even moderately sized datasets.
We therefore employ the linear-time subset scanning property (LTSS) neill-ltss-2012, which allows for efficient and exact maximization of any function satisfying LTSS over all subsets of the data. We formally define the LTSS property below, but intuitively it guarantees that the optimization over all subsets $S$ can be done by ranking data elements (according to a specific “priority function”) and then only considering the top-$t$ subsets as candidates.\\[1ex] LTSS Property Definition: Given a set of data elements $R=\{R_1, \ldots, R_n \}$, a score function $F(S)$ mapping $S\subseteq R$ to a real number, and a priority function $G(R_i)$ mapping a single data element $R_i \in R$ to a real number. If $F(S)$ satisfies the LTSS property with priority function $G(R_i)$, then the only subsets with the potential to be optimal are those consisting of the top-$t$ highest priority records, $S \in \{\{R_{(1)}, \ldots, R_{(t)}\} \}_{t\in\{1,2,\ldots,n\}}$. In other words, there exists some $t \in \{1,2,\ldots,n\}$ such that $\arg\max_S F(S) = \{R_{(1)}, \ldots, R_{(t)}\}$. We also formally restate the original LTSS theorem:
In this work, we use Theorem (ref) to optimize $F_\alpha(S) = \Delta(\alpha,N_\alpha(S),N(S))$ with a fixed value of $\alpha$; therefore, $X(S) = N_\alpha(S)$ and $Y(S) = N(S)$ are “additive sufficient statistics”, i.e., both $N_\alpha(S)$ and $N(S)$ are additive statistics of $S$, from (ref), and $F_\alpha(S)$ can be written as $F_\alpha(N_{\alpha}(S), N(S))$. Moreover, for $F_\alpha(S)$ (with $\alpha$ fixed) to satisfy LTSS, we also assume: (A1) $\Delta$ is monotonically increasing w.r.t. $N_{\alpha}$, (A2) $\Delta$ is monotonically decreasing w.r.t. $N$, and \emph{(A3)} $\Delta$ is \emph{\textbf{convex}} w.r.t. $N_{\alpha}$ and $N$. These properties are intuitive because the ratio of observed to expected number of significant $p$-values $\frac{N_\alpha}{\alpha N}$ increases with the numerator (A1) and decreases with the denominator (A2). Also, a fixed ratio of observed to expected is more significant when the observed and expected counts are large (A3). In Appendix (ref), we show that nonparametric scan statistics using a large class of goodness of fit functions, including the Berk-Jones scan statistic utilized in this work, exhibit these properties.
We now extend Theorem (ref) to the (potentially high-dimensional) tensor context using Corollary (ref) below. Essentially, the corollary demonstrates that the nonparametric scan statistic satisfies LTSS in the context of TESS, and therefore a single mode of a tensor can be efficiently optimized over subsets, conditioned on the subsets of values for the other modes. Let $U_{\alpha}(S)$ be the set of unique $p$-values between $\alpha_{\text{min}}$ and $\alpha_{\text{max}}$ contained in subpopulation $S$. Then the quantity $\max_S F(S) = \max_{\alpha \in U_{\alpha}(S)} \max_S F_\alpha(S)$ can be efficiently and exactly computed over all subsets $S = v^j \times v^{-j}$, where $v^j \subseteq V^j$, for a given subset of values for each of the other modes $v^{-j}$.\footnote{Note that for convenience of notation we define $S = v^j \times v^{-j}$; however, the elements of the set $v^j$ still appear at the $j^{\text{th}}$ position of the covariate profiles in $S$.} To do so, consider the set of distinct $\alpha$ values, $U = U_{\alpha}(V^j \times v^{-j})$. For each $\alpha \in U$ we employ the logic described in Corollary (ref) to optimize $F_{\alpha}(S)$: we compute the priority $G_\alpha(v^{j}_m)$ for each value ($v^{j}_m \in V^j$), sort the values based on priority function $G_{\alpha}(v^{j}_m)$, and evaluate subsets of the form $S=\{v^{j}_{(1)}, \ldots, v^{j}_{(t)}\} \times v^{-j}$ consisting of the top-$t$ highest priority values, for $t=1, \ldots,|V^j|$.
TESS iterates over modes of the tensor, using the efficient optimization steps described above to optimize each mode: $v^j = {\arg\max}_{v^j \subseteq V^j} F(v^j \times v^{-j}), j = 1 \ldots d$. The cycle of optimizing each mode continues until convergence, at which point TESS has reached a conditional maximum of the score function, i.e., $v^j$ is conditionally optimal given $v^{-j}$ for all $j = 1 \ldots d$. This ordinal ascent approach is not guaranteed to converge to the joint optimum, but with multiple random restarts the combination of subset scanning and ordinal ascent has been shown to locate near globally optimal subsets with high probability neill-mvltss-2013, mcfowland-fgss-2013. We further provide asymptotic guarantees for TESS to recover the precisely correct subpopulation that is also shown to be the globally optimal subset (Theorems (ref) and (ref)). Moreover, if $\sum_{j=1}^d |V^j|$ is large, this iterative procedure makes the ability to detect anomalous subpopulations computationally feasible, without excluding potentially optimal subpopulations from the search space (as a greedy top-down approach may). A single iteration (optimization of mode $j$ of the tensor) has a complexity of $O\left(|U| \left(n_t + |V^j| \log |V^j|\right)\right)$, where the $n_t$ term---the number of treatment units---results from collecting the $p$-values for all units in $V^j \times v^{-j}$ over our sparse tensor; $U = U_{\alpha}\left(V^j \times v^{-j}\right)$, with $|U| \le n_t$ mcfowland-fgss-2013; and $O\left(|V^j| \log |V^j|\right)$ is required to sort, based on the priority, the values of tensor mode $j$. Therefore a step in the procedure (a sequence of $d$ iterations over all modes of the tensor) has complexity $O\left(\bar{U} d \left(n_t + \bar{V} \log \bar{V}\right)\right)$, where $\bar{U}$ and $\bar{V}$ are the average numbers of $\alpha$ thresholds considered and covariate arity, respectively. Thus the TESS search procedure has a total complexity of $O\left(I \bar{Z} \bar{U} d \left(n_t + \bar{V} \log \bar{V}\right)\right)$, where $I$ is the number of random restarts and $\bar{Z}$ is the average number of iterations required for convergence. We note that $\bar{Z}$ is typically very small; $\bar{Z} \le 5$ across all simulations discussed in \S(ref).
Inputs: randomized experiment dataset, $\alpha_{\text{min}}$, $\alpha_{\text{max}}$, number of iterations $I$.
In the above sections we outline a procedure to efficiently estimate $\max_{S \in Rect} F(S)$, where $Rect$ represents the space of all rectangular subsets of $D$. In this section we treat $\max_{S \in Rect} F(S)$ as a statistic of the data, and aim to show that it has desirable statistical properties. It is known that for data $X_1,\ldots,X_n \overset{iid}{\sim} \mathbb{F}$ and the corresponding empirical distribution function $\mathbb{F}_n$, $\|\mathbb{F}_n - \mathbb{F}\|_\infty \xrightarrow{a.s.}\ 0$. Many goodness-of-fit statistics $GoF(\mathbb{F}_n, \mathbb{F})$ are equivalent to an empirical process over centered and scaled empirical measures; and empirical process theory provides tools to control Type I and II error gaenssler-glivenko_cantelli-2004,dvoretzky-dkw-1956,shorack-emp_process-1986. However, in a general sense our goal is to control the behavior of $\max_{S \subseteq \{X_1,\ldots,X_n\}}GoF(\mathbb{F}_S, \mathbb{F})$, where $\mathbb{F}_S$ is the empirical distribution given by the subset $S$. It is not obvious whether the desirable properties present for $\mathbb{F}_n$ will persist when considering the empirical distribution of $\mathbb{F}_S$, a non-random subset of the data chosen by our optimization procedure. Given that this context of optimization over subsets is not considered in the current goodness-of-fit literature, we provide various theoretical results in support of our subset scanning algorithm. In the remainder of the section we present the key statements necessary to show our desired properties below, while additional results and all proofs can be found in Appendix (ref). We begin with
the fact that our score function can be considered a test statistic for a hypothesis test analogous to that described in (ref):
where $D$ is our dataset (or tensor) of treatment units and $Rect$ is the set of all rectangular subsets of $D$.\footnote{The null hypotheses defined in (ref) and (ref) are analogous because if $\tau_{\text{CQTE}_{\alpha}}(x) = \alpha~ ~ \forall x\in S, \forall \alpha$, then the distributions for $Y(1) | X$ and $Y(0) |X$ are the same. Similarly, the alternative hypotheses are analogous because if $\exists \beta > \alpha$ such that $\tau_{\text{CQTE}_{\alpha}}(x) = \beta~ \forall x \in S$, then $Y(1) | X$ is different from $Y(0) | X$, and if $\tau_{\text{CQTE}_{\alpha}}(x) = \alpha~~ \forall x \notin S,~ \forall \alpha$, then again the distributions for $Y(1) | X$ and $Y(0) |X$ are the same.} The null hypothesis is that all of the observed outcomes of treatment units are drawn from the same conditional outcome distribution (given the observed covariates) as their control group.
Recall that $U_{X}(D)$ is the set of unique covariate profiles (non-empty tensor cells) in our data, with cardinality $|U_{X}(D)| = M$; while $S^{\ast}=\arg \max_{S \in Rect} F(S)$ and $S^{\ast}_{u} = \arg \max_S F(S)$ represent the most anomalous rectangularly constrained subset and the most anomalous unconstrained subset respectively. For mathematical convenience, we assume $N(x) = n$ for all $x \in U_{X}(D)$, i.e., $n$ units belong to each unique covariate profile (non-empty cell) in the data and treatment condition.\footnote{We can redefine $n = \min_{x} N(x)~\forall x \in U_{X}(D)$ and our results can be extended.} We consider the case where $n \longrightarrow \infty$, maximizing $F(S)= \max_{\alpha\in[\alpha_{\min},\alpha_{\max}]} F_\alpha(S)$ for $0 < \alpha_{\min} < \alpha_{\max} < 1$. We can therefore demonstrate:
Thus, when the null hypothesis is true, the most anomalous unconstrained subset's score distribution can be upper bounded. Our ability to understand the limiting behavior of the $F\left(S^{\ast}_{u}\right)$ exploits its structure, which we get from LTSS theory: the optimal unconstrained subset will be $S^{\ast}_{u} \in \{\{x_{(1)}, \ldots, x_{(t)}\}\}_{t\in\left\{1,2,\ldots,M\right\}}$, where $x_{(t)}$ has the $t^{th}$ largest value of the random variable $\frac{N_{\alpha}(x)}{N(x)}~\forall x \in U_{X}(D)$. Next we note that the score maximized over the space of unconstrained subsets upper bounds the score maximized over the subspace of rectangular subsets, i.e., $F\left(S^{\ast}\right) \le F\left(S^{\ast}_{u}\right)$. We use this fact to obtain the following result:
Theorem (ref) indicates that $\max_{S \in Rect} F(S)$ provides a statistic to quantify the evidence to reject $H_{0}$, enabling an (asymptotically) valid $\delta$-level hypothesis test such that $P_{H_0}(\text{Reject }H_0) \le \delta$, for any fixed Type I error rate $\delta > 0$. From the proof of Theorem (ref), in Appendix (ref), we derive that $h(\delta)= \left(0.45\sqrt{M} + \frac{w(\delta)}{\sqrt{2}}\right)^2$, where $w(\delta)$ returns $w$ such that $P\left(\mathbb{W}(\alpha_{\min},\alpha_{\max}) > w\right) = \delta$. For intuition, $w(\delta)$ is typically small, e.g., $w(\delta) \approx 5.81$ for $\delta = 10^{-6}$ and $(\alpha_{\min}, \alpha_{\max}) = (.01, .99)$. We note that because we are maximizing both over subsets $S$ and thresholds $\alpha$, these results are distinct from the straightforward application of known results from empirical process theory or Dvoretzky-Kiefer-Wolfowitz bounds, which would give us $\max_\alpha | \frac{N_\alpha(S)}{N(S)} - \alpha| \longrightarrow 0$ for a given $S$.
Next, we turn our attention to the alternative hypothesis, where $S^T \in Rect$ represents the truly affected subset, $k=\frac{|U_{X}(S^T)|}{|U_{X}(D)|}$ is the proportion of covariate profiles included in $S^T$, and $H_{1}\left(S^T\right)$ implies that there exist constants $\alpha$ and $\beta(\alpha) > \alpha$ such that $\beta(\alpha) = \mathbb{F}_{Y(1) \mid X=x}\left(\mathbb{F}_{Y(0) \mid X=x}^{-1}(\alpha)\right)$ for all $x \in U_{X}(S^T)$. We then have the following results:
Thus, when the null hypothesis is false, the expected value of the true subset's score at the $\alpha^{\ast}$ quantile, $F_{\alpha^{\ast}}(S^T)$, is increasing with $n$. This result, and the fact that $F_{\alpha}(S) \le F(S)~\forall S,\alpha$, are used to obtain the following result:
As a consequence of Theorem (ref), the $\delta$-level hypothesis test based on $\max_{S \in Rect} F(S)$ has full asymptotic power $P_{H_1}(\text{Reject }H_0) \longrightarrow 1$. Note that in this context we consider a fixed alternative $\beta(\alpha)$, as opposed to a local alternative where $\beta_{n}(\alpha) \longrightarrow \alpha$ as $n \longrightarrow \infty$.
In practice, when we consider an experiment with finite $M$ and $n$, permutation testing can be used to control the Type I error rate of our scanning procedure, and conditions have been shown where permutation calibrations achieve the Type II error rates of an oracle scan test arias-castro-calibration-2017. These theoretical and practical results intuitively capture our statistic's ability to conclude that the null hypothesis is false--i.e., that there exists some subset that follows $H_1$, and therefore invalidates $H_{0}$. However, this does not necessarily provide a guarantee that the statistic will exactly capture the true subset. Therefore, next we derive finite sample sufficient conditions under which our framework achieves subset correctness: $S^{\ast} = S^T$. We then show asymptotic convergence of $P(S^{\ast} = S^T ) \xrightarrow[]{} 1$ as $n\xrightarrow[]{}\infty$.
We introduce additional notation for this discussion: $r_{\text{mle}}(x) = \frac{N_{\alpha}(x)}{N(x)} - \alpha$, $r^{\text{aff}}_{\text{mle}-h} = \max_{x \in U_{X}(S^T)} r_{\text{mle}}(x)$, $r^{\text{aff}}_{\text{mle}-l} = \min_{x\in U_{X}(S^T)} r_{\text{mle}}(x)$, $r^{\text{unaff}}_{\text{mle}-h} = \max_{x \not\in U_{X}(S^T)} r_{\text{mle}}(x)$, and $\eta = \left( \frac{\sum_{x \in U_{X}(S^T)}{N\left(x\right)} }{ \sum_{x \in U_{X}(D)}{N\left(x\right)} } \right)$. We also introduce the concepts of $\nu-homogeneous$, which means that $\frac{r^{\text{aff}}_{\text{mle}-h}}{r^{\text{aff}}_{\text{mle}-l}} < \nu$, and $\delta-strong$, which means that $ \frac{r^{\text{aff}}_{\text{mle}-l}}{r^{\text{unaff}}_{\text{mle}-h}} > \delta$. Intuitively, the concept of homogeneity measures how similarly the treatment affects each $\mathbb{F}_{Y| X=x}$ for $x \in U_{X}(S^T)$, while strength measures how large of an effect the treatment exhibits across all $\mathbb{F}_{Y | X=x}$ for $x \in U_{X}(S^T)$. More specifically, these concepts respectively imply that for any pair of the affected covariate profiles $\left(x_i, x_j \in U_{X}(S^T)\right)$, the anomalous signal (i.e., treatment effect) observed in $x_i$ is less than $\nu$ times that which is observed in $x_j$, and the treatment effect observed in every affected covariate profile is more than $\delta$ times that of the unaffected profiles. Using these concepts we have the following results:
While Theorem (ref) shows that $S^\ast = S^T$ with high probability as $n\rightarrow\infty$---i..e, the most anomalous (rectangular) subset is the true subset---it does not guarantee that the TESS algorithm presented in Section (ref) will converge to the true subset, because iterative ascent algorithms converge to a local maximum. For example, if Step 2(a) of the algorithm chooses an initial subset $S_0$ that is disjoint from $S^T$, then it is possible that no maximization over subsets of values for any single mode in Step 2(b).i.A. will improve the score, and TESS will fail to identify $S^T$ on that iteration. However, we can show the following:
Thus as $n\rightarrow\infty$ for fixed $M$, TESS will identify the correct rectangular subset $S^T$ w.h.p., as long as it is initialized to some rectangular subset $S_0$ that overlaps $S^T$, for at least one of the its iterations. For example, $S_0$ could be the entire dataset $D$, thus guaranteeing $S_0 \cap S^T \ne \emptyset$.
Together, these results demonstrate that the test statistic $F^\ast = \max_{S \in Rect} F(S)$ and corresponding subset $S^\ast = \arg\max_{S \in Rect} F(S)$ possess desirable statistical properties. Theorems (ref) and (ref) imply that the asymptotic Type I and II errors of our procedure can be controlled, with implications for maximization over subsets of empirical processes more generally. Theorems (ref) and (ref) indicate that for a score function there exist constants $\nu$ and $\delta$ that define how similar and strong the treatment effect must be in the affected subpopulation, to ensure that the highest-scoring subset corresponds exactly to the true affected subset $S^{\ast} = S^T$. Finally, Theorem (ref) shows asymptotic convergence for $P(S^{\ast} = S^T) \xrightarrow[]{} 1$ as $n\xrightarrow[]{}\infty$, and Theorem (ref) shows that the TESS algorithm will identify $S^T$ w.h.p. as $n\rightarrow\infty$. To our knowledge, this is the first work on heterogeneous treatment effects that provides conditions on the exactness of subpopulation discovery.
There has been a growing literature using statistical learning methods to provide data-driven approaches for estimating heterogeneous treatment effects in randomized experiments. Recent work has adapted regularized regression for treatment effect heterogeneity imai-hte_lasso-2013,tian-hte_lasso-2014,weisberg-subgroup-2015. These regularized regression approaches, however, require the researcher to select which covariate and treatment interactions to include in the model specification, compromising their ability to discover unexpected treatment patterns in subpopulations.
Regression tree based methods su-subgroup-2009,athey-hte-2016 select subpopulations and estimate treatment effects by recursively partitioning the data into homogeneous subpopulations that share a subset of covariate profile values and have similar outcomes. The effectiveness of tree methods can be severely compromised in many settings as a result of their greedy partitioning.\footnote{Though there are non-greedy tree based learning methods, these methods are used for optimal treatment assignments for individual units in observational data zhou-offline-policy-2018.} Tree models can be unstable; they can provide extremely discontinuous approximations of an underlying smooth function, limiting overall accuracy; and they can struggle to estimate functions which exhibit specific properties, including when a small proportion of the covariates constitute the influential interactions friedman-mars-1991.
Other treatment effect estimation approaches use ensemble methods, including the use of Bayesian Additive Regression Trees hill-hte_bart-2011, green-hte_bart-2012, Random Forests foster-hte_randomforest-2011, wager-causalforest-2018, and ensembles of strong learners grimmer-hte_ensemble-2017. Ensembles provide more stable and smooth function estimates wager-causalforest-2018; however, they lose the interpretability of natural groupings (e.g., specific combinations of covariates or clearly defined leaves) which is important for identifying affected subpopulations.
Finally, chernozhukov-generic-2018 propose approaches to test if there is detectable heterogeneity in treatment effects, finding the quantiles exhibiting heterogeneous treatment effects induced by the machine learning proxy predictors, and then identifying the covariates that appear associated with the heterogeneity. Therefore this approach is a post-hoc analysis of existing machine learning predictors (e.g., regression tree or ensemble methods) which are optimized for overall risk minimization and not the discovery of subpopulations with significant treatment effects. Furthermore, there are orthogonal research streams focused on (heterogeneous) treatment assignment that include policy learning zhou-offline-policy-2018, athey-policy-learning-2020 and welfare maximization toru-welfare-max-2018, mbakop-welfare-max-2020 in observational and experimental data. These methods assign a personalized treatment for each individual unit. Instead, TESS takes as input a set of units which have already been randomly assigned a treatment.
Although this literature contains a growing set of novel statistical learning methods for causal inference, at the core of the majority of these approaches are objective functions designed for flexible estimation (and risk minimization) instead of subpopulation discovery (and significant effect maximization). TESS is therefore unique as it is optimized to discover interpretable subpopulations that exhibit significant evidence of treatment effects. When necessary, TESS can use flexible (risk minimizing) statistical learning models for a purpose that is aligned with their objective: as an accurate estimator of the conditional outcome distribution for the control group as in Section (ref).
In this section we empirically demonstrate the utility of the TESS framework as a tool to identify subpopulations with significant treatment effects. We use data from the Tennessee Student/Teacher Achievement Ratio (STAR) randomized experiment word-star-1990 in order to provide representative performance in real-world policy analysis. We review the original STAR data (\S(ref)), and describe our procedure for simulating affected subpopulations (\S(ref)).
Through the simulation results described in \S(ref), we compare the ability of TESS to detect significant subpopulations to three recently proposed statistical learning approaches: Causal Tree athey-hte-2016, Interaction Tree su-subgroup-2009,athey-hte-2016, and Causal Forest wager-causalforest-2018. Specifically, we evaluate each method on two general metrics: detection power and subpopulation accuracy. Detection power measures $P_{H_1}(Reject~H_0)$, or how well a method can detect the existence--not necessarily the location--of treatment effect heterogeneity in the experiment. Subpopulation accuracy, on the other hand, is specifically designed to measure how well a method can precisely and completely capture the subpopulation(s) with significant treatment effects.
Finally, we conduct an exploratory analysis of the STAR dataset, and in \S(ref) discuss the subpopulations identified by TESS as affected by treatments. In some cases, the identified subpopulation is consistent with the literature on the STAR experiment; in other cases, TESS uncovers previously unreported, but intuitive and believable, subpopulations. These empirical results demonstrate TESS's potential to generate potentially useful and non-obvious hypotheses for further exploration and testing.
The Tennessee Student/Teacher Achievement Ratio (STAR) experiment is a large-scale, four-year, longitudinal randomized experiment started in 1985 by the Tennessee legislature to measure the effect of class size on student educational outcomes, as measured by standardized test scores. The experiment started monitoring students in kindergarten (during the 1985-1986 school year) and followed students until third grade. Students and teachers were randomly assigned into conditions during the first school year, with the intention for students to continue in their class-size condition for the entirety of the experiment. The three potential experiment conditions were not based solely on class size, but also the presence of a full-time teaching aide: small classrooms (13-17 pupils), regular-size classrooms (22-25 pupils), and regular-size classrooms with aide (still 22-25 pupils). Therefore, the difference between the former two conditions is classroom size, and the difference between the latter two conditions is the inclusion of a full-time teacher's aide in the classroom. The experiment included approximately 350 classrooms from 80 schools, each of which had at least one classroom of each type. Each year more than 6,000 students participated in this experiment, with the final sample including approximately 11,600 unique students.
The Tennessee STAR dataset has been well studied and analyzed, both by the project's internal research team word-star-1990, folger-star-1989 and by external researchers krueger-star-1999,nye-star-2000. As indicated by krueger-star-1999, the investigations have primarily focused on comparing means and computing average treatment effects. krueger-star-1999 presents a detailed econometric analysis and draws similar conclusions to the previous research: students in small classrooms perform better than those in regular classrooms, while there is no significant effect of a full-time teacher's aide, or moderation from teacher characteristics. Moreover, the effect accumulates each year a student spends in a small classroom krueger-star-1999. Additionally, these conclusions are robust in the presence of potentially compromising experimental design challenges: imbalanced attrition, subsequent changes in original treatment assignment, and fluctuating class sizes krueger-star-1999.
The goal of our experimental simulation is to replicate conditions under which a researcher would want to use an algorithm to discover subpopulations with significant treatment effects, and to observe how capable various algorithms are at identifying the correct subpopulation(s). In order to replicate realistic conditions, we use the STAR experiment as our base dataset, and inject into it subpopulations (of a given size) with a treatment effect (of a given magnitude). More specifically, we treat each student-year as a unique record and for each record capture ten covariates: student gender, student ethnicity, grade, STAR treatment condition, free-lunch indicator, school development environment, teacher degree, teacher ladder, teacher experience, and teacher ethnicity. We note that each of these variables, other than teacher experience, is discrete; we discretize experience into five-year intervals: $[0,5), [5,10), \ldots, [30, \infty)$. The number of values a covariate can take ranges from two to eight. By preserving the overall data structure of the STAR experiment--number of covariates, covariate value correlations, subpopulations, sample sizes, etc.--our simulations are more able to replicate the structure (and challenges) faced by experimenters.
The process we follow to generate a simulated treatment effect begins with selecting a subpopulation $S_{\text{affected}}$ to affect. Recall that the dataset contains a set of discrete covariates $X = (X^{1}, \ldots, X^{d}$), where each $X^{j}$ can take on a vector of values $V^{j}=\{v^{j}_{m}\}_{m = 1...|V^j|}$ and $|V^{j}|$ is the arity of covariate $X^{j}$. Therefore, we define a subpopulation as $S = v^1 \times \ldots \times v^d$, where $v^j \subseteq V^j$. The affected subpopulation is generated at random based on two parameters: $num\_covs$, or the number of covariates to select, and $value\_prob$, or probability a covariate value is selected. We select $num\_covs$ covariates at random, and for each of these covariates we select each of their values with probability equal to $value\_prob$, ensuring that at least one value for each of these covariates is selected. The final affected subpopulation is then $S_{\text{affected}} = v^1 \times \ldots \times v^d$, where $v^j$ is the selected values if $X^{j}$ is one of the $num\_covs$ covariates, and otherwise $v_j = V^j$. In other words, for a random subset of covariates, $S_{\text{affected}}$ only includes a random subset of their values, and for all other covariates $S_{\text{affected}}$ includes all of their values. This treatment effect simulation scheme allows for variation in the size of the subpopulation that is affected: instances of $S_{\text{affected}}$ can constitute a small subpopulation (a challenging detection task), a large subpopulation (a relatively easier detection task), or something in between. Therefore a set of simulations, with varying parameter values, captures the spectrum of conditions a researcher may face when analyzing an experiment to identify subpopulations with significant treatment effects.
The next step in the process involves partitioning the dataset into treatment and control groups, and generating outcomes for each record. Outcomes are drawn randomly from one of two distributions: the null distribution ($f_0$) or the alternative distribution ($f_1$). Any record in the treatment group that has a covariate profile $x \in U_{X}\left(S_{\text{affected}}\right)$ has outcomes generated by $f_1$; all other records have outcomes drawn from $f_0$. Therefore only $S_{\text{affected}}$ has a treatment effect, whose effect magnitude is the distributional difference between $f_0$ and $f_1$, represented by the parameter $\delta$.
Each of the methods we consider in these experiments has a unique approach to identifying potential subpopulations with differential treatment effects. Furthermore, as mentioned in \S(ref), most methods in the literature do not provide a process for identifying extreme treatment effects. Therefore, we devise intuitive post-processing steps in an attempt to represent how researchers would use each method to identify potential subpopulations that have significant treatment effects. Each method returns identified subpopulations and corresponding scores (measures) of the treatment effect. For the single tree-based methods athey-hte-2016,su-subgroup-2009 we follow the suggestion of athey-hte-2016 to perform inference (via a two-sample Welch T-Test) in each leaf of the tree, and we then sort the leaves based on their statistical significance. The final subpopulation returned by the tree is the leaf with the most statistically significant treatment effect, and the final treatment effect measure is this leaf's statistical significance ($p$-value). For a method that provides an individual level treatment effect (and estimate of variance) wager-causalforest-2018, we propose to perform inference for each unique covariate profile, and return those that are statistically significant. The final treatment effect measure is the smallest $p$-value of the covariate profiles. We empirically compare these prominent methods from the literature to our TESS algorithm, selecting the Berk-Jones nonparametric scan statistic to score a subset (\S(ref)) and the most general empirical estimation approach to model the reference distribution (\S(ref)). The TESS algorithm, by design, provides the subpopulation it determines to have a statistically significant distributional change (treatment effect) and a measure of this change, so no post-processing is necessary.
For any given combination of simulation parameter values ($\delta, num\_covs, value\_prob$), detection power measures $P(Reject~H_0 \mid H_1(S_{\text{affected}}))$, or how well a method is able to identify the presence of $S_{\text{affected}}$. This is accomplished by comparing the treatment effect measure (score of the detected subset) found under $H_1(S_{\text{affected}})$ to the distribution of the treatment effect measure under $H_0$. More specifically, for a given set of parameter values, we generate a random dataset which only exhibits a treatment effect in the randomly selected subpopulation $S_{\text{affected}}$; each method attempts to detect this subpopulation. As described in \S(ref), each method returns a final treatment effect measure for the subpopulation it detects in this affected dataset. For the same dataset, we then conduct randomization testing to determine how significant this treatment effect measure is under $H_0$. We make many copies of the dataset (1000 in our experiments) and in each copy, we generate new outcomes (drawn from $f_0$) such that no subpopulation has a treatment effect. Each method then generates a detected subpopulation and corresponding treatment effect measure for each of these null datasets. These treatment effect measures from the null datasets together provide an empirical estimate of the distribution of the treatment effect measure under $H_0$ for that method. Subsequently, a $p$-value is computed for the treatment effect measure captured under $H_1(S_{\text{affected}})$. This process is repeated many times (300 in our experiments), where each time we 1) generate a random $S_{\text{affected}}$, 2) generate a random dataset under $H_1(S_{\text{affected}})$ and compute each method's treatment effect measure, and 3) generate 1000 copies of the dataset with no treatment effect to compute each method's treatment effect measure distribution under $H_0$. This process creates 300 $p$-values for each method which describe how extreme each of the $S_{\text{affected}}$ appear under $H_0$. A method rejects $H_0$ for a given $p$-value if it is less than or equal to some test-level $\gamma$, corresponding to the $1-\gamma$ quantile of the null distribution ($\gamma = 0.05$ in our experiments). Therefore, the detection power $P\left(Reject~H_0 \mid H_1(S_{\text{affected}})\right)$ is captured as the proportion of $p$-values that are sufficiently extreme that they lead to the rejection of $H_0$ at level $\gamma$.
While detection power measures how well a method identifies the presence of a subpopulation with a treatment effect $S_{\text{affected}}$, as compared to datasets with no treatment effect, detection accuracy measures how well a method can precisely and completely identify the affected subpopulation $S_{\text{affected}}$. Accurately identifying in which subpopulation(s) a treatment effect exists can be crucial, particularly when there is no prior theory to guide which subpopulations to inspect, or when the goal itself is to develop intuition for new theory. As described in \S(ref), each of the methods we consider is able to return the subpopulation that it determines as having the most statistically significant treatment effect $S_{\text{detected}}$. Each method will pick out a set of covariate profiles, which could have coherent structure (as with TESS, Causal Tree, and Interaction Tree), or be an unstructured collection of individually significant covariate profiles (as with Causal Forest). To accommodate both types of subpopulations, we therefore define detection accuracy as
where $R_i$ are records in the treatment group. This definition of accuracy, commonly known as the Jaccard coefficient, is intended to balance precision (i.e., what proportion of the detected subjects truly have a treatment effect) and recall (i.e., what proportion of the subjects with a treatment effect are correctly detected). We note that $0 \le \text{accuracy} \le 1$; high accuracy values correspond to a detected subset $S_{\text{detected}}$ that captures many of the subjects with treatment effects and few or no subjects without treatment effects.
Our first set of results involve a treatment effect that is a mean shift in a normal distribution: the null distribution $f_0 = N(0,1)$ and the alternative $f_1= N(\delta,1)$, where $\delta$ captures the magnitude of the signal (treatment effect). Recall from \S(ref) that there are three parameters that we can vary to change the size and magnitude of the signal. For our simulation, we specifically consider $\delta \in \{0.25, 0.5, \ldots, 3.0\}$, $num\_covs \in \{1,2, \ldots, 10\}$, and $value\_prob \in \{0.1, 0.2, \ldots, 0.9\}$; the former controls magnitude of the treatment effect, while the latter two control the concentration of the treatment effect (i.e., the expected size of the affected subpopulation). Instead of considering every combination, we select the middle value of each parameter interval as a reference point ($\delta= 1.5, num\_covs = 5, value\_prob = 0.5$) and measure performance changes for one parameter, while keeping the others fixed.
Figure (ref) shows the changes in each method's detection power performance as we vary each of the three parameters that contribute to the strength of the treatment effect. From each of the three graphs we observe that TESS consistently exhibits more power than (or equivalent to) the other methods. More importantly, TESS exhibits statistically significant improvements in power for the most challenging ranges of parameter values (i.e., more subtle signals). The top plot varies effect size (or $\delta$), which is positively associated with signal strength and negatively associated with detection difficulty; for values $2.0$ and below TESS has significantly higher detection power than the competing methods. The middle plot varies the number of covariates selected to have only a subset of values be affected ($num\_covs$). This parameter is negatively associated with signal strength and positively associated with detection difficulty; for values $5$ and above, TESS has significantly higher detection power. The bottom plot varies the expected proportion of values, for the selected covariates, which will be affected ($value\_prob$). This parameter is positively associated with signal strength and negatively associated with detection difficulty; for values $0.5$ and below TESS exhibits significantly higher detection power. We see that, for sufficiently strong signals (based on both signal magnitude and concentration), all methods are able to distinguish between experiments with and without a subpopulation exhibiting a treatment effect, while TESS provides significant advantages in detection power for weaker signals.
Figure (ref) shows the changes in each method's detection accuracy as we vary each of the three parameters that contribute to the strength of the treatment effect. From each of the three graphs we observe that TESS consistently exhibits significantly higher accuracy than any other method. Recall that we measure subpopulation accuracy as in (ref), which captures both precision and recall of the subpopulation returned by a method. The single tree methods tend to have high precision but low recall, resulting in compromised overall accuracy. Intuitively, these results indicate that the truly affected subpopulation is being spread over multiple leaves of the tree, despite its goal of partitioning the data into subpopulations with similar outcomes. This phenomenon may be caused by the greedy search aspect of tree learning: if the tree splits the affected subpopulation between two branches of the tree, the recall of any leaf will be compromised, especially when this split occurs close to the root of the tree. The Causal Forest ensemble method, on the other hand, exhibits relatively higher recall than precision. These results indicate that it is difficult for Causal Forest to distinguish between the covariate profiles that do and do not make up the truly affected subpopulation, as profiles from both sets appear to have statistically significant treatment effects. This inability stems from the fact that ensemble methods are designed to provide individual level predictions, therefore their conclusions regarding the statistical significance of a covariate profile are made in isolation from the other covariate profiles that also make up the affected subpopulation. Unlike single-tree methods, ensemble methods do not provide coherent and natural groupings of subpopulations. TESS, however, does provide a coherent subpopulation, which seems to balance precision and recall, maintaining a significantly higher subpopulation accuracy.
It is also important to note that the data generating process for these simulations (a treatment effect that occurs as a mean shift between treatment and control distributions) corresponds to the modeling assumptions of the current methods in the literature, which specifically attempt to detect mean shifts, while TESS is designed to detect more general distributional changes. TESS's improved performance, as compared to the competing methods, in these adverse conditions may be due to its subset-scanning based approach, which combines information across groups of data in an attempt to find exactly and only the affected subset of data. Even if each individual covariate profile that is truly affected exhibits small evidence of a treatment effect, TESS can leverage the group structure and signal of all the affected covariate profiles, and correctly conclude that collectively the subpopulation exhibits significant evidence of a treatment effect. Additionally, the fact that TESS executes its optimization iteratively, unlike the greedy search of tree-based methods, enables it to rectify initial choices of subset that are later determined to be inferior.
Our second set of results considers treatment effects that do not align with the mean shift assumption that pervades the literature. Therefore, the null distribution is still $f_0 = N(0,1)$; however, the alternative is a mixture distribution $f_1= \frac{1}{2}N(-\delta,1)+\frac{1}{2}N(\delta,1)$. Here $\delta$ still captures the magnitude of the signal (treatment effect), and the remainder of the simulation process remains unchanged. This mixture distribution alternative, however, changes the detection task dramatically: while the average treatment effect is zero, there is still a clear difference in the outcome distribution between treated and control individuals.
Figure (ref) shows how each method's detection power and accuracy change as we vary each of three parameters that contribute to the strength of the treatment effect. If we compare these simulations to those above with a mean shift, TESS exhibits a consistent pattern of high performance, while the performance of the competing methods is dramatically lower. The detection power results indicate that, for the competing methods, it is hard to distinguish even strong distributional changes from random chance, while the accuracy results indicate that their pinpointing of the affected subpopulation is little better than random guessing. Given that there is no observable mean shift in these simulations, these results are consistent with what we expect: TESS is designed to identify more general distributional changes, while the other methods are unable to identify distributional changes without corresponding mean shifts.
There appears to be a consensus in the literature that the presence of a teaching aide in a regular-size classroom has an insignificant effect on test scores word-star-1990,krueger-star-1999,folger-star-1989,stock_watson-econ-2nd. (One significant effect was observed in first grade, but this effect was largely considered to be a false positive.) Therefore, we want to use TESS to compare regular classrooms with an aide to regular classrooms without an aide, to determine if there appears to be a subpopulation that was significantly and positively affected by the treatment. To do so, we replicated the analysis of the internal STAR team, using TESS to extend the results, with the goal of demonstrating what the STAR team could have surmised with present-day tools for uncovering heterogeneity. We replicate the original STAR analysis from word-star-1990,stock_watson-econ-2nd which includes the sum of the Stanford math and reading scores as the outcome of interest. For the data provided to TESS for detection, we combine the panel data across years and include student's grade level as a covariate.
We would also like to obtain an unbiased estimate of the average treatment effect in the subpopulation identified by TESS. Therefore, we follow a cross-validation paradigm, where the entire dataset is partitioned into ten folds, and iteratively each fold is held out as a validation set (to obtain an estimate of the treatment effect) while the remaining nine folds are provided to TESS (for detection). We further partition the data into records corresponding to students observed in a regular classroom with an aide and a regular classroom without an aide, which serve as treatment and control groups respectively. In three of the ten folds, TESS identified exactly the same subset, which we will call the “detected subpopulation”. Essentially, this detected subpopulation is composed of students in second or third grade, who attended an inner-city or urban school, receiving instruction from a teacher with 10 or more years of experience\footnote{The detected subpopulation excluded teacher experience between 25 and 30 years. Including this range yields qualitatively the same results and conclusions.}. Therefore, it appears that the presence of an aide raised the test-scores of students exhibiting the selected covariate values described above for grade, school type, and teacher experience, in addition to any values for gender, free-lunch status, teacher ethnicity, and teacher degree. The subpopulations that were returned in each of the ten folds exhibited a large amount of agreement with the detected subpopulation: the fold subpopulations exhibited 88% agreement (on average) with the detected subpopulation on the detection status of a record. The estimated average treatment effect for this detected subpopulation, averaged across all validation folds, is approximately a 34.19 point increase in total test score (36.45 and 22.28 for second and third grades respectively).
Given this consistency across folds, we use the full data to better understand the effect in the detected subpopulation generally. Table (ref) shows the evaluation of the treatment effect for all second-grade students (column 1), second-grade students in the detected subpopulation (column 2), and second-grade students in the complement of the detected subpopulation (column 3). Additionally, Figure (ref) shows the kernel density plots of the cumulative scores for all second-grade students and students in the detected subpopulation respectively. Figure (ref) depicts a strong similarity in the distribution of all second graders' scores with and without a full-time aide; there is a slight difference around the center of the distribution, but its magnitude is not sufficiently large to be significant, as seen by column 1 of Table (ref). Conversely, Figure (ref) depicts a difference in test scores for the detected subpopulation of second graders: there appears to be a clear effect of the treatment (dominated by a large mean shift), supported by column 2 of Table (ref). We conduct a similar analysis with third graders, and observe similar results in Figure (ref) and Table (ref). However, the effect of the treatment in third grade appears to result in less of a mean shift, and is better characterized by a change in the skew (third moment) and therefore, the overall form of the distribution (Figure (ref)). We note that because TESS is able to identify effects that change the distribution (and therefore higher order moments) of test scores, even if the difference in mean score between treatment and control students in third grade was smaller, TESS could potentially still identify the existence of a treatment effect.
There appears to be another consensus in the literature that small classrooms have a consistent, positive, and significant effect word-star-1990,krueger-star-1999,folger-star-1989; therefore, we also compare small classrooms to regular classrooms, and determine whether there appears to be a subpopulation which is the main driver of this effect. We conduct an analysis as above but with STAR data records corresponding to students observed in a small classroom (treatment group) and a regular classroom (control group). For this analysis, TESS identified the entire population, which is congruent with the previous literature's analysis of the consistent and significant average treatment effect in each grade. This result from TESS appears to indicate that the effect of small classroom size was not limited to a specific subpopulation. For both TESS analyses, we also conducted permutation testing to compensate for multiple hypothesis testing. Based on these results, we conclude that there is a less than 0.01% chance we would obtain a subpopulation with a score as extreme under the null hypothesis.
The detected subpopulation in the classrooms with aides is not only statistically significant, but may also provide domain insight into the efficacy of full-time aides. A possible explanation for the effect we observe in the detected subpopulation is the fact that 13 schools were chosen at random to have teachers participate in an in-service training session, which the literature has also deemed ineffective word-star-1990. More specifically, 57 teachers were selected each summer from these schools to participate in a three-day in-service to help them teach more effectively in whatever class type they were assigned to; part of the instruction focused on how to work with an aide and also had the aides present. We note that the in-service only occurred during the summers prior to 2nd and 3rd grade, which are the grades identified by TESS. Therefore, it is possible that when provided proper training, the combination of an aide and an experienced teacher can provide a significantly enhanced education environment even in the challenging teaching environments that exist in inner-city and urban schools. An additional explanation is that the educational benefits may be cumulative--i.e., in each additional year a student in this subpopulation has access to the combination of an aide and experienced teacher, the treatment effect compounds--similar to what has been demonstrated in small classrooms for the overall population krueger-star-1999. However, unlike in small classrooms, for this subpopulation in regular classrooms with an aide, the effects were not large enough to be distinguishable from zero (given the much smaller sample size of the affected subpopulation and smaller treatment effect) until after two years. While a more detailed follow-up analysis of these hypotheses might reveal other causal mechanisms at work, we believe that these results do present evidence that a treatment previously believed to be ineffective may actually have been effective for a particularly vulnerable subpopulation. Therefore, this analysis provides a sense of how TESS can be used as a tool for data-driven hypothesis generation in real-world policy analysis.
This paper has presented several contributions to the literature on statistical machine learning approaches for heterogeneous quantile treatment effects. Specifically, we detect the existence of a subpopulation for which the conditional quantile treatment effect (CQTE) is non-zero. This allows detection of treatment effects that manifest as arbitrary effects on the potential outcome distributions (or specific quantiles), rather than being limited to detection of mean shifts. Furthermore, we consider the challenge of identifying whether any subpopulation has been affected by treatment, and precisely characterizing the affected subpopulation, as opposed to the more typical problem setting of estimating individual-level treatment effects. We formalize the identification of subpopulations with significant treatment effects as an anomalous pattern detection problem, and present the Treatment Effect Subset Scan (TESS) algorithm, which serves as a computationally efficient test statistic for the maximization of CQTE over all subpopulations. We demonstrate that the estimator used by TESS satisfies the linear-time subset scanning property, allowing it to be efficiently and exactly optimized over subsets of a covariate's values, while evaluating only a linear rather than exponential number of subsets. This efficient conditional optimization step is incorporated into an iterative procedure which jointly maximizes over subsets of values for each covariate in the data: the result is a subpopulation, described as a subset of values for each covariate, which demonstrates the most evidence for a statistically significant treatment effect. In addition to its computational efficiency, we derive desirable statistical properties for the TESS estimator: bounded asymptotic probability of Type I and Type II errors under the sharp null hypothesis of no treatment effect, as well as providing sufficient conditions under the alternative hypothesis that will result in TESS exactly identifying the affected subpopulation. These properties apply more generally to the class of nonparametric scan statistics upon which TESS is built; therefore, this theory also provides additional contributions to the anomalous pattern detection, scan statistics, and goodness-of-fit literatures.
In addition to proposing a novel algorithm with desirable properties, we provide an extensive comparison between TESS and other recently proposed statistical machine learning methods for heterogeneous treatment effects (Causal Tree, Interaction Tree, and Causal Forest) through semi-synthetic simulations. Our results indicate that TESS consistently outperforms the other methods in its ability to identify and precisely characterize subpopulations which exhibit treatment effects. TESS significantly outperforms competing methods in the challenging scenarios where the treatment effect signal is weak (i.e., the signal magnitude is low or the affected subpopulation is small) because the subset scanning approach allows it to combine subtle signals across various dimensions of data in order to identify effects of interest. Moreover, TESS's detection performance is consistent even when the treatment outcome distribution in the affected subpopulation has the same mean as the control outcome distribution, while the competing methods demonstrate essentially no ability to identify the affected subpopulation in the absence of a mean shift.
After demonstrating TESS's performance through simulation, we explore the well-known Tennessee STAR experiment, searching for previously unidentified subpopulations with significant treatment effects. As a result of this analysis, TESS uncovered an intuitive subpopulation that seems to have experienced extremely significant improved test scores as a result of having a teacher's aide in the classroom, a treatment that has consistently been considered ineffective (as measured by the average treatment effect) by the literature on the Tennessee STAR. This provides a sense of how TESS can be utilized as a tool for generating hypotheses to be further explored and tested. We do however caution researchers to view algorithms like TESS not as a replacement, but rather an assistive tool, for developing scientific and behavioral theory. Results discovered by these methods should be investigated further and evaluated to develop a deeper theoretical understanding of the phenomena they uncover. When used to this end, these tools fill a critical void: in many contexts it is rare to know a priori which hypotheses are relevant and supported by data, and the use of traditional methods (e.g., regression) puts the onus on the researcher to know which hypothesis to test. This process necessitates that theory comes first, and subsequent investigation is a form of confirmatory analysis. However, such a process can become an impediment to data-driven discovery: there is an increasing need for scalable methods to use (big) data to generate new hypotheses, rather than just confirming pre-existing beliefs.
In the late 1970s, John W. Tukey began to outline his vision for the future of statistics, which included a symbiotic relationship between exploratory and confirmatory data analysis. He argues these two forms of data analysis “can--and should--proceed side by side” tukey-eda-1977 because he believed ideas “come from previous exploration more often than from lightning strokes” tukey-eda_cda-1980. To this end Tukey advocates for using data to suggest hypotheses to test, or what we now call data-driven hypothesis generation. We see our work as the natural evolution of Tukey's vision of data analysis: we develop an approach--rigorously conducted and theoretically grounded--to conduct exploratory analysis in randomized experiments, with the hope of catalyzing “lightning strokes” of discovery and the advancement of science.