{'fraction_non_alphanumeric': 0.032870970410611364, 'fraction_numerical': 0.009541405840664038, 'mean_word_length': 4.854698094930578, 'pattern_counts': {'":': 1, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Assessing effect heterogeneity of a randomized treatment using conditional inference trees

Treatment effect heterogeneity occurs when individual characteristics influence the effect of a treatment. We propose a novel, intuitive approach that combines prognostic score matching and conditional inference trees to characterize effect heterogeneity of a randomized treatment. One key feature that distinguishes our method from alternative approaches is that it controls the Type I error rate, i.e., the probability of identifying effect heterogeneity if none exists. This feature makes our technique particularly appealing in the context of clinical trials, where there may be significant costs associated with erroneously declaring that effects differ across population subgroups.


Introduction
Under mild assumptions, randomized experiments estimate the average causal effect (ACE) of an intervention. However, individuals may vary in their response to intervention so that the ACE is a poor representation of some people's expected benefit (or harm) from the intervention, a phenomenon often referred to as treatment effect heterogeneity [1,2,3]. Characterizing treatment effect heterogeneity can reveal both "weak responder" and "strong responder" subpopulations, leading to greater tailoring of intervention strategies.
For a given intervention, two heterogeneity-related questions arise: 1) Is the effect of the intervention heterogeneous? 2) If the intervention effect is heterogeneous, how does it vary across individuals? Most approaches address both questions using a single model or procedure. Traditional approaches to characterizing treatment effect heterogeneity have been primarily centered around regression modeling with interaction terms between the treatment and covariates. In such models, the interaction term can be used to assess whether treatment effect heterogeneity exists and also to describe heterogeneity. Alternatively, formal nonparametric tests have been developed to test a null hypothesis of zero average treatment effect for any subpopulation defined by covariates and whether the average treatment effect is identical for all subpopulations [4]. More recently, intense interest in personalized medicine combined with a rapidly growing toolkit of flexible machine learning techniques has produced several data-driven methods for characterizing treatment effect heterogeneity. Athey and Imbens [5] proposed the Causal Tree approach, which modifies the traditional CART, to determine subgroups with above-average or below-average treatment effects. The Causal Tree uses MSE-based criterion expressions to build (both splitting and cross validation) the tree and estimate the conditional average treatment effect at each leaf node over covariate defined subpopulations. A key feature of the Causal Tree approach is that it partitions the data into two subsets, one of which is used for subgroup construction and the other of which is used to estimate treatment effects within these subgroups. Wager and Athey [6] extended Causal Trees to Causal Forests, which averages treatment effect estimates over many Causal Trees. Tran and Zheleva [7] developed a method based on the Causal Tree framework to identify thresholds over ordinal (or continuous) treatments that trigger an effect. [8] discussed and proposed multiple treatment effect estimation methods that specifically handle high dimensional and observational datasets. [9] described several meta-algorithms for the estimation of conditional average treatment effects and introduced a new meta-algorithm in a binary treatment setting that is particularly effective when there are a greater number of units in one treatment group as opposed to another.
Within a Bayesian framework, Green and Kern [10] and Hill [11] proposed methods that apply Bayesian Adaptive Regression Trees (BART). Anoke et al. [12] offered a review of three classes of modeling approaches, relevant representative methods for each class (BART, propensity scores estimated with generalized boosted models (GBM) and the facilitating score (FS)) and evaluated these methods' ability to detect treatment effect heterogeneity in an exploratory manner.
While many of the aforementioned techniques have shown impressive abilities to identify heterogeneous subgroups in situations where heterogeneity exists, they are often overly aggressive in identifying treatment effect heterogeneity when, in truth, there is none. Put more simply, most existing procedures for detecting treatment effect heterogeneity do not control Type I error. This lack of control of Type I error is problematic in the context of randomized trials, where false declarations of treatment effect heterogeneity for a therapeutic agent could lead to wasteful follow-up studies and/or inappropriate off-label use in specific subpopulations. Two recently-proposed approaches have taken initial steps to address this issue: [13] proposed two new splitting criterions within a CART-like framework to maintain a balance between minimizing the error in estimating the treat-ment effect and maximizing heterogeneity, and [14] introduced a matching plus classification and regression tree (mCART) that reduces the potential for falsely detecting treatment effect heterogenity. In this paper, we build on some of these ideas to propose a novel approach to characterizing heterogeneity of a binary treatment effect on a continuous outcome, while explicitly controlling the Type I error rate. Our Treatment Effect Heterogeneity Tree (TEHTree) method involves building a conditional inference tree [15] using pairs of individuals matched on the prognostic score [16]. After describing the TEHTree method and providing theoretical motivation for matching based on prognostic scores, we present the results of a substantial simulation study demonstrating TEHTree's power under both null and alternative hypotheses and comparing its performance to the recentlydeveloped Causal Tree technique [5]. We also offer a comparison of the real-world performance of the two methods using data from a recent randomized trial.

Setup and Notation
Let Y = {Y 1 , . . . , Y n } be a continuous response vector for k subjects randomized to treatment Z = 0 and n − k subjects randomized to treatment Z = 1. The treatment assignments for all subjects are denoted by Z = {Z 1 , . . . , Z n }. An accompanying n × p-dimensional matrix X = {X 1 , . . . , X p } contains the p covariates for each of the n subjects with X i = {X i1 , . . . , X ip }.
In the counterfactual framework, each individual has a pair of counterfactual outcomes where Y 1i is the outcome of subject i if assigned to treatment Z = 1 and Y 0i the outcome if assigned to Z = 0. Hence, every individual has a (counterfactual, causal) treatment effect Y 1i − Y 0i , and the Average Treatment Effect (ATE) is defined as the mean of these within-individual differences, ATE = E[Y (1) − Y (0)] (we have switched to the parenthetical counterfactual notation Y (Z) to denote the observation Y Zi for arbitrary i). One of the benefits of randomization is that, under often plausible assumptions, the difference in means of randomized groups E(Y|Z = 1) − E(Y|Z = 0) estimates the ATE. The Stable Unit Treatment Value Assumption (SUTVA) [17,18] bundles together two assumptions: 1) the treatment assigned to an individual affects only the outcome for that individual, and 2) there is only one "version" of treatment. In many studies, SUTVA is plausible; a notable exception includes studies of infectious diseases in closed populations.
The other key assumption is ignorability, i.e., that treatment assignment Z is independent of the counterfactual pair (Y (0), Y (1)). While this "no unmeasured confounding" assumption is nontrivial in observational studies, it is satisfied by design in a randomized trial.
The counterfactual framework allows every individual i to experience a different effect of treatment, but because study participants are typically assigned to either Z = 1 or Z = 0, these individual-level effects Y 1i − Y 0i are unobserved. However, it is possible to estimate the Conditional Average Treatment Effect [19], Because Z is randomized, ignorability holds within any subset defined by X = x, and hence Therefore, the key challenge to characterizing treatment effect heterogeneity in randomized studies is to identify distinct subgroups defined by X with different CATEs. In the absence of treatment effect heterogeneity, the null hypothesis H 0 : CATE(x) = ATE ∀x holds. Most methods that seek to characterize how CATEs vary employ flexible semi-and non-parametric techniques in an attempt to identify regions of heterogeneity, and do not control the Type I error probability. In contrast, our approach embeds a classical parametric regression framework within a flexible tree model, allowing for both explicit control of the Type I error rate and characterization of treatment effect heterogeneity when it exists. The following two sections introduce the matching and conditional inference tree techniques that form the basis of our method.

Matching
If we observed Y 1i , Y 0i , and X i for all i then standard approaches to characterizing variability in a continuous outcome with respect to covariates could be used to estimate CATEs; for example, we could fit a regression tree using the differences Y 1i −Y 0i as outcomes and X i as predictors. However, in most trials an individual's outcome is observed under only one treatment, and hence Y 1i − Y 0i is unobserved. So, we propose to impute it by matching each individual i assigned to Z i = 1 with a "similar" individual j having Z j = 0 and using (Y i − Y j ) to approximate Y 1i − Y 0i . If j is an "exact" match for i in the sense that X i = X j , then we can use When the number of covariates p is even moderately large and/or elements of X are continuous, it will typically be impossible to find exact matches for most individuals. One way of overcoming this problem is by deriving a single measure that characterizes the "distance" between individuals.
If two individuals i and j with Z i = 1 and Z j = 0 have distance d ij = d between them, Hence, pairs matched according to d ij can be used to estimate CATEs provided ∆ 0 is small.
Note that, in general, ∆ 0 may be non-zero even if d = 0; it is the price paid for reducing the multidimensional vectors x and x to the scalar distance d.
A number of distance measures for matching have been proposed, some of which we review briefly here. Broadly speaking, these measures can be broken down into three categories according to how they define "similarity": based on the distance between covariate vectors, based on the probability of being treated (propensity score), and based on the predicted value of the outcome (prognostic score).
Covariate-based distances such as the Mahalanobis distance are vulnerable to the dimensionality reduction penalty alluded to above; when the number of covariates is large, ∆ 0 (x, x , d = ) may be relatively large even if the scalar distance is small (or zero) [20]. The propensity score [21,22], defined as the probability of receiving the treatment given the observed covariates, is often used for matching. However, when treatment is randomized, the propensity score is unhelpful since by design the covariates are independent of treatment assignment and as a result propensity score matching does not make ∆ 0 small.
Another way of summarizing the similarity between individuals is via the prognostic score [16].
If φ(X) is sufficient for Y 0 , in the sense that Y 0 is independent of X given φ(X), we call φ(X) a prognostic score. Individuals with similar prognostic scores have similar predicted values of the outcome under treatment Z = 0 (typically a control condition). If Y|X follows a homoscedastic linear model, then E(Y 0 |X) is a prognostic score, which can be estimated by fitting a regression of Y on X among those with Z = 0, then using that model to obtain predictions of the outcome under the control condition for all individuals [22]. Matching on prognostic scores is appealing in our context where the goal is approximate individual causal treatment effects Y 1i − Y 0i . As noted above, when Y 0 |X follows a homoscedastic linear model, E(Y 0 |X) ≡ φ(X) is a prognostic score. If two individuals i and j with Z i = 1 and Z j = 0 have the same prognostic score φ( This result immediately implies that if d ij = |φ(X i ) − φ(X j )| = 0, then ∆ 0 = 0 in (1). In other words, matching on the prognostic score yields pairs that can be used to estimate conditional average treatment effects. Note that this result holds if X is replaced by any measurable function m(X), so that if m captures the way in which X modifies the effect of treatment, then pairs matched on φ(X) retain information about effect modification.

Conditional Inference Trees
With matched pairs in hand that can be used to estimate CATE(x), the next step is to characterize how the CATE varies with x. Our approach uses conditional inference trees, a variant of decision trees which we briefly introduce here.
Decision trees [23] are non-parametric models based on recursively partitioning a sample into distinct subgroups that share similar outcome values. Sample partitions (or "splits") are defined by thresholding covariate values. The most popular and commonly used technique for building decision trees, the Classification and Regression Tree (CART) technique, was introduced by [24]. CART has been applied in many areas including economics [25,26], political science [27,10] and public health [28,29]. Because of its interpretability and flexibility, CART has also been incorporated into several methods for assessing treatment effect heterogeneity, notably CausalTree [5]. One of CART's drawbacks is that, because it considers many possible thresholds on all possible variables when searching for an optimal split, it has a tendency to overfit the data on hand and produce overly complex models. This overfitting tendency can be controlled somewhat by "pruning" trees based on a complexity parameter. However, as we show in our simulation study, even pruned CARTs do not control the Type I error for effect heterogeneity.
One alternative to CART is the Conditional Inference Tree (CTree), proposed by [15]. The main difference between CTrees and CARTs is in the splitting process: in CTrees, the processes for choosing a variable to split on or to stop splitting (the "variable selection" step) and choosing an optimal splitting threshold for the selected variable (the "splitting" step) occur sequentially, while in CART they happen simultaneously. In the variable selection step of CTree, the decision of whether not to continue splitting is based on a test of the global null hypothesis H 0 : E(Y|X) = E(Y), which is tested by considering all marginal null hypotheses H m 0 : E(Y|X m ) = E(Y) for m = 1, . . . , p. In the simplest case, each H m 0 can be assessed by calculating the p-value for the slope term from a univariate linear regression model of Y on X m . More generally, this step can accommodate a wide variety of models and test statistics; even if a statistic's sampling distribution is unknown, permutation tests can be used to calculate p-values for each partial null hypothesis. In our method, we calculate p-values associated with the (fixed) slope term from univariate linear mixed models.
Since H 0 = p m=1 H m 0 , the global null H 0 is rejected if the minimum p-value for all of the partial null hypotheses is less than a pre-specified level of significance. Control of Type I error can be achieved by setting this level of significance using an appropriate multiplicity adjustment to account for the testing of the p partial null hypotheses (see Section 2.4.3). If the minimum p-value exceeds the threshold, the tree does not split the given subset further. Otherwise, the partial null hypothesis H s 0 that results in the smallest p-value will indicate the covariate X s that is most strongly associated with the outcome Y and the algorithm proceeds to the next step to determine how to optimally threshold X s .
Once covariate X s has been selected for splitting, the second step of the CTree algorithm is to find the threshold c that maximizes the discrepancy |E(Y|X s ≤ c) − E(Y|X s > c)|. Typically, the process is performed over a finite set of candidate splits c = {c 1 , . . . , c k } which may be restricted to ensure properties such as a minimum number of observations in the partitions defined by (X s ≤ c) and (X s > c). This two-phase splitting procedure is repeated on the resulting partitions until no more subsets are eligible for splitting.

Treatment Effect Heterogeneity Trees (TEHTrees)
We propose a two-stage approach to assessing treatment effect heterogeneity in randomized studies. In the first stage, prognostic scores are calculated and every treated subject is matched to a control subject (with replacement) based on the prognostic score. In the second stage, within-pair differences in the outcome along with the covariate values of the treated member of each pair are used as inputs to a conditional inference tree. The nodes of the fitted conditional inference tree identify subgroups across which the causal effect of treatment varies. The full algorithm is provided below; in the sections that follow, we provide details about its key steps.

Full TEHTrees algorithm
1. Separate the dataset into a training and holdout sample for the purpose of constructing a tree and for estimating treatment effects over the determined terminal nodes.
2. Fit a model to calculate prognostic scores φ using the available training data, and obtain estimated prognostic scoresφ for each individual in the sample. Model details for prognostic score estimation are provided in Section 2.4.2.

For each pair
pair can now be viewed as a single "pseudo-individual" represented by the scalar continuous outcome δ l and the covariate vector X l ≡ X i . 5. Use the pseudo-individual data (δ l , X l ) created in the previous step and the desired Type I error rate as inputs to create a Treatment Effect Heterogeneity Tree (TEHTree), as described in Section 2.4.3.
6. Estimate the treatment effect within each terminal node of the fitted TEHTree as described in Section 2.4.4.
The preceding algorithm assumes that a sufficient amount of data is available to create a holdout test set of sufficient size to accurately estimate treatment effects within subgroups defined by each terminal node. If the sample size is limited, steps 2-6 can be carried out on the entire dataset with an in-sample estimation approach for step 6 as described in Section 2.4.4.

Estimating prognostic scores
We apply the Super Learner [30] to estimate the prognostic score φ(X) ≡ E(Y 0 |X) = E(Y|Z = 0, X). The base learners in our application consist of the sample mean, a linear model (with and without interaction terms), a generalized additive model, a Random Forest, stepwise regression (with and without interaction terms), and "polymars" (multivariate adaptive polynomial spline regression) as base learners. The Super Learner estimates the performance of these base learners using cross-validation and develops a weighted average over a combination of these methods [30].

Testing partial null hypotheses
Because the outcome values δ i,j are derived from pairs formed by matching with replacement, inputs to the conditional inference tree are correlated and hence a standard univariate linear regressionbased approach to evaluating the partial null hypotheses {H 0 m : E(Y|X m ) = E(Y)} will produce invalid p-values. Instead, we tested partial null hypotheses by fitting univariate linear mixed models of the form: is a random intercept corresponding to the control subject in each pair. Similar models are used to determine the optimal splitting for selected covariate X s , replacing X im in (3) by To establish proof of concept for our method, we used the Bonferroni method to adjust the marginal hypothesis test p-values for multiple comparisons, which sets the significance threshold at α p for desired Type I error rate α. Other less conservative adjustment methods could also be applied.

Treatment effect estimation
A Treatment Effect Heterogeneity Tree assigns every observation to a terminal node in the tree; if the tree has more than one terminal node, heterogeneity in treatment effects results because the estimated treatment effect differs between these terminal nodes. The algorithm in Section 2.4.1 uses an out-of-sample (test/holdout set) approach to estimating heterogeneous treatment effects that parallels the one used to develop the CausalTree method. After having generated the TEHTree based on the training set, we determine the TEHTree terminal node that each individual in the holdout set belongs to. Each terminal node T consists of the union of two subsets T 1 = {i : We compute the treatment effect estimate as Note that, for a 1:1 randomized treatment, large discrepancies between |T 1 | and |T 0 | are unlikely, and hence the precision of ∆(T ) will be proportional to 1 If limited sample size precludes carving out a holdout set from the original data, a straightforward in-sample estimation approach can be used. Let T denote the set of matched pairs from the data belonging to each terminal node. Then, the in-sample treatment effect estimate is simply Estimation of the precision of the in-sample estimate of ∆(T ) is complicated (relative to the out-of-sample approach) by the need to consider the correlation between matched pairs. Due to the limited sample, we apply this in-sample approach to the data illustration in Section 4.

Implementation
The TEHTree method was developed using a modified conditional inference tree framework and implemented in R [31] using relevant functions in the partykit package [32]. Continuous variables were generated using the mvtnorm package [33], matching was conducted using the Matching package [34], and all linear mixed models were fit using the nlme package [35]. The Super Learner was used to estimate the prognostic score; base learners were the sample mean, a linear model (with and without interaction terms), a generalized additive model, a random forest, stepwise regression (with and without interaction terms), and "polymars" (multivariate adaptive polynomial spline regression). Code for implementing TEHTree can be found at https: //github.com/AshwiniKV/TEHTree.

Simulation Study
We conducted simulations to evaluate the TEHTrees and compare its performance to Causal Trees [5], implemented using the package causalTree [36] with the default settings. We evaluated the Type I error, power and other statistical properties of estimators (e.g., bias and MSE) for TEHTrees and Causal Trees over different data generating scenarios that are documented in the Supplementary Tables 1 and 2. Factors that were varied over the scenarios included sample size (N = 100, 200,   500, 1000 and 2000), number of covariates, type of covariates (binary and continuous), coefficients, and pairwise correlation among covariates. The treatment variable Z was generated such that N/2 subjects received treatment (Z = 1) and the N/2 subjects received control (Z = 0). All results described in this simulation study are based on 1,000 simulations per scenario. Continuous covariates were generated from multivariate normal distributions with mean zero, unit variance, and varying pairwise correlations. Binary covariates were generated as independent Binomial(1, 0.5). Continuous outcomes were generated as independent N (µ, 1) with linear predictors varying according to the scenarios described in the Supplementary Materials.

Type I Error
For a case when there is no treatment effect heterogeneity, a Type I error occurs when a tree produces more than one terminal node. We therefore generated data under two sets of scenarios where there was no treatment effect heterogeneity. In the first set of scenarios, outcomes were generated from a linear model with main effects for treatment and covariates (Model (M1) in Supplementary   Table 1). In these scenarios, a simple linear model including treatment and covariates correctly specifies the prognostic score; in other words, the SuperLearner ensemble used to estimate the prognostic score contains the correct model. Figure 1   does not contain the true model. Figure 1(b) displays the Type I error (in logarithmic scale) for these scenarios. The Type I error rate of TEHTrees tends to increase with sample size and is no longer below the desired 0.05 in every scenario. However, the Type I error rate using TEHTrees is still much lower than the Type I error rate using Causal Trees, and when the correct prognostic score model is used, the Type I error rate is controlled. We note that this is a particularly challenging scenario for an approach based on prognostic score matching; even modest misspecification of the prognostic score could markedly increase the proportion of matched pairs where one individual has X > 0 and the other has X ≤ 0, leading to the (erroneous) conclusion that treatment effects are heterogeneous in X.

Power
We performed additional simulations to characterize the performance of TEHTrees and Causal Trees under a number of different data generating scenarios where treatment effect heterogeneity is driven by binary (dichotomized) covariates only, and a mixture of binary and continuous covariates.
To avoid crediting trees with incorrect splits, we defined the power for detecting treatment effect heterogeneity as the probability that a tree produced a split on a variable having a non-zero interaction with treatment (i.e., one that is responsible for producing treatment effect heterogeneity).
Note that this definition of power is somewhat generous; in cases where the treatment effect interacts with several different covariates, the methods are credited with rejecting the null hypothesis if the tree splits on any of these covariates.  (Figure 2(a)) and by continuous X 1 (Figure 2(b)); the power of TEHTrees is modestly lower for small sample sizes when both the indicator and continuous value contribute to heterogeneity (Figure 2(c)). The power of TEHTrees is substantially lower across most sample sizes for scenarios where heterogeneity is driven by I(−0.5 < X 1 < 0.5) (Figure 2(d)) and sin(X 1 ) ( Figure 2(e)).

Figure 3 displays the power of TEHTrees and Causal Trees when heterogeneity in the treatment
effect is due to two variables, X 1 and X 2 . For these relevant models (Models M8-M10), we considered two definitions of power: splitting on either X 1 or X 2 , and splitting on both X 1 and X 2 . In general, both types of power are higher for Causal Trees than TEHTrees.     Table 1 and Table 2. Power is defined to be the probability of the tree making a split on any variable with treatment effect heterogeneity ("Power Any") or on all variables with treatment effect heterogeneity ("Power All").

Treatment Effect Estimation
Model M6 in Figure 5(b) where heterogeneity is due to I(−0.5 < X 1 < 0.5); this is the same scenario displayed in Figure 2(d) where TEHTrees lacked power.

Illustration
In this section, we evaluate the proposed TEHTree's ability to detect treatment effect heterogenity against the Causal Tree approach. Both these methods are illustrated by an application to a real world dataset from the Box Lunch Study (BLS). The Box Lunch Study [37,38]  inference-based framework (as opposed to CART's greedy search approach) for growing a tree structure leads to a reduced risk of overfitting.
Both Causal Trees and TEHTrees are data-driven approaches that seek to identify subgroups that represent above-average or below-average benefit from treatment. In other words, these approaches determine estimates of the conditional average treatment effect (CATE). The trees are utilized to identify groups of study subjects that exhibit significant differences in the caloric intake due to a portion size intervention as opposed to no intervention. Relevant to this analysis, the caloric intake through the day is utilised as an outcome of interest and covariates include hunger and wanting. In a different scenario, the association between caloric intake and the covariates hunger1 and wanting1 is removed by permuting along caloric intake while maintaining the correlation between hunger1 and wanting1. In this modified dataset, to distinguish Scenario II from Scenario I, the covariates are labeled as hunger2 and wanting2. In Figure 6(c) and Figure 6(d), the tree frameworks are applied to the modified dataset. Using the same settings (as for Figure 6(a)), Causal Tree continues to find splits at wanting2 even with the removal of any association between the caloric intake and the covariates. A split at wanting ≥ 6 results in a subgroup with a positive treatment effect. However, at the same settings of the TEHTree (as for the TEHTree displayed in Figure   6(b)), there are no splits and a single root node with the corresponding average treatment effect is displayed in Figure 6(d). This difference in results suggests the ability of TEHTree to control for the false discovery of subgroups.  : Two scenarios are formulated to demonstrate differences in estimated treatment effects and ability to detect treatment effect heterogeneity. These tree frameworks (Causal Trees and TEHTrees) in Scenarios (I and II) are illustrated by an application to the Box Lunch Study (BLS) dataset. In both scenarios, the outcome of interest is energy intake through the day. Scenario I uses covariates hunger1 and wanting1, where there is no prior knowledge about correlation between the covariates and the outcome of interest. In Scenario II, the covariates hunger2 and wanting2 are permuted together to maintain their correlation but to remove existing correlations with the outcome.

Discussion
Characterizing treatment effect heterogeneity is becoming a common target of secondary analyses of data from randomized controlled trials. As an alternative to methods that require that the nature of potential subgroups be pre-specified (e.g., via covariate interactions with treatment), several methods have been recently proposed to detect treatment effect heterogeneity in a more data-driven manner [5,6,8,9]. However, while most of these methods incorporate procedures for preventing overfitting, they do not offer any guarantees about Type I error, i.e., the probability of identifying heterogeneous subgroups in the absence of treatment effect heterogeneity. Particularly in the context of randomized trials, explicit control of Type I error may increase the willingness of researchers to apply treatment effect heterogeneity techniques. In this paper, we propose TEHTrees, a novel method that uses a conditional inference tree framework to characterize effect heterogeneity of a binary treatment while controlling the Type I error rate.
The ability to tailor intervention strategies to particular population subgroups is aided by the detection of treatment effect heterogeneity. More specifically, differences from the average treatment effect indicate patient subpopulations that receive benefit/harm from a given treatment allocation.
However, traditional approaches to treatment effect heterogeneity that identify such subgroups do not address the issue of preserving or explicitly controlling Type I error rate. Our approach offers the ability to explicitly control the Type I error rate and is able to preserve balance (i.e., between the treatment and control groups) in patient characteristics by matching patients over computed prognostic scores.
As shown in our simulation study, popular existing methods (e.g., CausalTree) yield much higher than nominal Type I error rates, with this error rate generally increasing with sample size. In contrast, TEHTree maintains the specified Type I error rate across a majority of scenarios; however, for prognostic score models that are misspecified, the error rate increases with sample size (as displayed in Figure 1(b)). In some scenarios, the power of TEHTrees to detect true heterogeneity is competitive with CausalTree; in other scenarios, TEHTrees has lower power, but these discrepancies arise in scenarios where CausalTree has very high Type I error rates, i.e., its power curve lies above that of TEHTrees for both null and alternative hypotheses.
In simulation experiments, the Type I error rate using TEHTrees was below 0.05 (the pre-specified significance level) in all the considered scenarios, while the Type I error using Causal Trees was at least 0.15 in every scenario and even greater than 0.9 in many cases. Though Causal Trees have slightly greater power when there are continuous covariates, compared to TEHTrees, the power when using TEHTrees is actually greater than the power when using Causal Trees with binary covariates. When there is 'known' treatment effect heterogeneity, Causal Trees also tend to grow larger trees with more terminal nodes, particularly with larger sample sizes. This makes it more difficult to infer the characteristics of groups that truly respond more (or less) to treatment when using Causal Tree compared to TEHTrees. Causal Trees displayed lower MSE than TEHTrees in multiple scenarios, which was most likely due to greater variability in split points for TEHTree structures. We conjecture that the variability in split points is larger for TEHTrees than with Causal Tree and this can be attributed to the bias inroduced by the matching estimator. Decreasing bias in the matching estimator, or using an alternative approach to estimating the outcomes that are used as inputs in the conditional inference tree of TEHTrees, may therefore improve estimation of treatment effects with TEHTrees when there are continuous covariates.
TEHTrees offer a flexible approach to detecting effect heterogeneity and allows for numerous mod- however, additional assumptions and steps in the algorithm would be required to achieve balance in the covariates. Other modifications can also include extensions to account for multiple treatment allocations as opposed to a binary treatment assignment.

Supplementary Material Simulation Scenarios
Continuous covariates were generated from multivariate normal distributions with mean zero, unit variance, and varying pairwise correlations. Binary covariates were generated as independent Binomial(1, 0.5). Continuous outcomes were generated as independent N (µ, 1) with linear predictor µ as defined below in Suppmentary Table 1. We set α to 0.8, θ to 0. .  Evaluating Type I Error Figure 1(a) (M1) (C1)-(C3) Correctly specified prognostic score Figure 1 Incorrectly specified prognostic score Evaluating Power Heterogeneity according to I(X 1 > 0) and Heterogeneity according to I(−0.5 < X 1 < 0.5) Heterogeneity according to sin(X 1 ) Heterogeneity according to I(X 1 > 0) and I(X 2 > 0) Figure 3 Heterogeneity according to X 1 and I(X 2 > 0)

Additional scenarios for Models M1 -M10
The results evaluated using Model M1 and Model M2 were displayed in Figures 1(a) Table 4 presents results for two additional vectors of coefficients.           Table 11: Power and MSE of the estimated average treatment effect using TEHTrees (TT) and Causal Trees (CT) when data are generated according to Model M8. Power is defined to be the probability of making a split on X 1 or X 1 ("Power Any") or on X 1 and X 2 ("Power All") γ 1 γ 2 N ρ Power Any Power All  # nodes  N  MSE  TT  CT TT CT TT  CT  TT CT  3  3 Table 12: Power and MSE of the estimated average treatment effect for TEHTrees (TT) and Causal Tree (CT) when data are generated according to Model M9. Power is defined to be the probability of making a split on X 1 or X 1 ("Power Any") or on X 1 and X 2 ("Power All").  Table 13: Power and MSE of the estimated average treatment effect for TEHTrees (TT) and Causal Tree (CT) when data are generated according to Model M10. Power is defined to be the probability of the tree making a split on any variables ("Power Any") or on all variables with treatment effect heterogeneity ("Power All"). γ is the coefficient of the interaction term.
γ Vars N ρ Power Any Power All # nodes MSE int.
TT CT  TT CT TT  CT  TT  CT  6