# Comparing Machine Learning Methods for Estimating Heterogeneous Treatment
Effects by Combining Data from Multiple Randomized Controlled Trials

Carly Lupton Brantner
Department of Biostatistics
Johns Hopkins Bloomberg School of Public Health
<EMAIL_ADDRESS>
&Trang Quynh Nguyen
Department of Mental Health
Johns Hopkins Bloomberg School of Public Health &Tengjie Tang
Department of Statistical Science
Duke University &Congwen Zhao
Department of Biostatistics and Bioinformatics
Duke University &Hwanhee Hong
Department of Biostatistics and Bioinformatics
Duke University &Elizabeth A. Stuart
Department of Biostatistics
Johns Hopkins Bloomberg School of Public Health

###### Abstract

Individualized treatment decisions can improve health outcomes, but using data
to make these decisions in a reliable, precise, and generalizable way is
challenging with a single dataset. Leveraging multiple randomized controlled
trials allows for the combination of datasets with unconfounded treatment
assignment to improve the power to estimate heterogeneous treatment effects.
This paper discusses several non-parametric approaches for estimating
heterogeneous treatment effects using data from multiple trials. We extend
single-study methods to a scenario with multiple trials and explore their
performance through a simulation study, with data generation scenarios that
have differing levels of cross-trial heterogeneity. The simulations
demonstrate that methods that directly allow for heterogeneity of the
treatment effect across trials perform better than methods that do not, and
that the choice of single-study method matters based on the functional form of
the treatment effect. Finally, we discuss which methods perform well in each
setting and then apply them to four randomized controlled trials to examine
effect heterogeneity of treatments for major depressive disorder.

_Keywords_ treatment effect heterogeneity, combining data, personalized
medicine, machine learning

##  1 Introduction

When tailoring treatment regimens to individual patients, one must strive to
understand how different treatment options might affect the specific patient
based on their characteristics or context. Rather than using a one-size-fits-
all approach, clinicians and researchers are turning more towards personalized
medicine with the goal of improving clinical outcomes. In this setting, the
focus of estimation becomes conditional average treatment effects, i.e., how
well the treatment is expected to work conditional on the person’s known
characteristics.

The benchmark for estimating treatment effects in an unbiased manner is most
often a randomized controlled trial (RCT). In an RCT, participants are
randomly assigned to treatment or control, therefore ensuring unconfounded
treatment assignment and unbiased treatment effect estimates in the given
sample. However, these trials often have sample sizes that are large enough to
detect main effects but lack power to estimate heterogeneous treatment effects
(Enderlein, 1988). To overcome these specific issues, researchers have started
combining information from multiple studies to improve treatment effect
estimation. Multiple studies allow for larger sample sizes and at times a more
representative sample of the target population. In the setting with multiple
RCTs, meta-analysis or hierarchical models are common techniques to combine
studies and estimate treatment effects (Debray et al., 2015; Seo et al.,
2021). These approaches often do not explicitly target conditional average
treatment effects though, and often only use aggregate-level data which makes
it challenging to estimate treatment effects conditional on individual-level
characteristics. Furthermore, meta-analysis is commonly applied within a
parametric framework, which is highly interpretable but requires
prespecification of effect moderators and distributional assumptions for
parameters. Non-parametric approaches are worth exploring in this setting
because they allow for high levels of flexibility in outcome and treatment
effect functions. Relationships between covariates and treatment effect can be
complex and non-linear in reality, and non-parametric machine learning methods
can better handle those scenarios.

Many non-parametric approaches exist to estimate heterogeneous treatment
effects (Künzel et al., 2019; Athey et al., 2019; Green and Kern, 2012;
Kennedy, 2020; Nie and Wager, 2021; Dandl et al., 2022); however, these
approaches have generally been developed only for the single-study setting.
Several of the common approaches are discussed in the section to follow (3.1),
and we subsequently extend these methods for use in multiple studies. Recent
research has investigated a few non-parametric approaches for the multiple
study setting, mostly geared towards combining data from one RCT with a large
observational dataset (Yang et al., 2020, 2022; Kallus et al., 2018; Rosenman
et al., 2020). In that work, the focus is often on estimating the bias present
in the observational data to determine the level at which the observational
study estimates can be combined with the RCT estimates. These methods are
therefore not as straightforward to use in the multiple RCT setting. With
multiple RCTs, each individual trial has the benefit of unconfounded treatment
assignment, but significant cross-trial heterogeneity could still exist due to
both observed and unobserved factors. The focus in this case is no longer de-
biasing one of the datasets, but instead determining the amount of
heterogeneity present and how to account for it.

Brantner and colleagues Brantner et al. (2023) wrote a comprehensive review of
methods geared towards combining datasets to estimate treatment effect
heterogeneity. That review included approaches for multiple RCTs; the most
common were individual participant-level data one-stage meta-analyses (Debray
et al., 2015). One alternative approach focuses on combining RCTs to estimate
conditional average treatment effects in a non-parametric framework (Tan et
al., 2021). That work by Tan and colleagues was done in the federated learning
setting, in which individual-level data could not be shared across study sites
and instead only aggregate results or models could be shared. In the sections
to follow, we tailor Tan et al.’s method to when individual-level data can be
shared across trials.

To our knowledge, this paper is the first to describe and compare machine
learning options for estimating heterogeneous treatment effects using data
from multiple RCTs, in the setting in which all data can be shared across
trials. Because not many methods exist to do this, we demonstrate several
options for extending current methods for single studies to the multiple-study
setting. We also build off of Tan et al.’s Tan et al. (2021) approach by
adapting it to the case when individual-level data can be shared across
trials. We conduct extensive simulations with varying data generation
parameters to determine which of the single-study and aggregation methods
perform best depending on different amounts of cross-trial heterogeneity in
the effects. We then apply the approaches to a set of four RCTs of depression
treatments and discuss the variability in estimates across the approaches and
potential substantive conclusions that can be made.

##  2 Notation

The estimand considered in this paper is the conditional average treatment
effect (CATE), defined under Rubin’s potential outcomes framework (Rubin,
1974). Let $A$ denote a binary treatment indicator (often treatment versus
control), $\boldsymbol{X}$ represent covariates, and $Y$ represent a
continuous outcome. Under Rubin’s framework, $Y(0)$ and $Y(1)$ denote the
potential outcomes under control and treatment, respectively. Let $S$ be a
categorical variable representing the trial in which the individual
participated and ranging from $1$ to $K$, where $K$ is the total number of
RCTs. Finally, represent the probability of receiving treatment given
covariates and trial membership (propensity score) as
$\pi_{s}(\boldsymbol{X})=P(A=1|\boldsymbol{X},S=s)$.

With a continuous outcome, the CATE is

$\tau(\boldsymbol{X})=E(Y(1)|\boldsymbol{X})-E(Y(0)|\boldsymbol{X}).$ (1)

In this paper, we note that the goal estimand is this “universal” CATE (1)
built off of potential outcomes that are not dependent upon study membership.
However, many methods in the following sections target a study-specific CATE:

$\tau_{s}(\boldsymbol{X})=E(Y(1)|\boldsymbol{X},S=s)-E(Y(0)|\boldsymbol{X},S=s).$
(2)

To identify the estimand when combining data across RCTs, many of the standard
causal inference assumptions are required, including the Stable Unit Treatment
Value Assumption (SUTVA) within each RCT. Other standard assumptions include:
unconfoundedness (Assumption 1), consistency (Assumption 2) and positivity
(Assumptions 3 and 4) (Tan et al., 2021). Assumption 2 varies slightly
depending on the estimand; under the universal CATE estimand (Equation 1), we
assume overall consistency, while under the study-specific estimand (Equation
2), we assume consistency within each study. Assumption 4, which requires that
any $\boldsymbol{X}$ is possible to be observed in all studies, can be relaxed
depending on the method.

###### Assumption 1

$\\{Y(0),Y(1)\\}\perp\\!\\!\\!\perp A\mid\boldsymbol{X},S=s$ for all studies
$s$.

###### Assumption 2

$Y=AY(1)+(1-A)Y(0)$ almost surely (in each study).

###### Assumption 3

There exists a constant $c>0$ such that $c<\pi_{s}(\boldsymbol{x})<1-c$ for
all studies $s$ and for all $\boldsymbol{x}$ values in each study.

###### Assumption 4

(Can be relaxed) There exists a constant $d>0$ such that
$d<P(S=s|\boldsymbol{X}=\boldsymbol{x})<1-d$ for all $\boldsymbol{x}$ and $s$.

##  3 Methods

This paper includes methods developed for treatment effect estimation in a
single study and aggregation approaches that apply these methods to multiple
studies. We chose some methods based off of those explored in Tan et al.’s Tan
et al. (2021) paper that have a similar goal to the present work. This section
discusses three single-study methods and five aggregation options that apply
the single-study methods to the multi-study setting.

###  3.1 Single-Study Methods

For a given RCT, many machine learning methods have been developed for CATE
estimation. The single-study methods that exist can be grouped into multiple
categories, as delineated by Brantner et al (2023). For ease of comparison,
three approaches are included that are common in practice, user-friendly, and
have been shown to be effective in previous literature: the S-learner,
X-learner (Künzel et al., 2019), and causal forest (Athey et al., 2019). The
first two approaches are multi-step procedures that involve first estimating
the conditional outcome mean under treatment or control and then combining the
two into one CATE function, while the causal forest involves tree-based
partitioning of the covariate space by treatment effect. In this paper, we use
random forests as the base learners for both the S-learner and the X-learner
to best compare with the causal forest, which is inherently forest-based.

####  3.1.1 S-Learner

The first single-study machine learning method used in this paper is called
the "S-learner" (Künzel et al., 2019). This method is classified as a "meta-
learner" in that it combines base learners (i.e., regression models) of any
form in a specific way (Künzel et al., 2019). The S-learner uses a base
learner (i.e., a random forest) to estimate a conditional outcome mean
function given observed covariates and assigned treatment:

$\mu(\boldsymbol{X},A)=E(Y|\boldsymbol{X},A).$

The conditional outcome mean function in this approach is not specific to
treatment group, but instead treatment is included together with the
covariates as features to be used by the random forest. The CATE can then be
directly estimated by plugging in $0$ and $1$ for the treatment indicator to
obtain predicted outcomes under treatment and control for each individual and
calculate

$\hat{\tau}(\boldsymbol{X})=\hat{\mu}(\boldsymbol{X},1)-\hat{\mu}(\boldsymbol{X},0).$

####  3.1.2 X-Learner

The second approach considered here is another meta-learner called the
"X-learner" (Künzel et al., 2019). The X-learner takes a similar approach as
the S-learner by modeling the conditional outcome mean functions before
estimating the CATE directly. However, rather than estimating one outcome mean
function for $Y(1)$ and $Y(0)$ simultaneously, the X-learner estimates two
functions separately and then imputes treatment effects for each treatment
group.

Specifically, the X-learner involves three steps. First, the conditional
outcome mean functions are estimated using base learners (in this case, random
forests) like in the S-learner, but separately by treatment group:

$\mu_{1}(\boldsymbol{X})=E(Y(1)|\boldsymbol{X})~{}~{}~{}\text{and}~{}~{}~{}\mu_{0}(\boldsymbol{X})=E(Y(0)|\boldsymbol{X}).$

Next, the unobserved potential outcomes for individuals in the treatment and
control groups are predicted using those models:

$\tilde{D}_{i:A=1}=Y_{i:A=1}-\hat{\mu}_{0}(\boldsymbol{X}_{i:A=1})~{}~{}~{}\text{and}~{}~{}~{}\tilde{D}_{i:A=0}=\hat{\mu}_{1}(\boldsymbol{X}_{i:A=0})-Y_{i:A=0}.$

Then $\tilde{D}$ is regressed on $\boldsymbol{X}$ to estimate
$\tau(\boldsymbol{X})$. This is done within each treatment group separately,
resulting in two estimates, labeled $\hat{\tau}_{1}(\boldsymbol{X})$ and
$\hat{\tau}_{0}(\boldsymbol{X})$. Finally, these are combined to obtain one
CATE estimate:

$\hat{\tau}(\boldsymbol{X})=g(\boldsymbol{X})\hat{\tau}_{1}(\boldsymbol{X})+(1-g(\boldsymbol{X}))\hat{\tau}_{0}(\boldsymbol{X}),$

where the weight $g(\boldsymbol{X})$ is often an estimate of the propensity
score (the case in this paper) or can be chosen otherwise (Künzel et al.,
2019).

####  3.1.3 Causal Forest

The third single-study approach is the causal forest (Athey et al., 2019). The
causal forest is similar to a random forest, but the focal estimand is the
treatment effect itself, rather than the outcome for a given individual. The
causal forest is based off of a causal tree, which involves recursive
partitioning of the covariates to best split based on treatment effect
heterogeneity. Here, the treatment effect is estimated as the difference in
average outcomes between the treatment and control group individuals within
leaves. From there, the causal forest is the weighted aggregation of many
causal trees.

One potential challenge with causal forests is that bias could occur when
there is overlap between the data used to form the trees and data used to
estimate the treatment effects within leaves. A solution to that problem,
called "honesty", has been proposed (Wager and Athey, 2018) .This concept
ensures that for every individual involved in creating the tree, their outcome
is used either for splitting the tree or estimating the treatment effect
within a leaf, but not both. Honesty has been used some in the literature, but
there is not a widespread conclusion as to whether trees should be fit with or
without honesty depending on the scenario. Dandl and colleagues compared
honesty versus adaptive (not honest) forests in their simulations including
causal forests and found that in their setting that was meant to represent an
RCT, the adaptive forests performed better (Dandl et al., 2022). Additionally,
honesty requires large sample sizes. Thus, we do not include honesty in the
causal forests in the primary simulations but do investigate it in a second
round of method comparisons.

###  3.2 Aggregation Methods

In many contexts, there are multiple RCTs available that compare the same two
treatments. It is then worth considering methods that allow combining across
trials. When aggregating to the multi-study level, the question becomes: how
much does the treatment effect vary based on study membership? This
variability can range along a continuum, where on one end is the possibility
that the trials are all very homogeneous in terms of the CATE, meaning that
participants in trial $j$ and in trial $k$ who have the same covariate values
would have the same treatment effect. At the other extreme, individuals with
the same covariates but in different trials could have completely different
treatment effects. These differences can be due to heterogeneity in the sites
in which the trials were conducted, heterogeneity in trial procedures
(including the treatment or control conditions themselves), heterogeneity in
trial samples, or other reasons. The aggregation methods to follow take
different approaches to incorporating trial membership into the treatment
effect estimation, ranging from assuming trial membership does not matter at
all, to allowing it to matter just as much as any other characteristic.

####  3.2.1 Complete Pooling

A complete pooling approach is very straightforward: the researcher simply
takes all data from each of the $K$ RCTs, creates a single dataset, and then
fits one of the three previously described methods (S-learner, X-learner, or
causal forest) to the pooled dataset. This approach is quick and easy to do,
but requires many assumptions. Namely, this approach assumes a high level of
homogeneity across trials and specifically that the CATE function is shared
across studies. This method is included because it represents a naive
comparison point and because it provides universal CATE estimates (i.e., not
study-specific).

####  3.2.2 Pooling with Trial Indicator

An alternative pooling approach is to incorporate trial membership in the
models but essentially still perform the pooling as before. Here, all of the
individual data from each RCT is combined into one comprehensive dataset, but
a categorical variable is included that represents the trial in which the
individual participated. Then, the researcher can apply one of the single-
study approaches to this full dataset, allowing for all of the covariates,
including trial membership, to be involved in the treatment effect function.
In this way, if trial membership is important for estimating effects,
estimates should be somewhat informed by trial membership; otherwise, the
treatment effect estimates should be similar across trials. While the previous
complete pooling approach gives estimates that were not trial-specific, this
approach yields trial-specific CATE estimates.

####  3.2.3 Ensemble Approach

The next approach is based off of Tan and colleagues’ (2021) methods for
federated learning, originally developed for scenarios in which individual
data cannot be shared across trial sites. Their original approach fits trial-
specific models and then applies those models to data from a single
coordinating site to derive an ensemble. We propose an adaptation of Tan’s
approach for settings where individual-level data from all trials are
available to the analyst.

This adaptation of Tan et al.’s approach involves three steps.

  1. 1.

First, the researcher builds localized models for the CATE within each trial,
using one of the three single-study methods previously discussed (S-learner,
X-learner, or causal forest).

  2. 2.

Next, they apply these localized models to each individual across all of the
RCTs to get for each individual their trial-specific CATE estimates, i.e., the
predicted effects had the individual been part of study 1, study 2 and so on.
For $K$ studies with a total of $N$ individuals in all studies combined, there
will be $K$ trial-specific CATE models. Once each of these models are applied
to all $N$ data points, every individual will have $K$ different estimates of
their CATE. So there will ultimately be $N*K$ CATE estimates in what Tan et
al. define as an "augmented" dataset. The difference between the second step
here and what Tan et al. did is that we apply the study-specific models to all
data points in all trials, rather than having to restrict to a single
coordinating site.

  3. 3.

The third and final step is to fit an ensemble model to the augmented dataset
that has CATE estimates for every individual crossed with every trial. In this
ensemble model, the response variable is the CATE estimate, and the predictors
are the individual covariates and a categorical variable indicating the local
model that had been used to compute the CATE estimate. We use three different
options for this final ensemble model fit to the augmented dataset: a
regression tree, a random forest, and a lasso regression. The regression tree
and random forest were explored in Tan et al.’s paper (2021), while the lasso
regression was added to provide a parametric comparison point.

The resulting functions from these ensemble approaches are trial-specific
estimates of the CATE; however, they have been adapted based on the CATEs from
the other trials. Therefore, this method allows for trial heterogeneity but
incorporates information across trials to hopefully improve the model from
each trial.

####  3.2.4 IPD Meta-Analysis

As a comparison point in the simulations to follow, we also include an
individual patient-level data (IPD) meta-analysis with a random intercept for
trial membership. This method serves as a parametric comparison to the
primarily non-parametric approaches outlined above. A meta-analysis does not
employ a single-study method like the S-learner, X-learner, or causal forest;
instead, all of the data is pooled together and trial-level relationships can
be included as fixed or random effects. The meta-analysis model can be
parametrized in many different ways; in this paper we set it up to mostly
mimic the setup of the first scenario in the simulation to follow:

$Y=(\alpha_{0}+a_{s})+(\alpha_{1}+b_{s})X_{1}+\alpha_{2}X_{2}+\alpha_{3}X_{3}+\alpha_{4}X_{4}+(\zeta+z_{s})A+(\theta+t_{s})X_{1}A+\epsilon.$

In this model, we allow the intercept to include a fixed component
($\alpha_{0}$) and a random component by study ($a_{s}\sim
N(0,\sigma_{a}^{2}))$, and our residual error is $\epsilon\sim
N(0,\sigma^{2})$. The fixed effects are
$\boldsymbol{\alpha}=\\{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\\}$, the
coefficients relating the covariates to the outcome; $\zeta$, the coefficient
for treatment; and $\theta$, the coefficient of the interaction between
treatment and a moderator $X_{1}$ (Seo et al., 2021). The random effects by
study are $b_{s}\sim N(0,\sigma_{b}^{2})$, the random slope for the covariate
$X_{1}$; $z_{s}\sim N(0,\sigma_{z}^{2})$, the random slope for treatment; and
$t_{s}\sim N(0,\sigma_{t}^{2})$, the random slope for the treatment-$X_{1}$
interaction term. From here, the estimate of the conditional average treatment
effect can be calculated as
$\hat{\tau}_{s}(\boldsymbol{X})=(\hat{\zeta}+\hat{z}_{s})+(\hat{\theta}+\hat{t}_{s})X_{1}$.

The meta-analysis framework assumes that the CATE function is shared across
studies, but that the mean potential outcome under control can differ across
studies. Notably, this functional form of the CATE assumes linear
relationships, and one must prespecify all variables that might be relevant to
the main effect of the covariates and to the treatment effect.

##  4 Simulation Setup

To compare both the single-study and aggregation methods, we performed a
simulation study, simulating data from multiple randomized controlled trials
and changing parameter values to compare which methods achieve the lowest mean
squared error (MSE) between the estimated and true individual CATEs. Because
there were three single-study methods and five aggregation methods being
compared along with meta-analysis, there were $3*5+1=16$ total combinations of
methods applied to each simulated dataset.

###  4.1 Data Generating Mechanism

The data generating mechanism is based somewhat off of Tan et al. Tan et al.
(2021) and Kunzel et al. Künzel et al. (2019) since both used methods similar
to those compared here. The potential outcomes are generated using the
following model:

$Y_{i}(a)=m(\boldsymbol{x}_{i},s_{i})+\frac{2a-1}{2}*\tau(\boldsymbol{x}_{i},s_{i})+\epsilon_{i}$
(3)

where $m(\boldsymbol{x}_{i},s_{i})$ represents the outcome mean conditional on
covariates and trial, and $\tau(\boldsymbol{x}_{i},s_{i})$ is the CATE. In the
main setting for the data generation, we employed two options for $m$ and
$\tau$. The first setup (1a) involves a linear $m$ and piecewise linear
$\tau$:

$m(\boldsymbol{x},s)=x_{1}/2+\sum_{j=2}^{4}x_{j}+\beta_{s}+\delta_{s}*x_{1}\text{
and }\tau(\boldsymbol{x},s)=x_{1}*I(x_{1}>0)+\beta_{s}+\delta_{s}*x_{1}.$

The second setup (1b) involves a more complicated non-linear function for
$\tau$:

$m(\boldsymbol{x},s)=0\text{ and
}\tau(\boldsymbol{x},s)=g(x_{1})g(x_{2})+\beta_{s}+\delta_{s}*x_{1}$

where $g(x)=\frac{2}{1+\exp(-12(x-1/2))}$.(Künzel et al., 2019) In both of
these, the coefficients $\beta_{s}$ represent trial-specific main effect
coefficients, and $\delta_{s}$ represent trial-specific interaction effect
coefficients (interaction between trial and the moderator $x_{1}$). In both
setups, $x_{1}$ is an effect moderator, and in the second setup, $x_{2}$ is as
well. If the coefficients $\beta_{s}$ and $\delta_{s}$ differ across $s$
(i.e., trial membership), then trial is making an impact in the moderation.

From this information, the components simulated are listed as follows:

  1. 1.

For each simulation, the number of trials was $K=10$.

  2. 2.

Each trial had a sample size of 500 individuals.

  3. 3.

Within each trial, we simulated five continuous covariates per person,
$\boldsymbol{X}_{i}\sim N(0,I_{5})$, with each following a standard normal
distribution.

  4. 4.

Each person was then assigned a treatment status, $0$ or $1$, according to a
propensity score of $\pi_{i}=0.5$ within each trial.

  5. 5.

Each person was also assigned an error term for their outcome function, so
$\epsilon_{i}\sim N(0,.01)$.

  6. 6.

We then sampled trial-specific main effect and interaction effect terms. Each
of the $K=10$ studies was assigned a main effect term according to
$\beta_{s}\sim N(0,\sigma_{\beta}^{2})$ and an interaction effect term
according to $\delta_{s}\sim N(0,\sigma_{\delta}^{2})$. The values of the
standard deviations were:
$(\sigma_{\beta},\sigma_{\delta})\in\\{(0.5,0),(1,0),(1,0.5),(1,1),(3,1)\\}$.
These standard deviation pairs are defined as (Low-Low), (Med-Low), (Med-Med),
(Med-High), (High-High), respectively.

  7. 7.

From this information, $m$, $\tau$, and $Y$ were calculated under either of
the two setups described above (1a and 1b).

From the above setups, there were five standard deviation pairs for the trial
effects, and two functional forms for $m$ and $\tau$, therefore yielding 10
total scenarios.

We then included one other scenario (2) to see how the methods would perform
when the functional form of the CATE itself differed across trials - a
particularly challenging situation for pooling. For this scenario, we used the
same form for $Y_{i}$ as in Equation (3), and now we set $m$ and $\tau$ to be
such that $m$ is linear and $\tau$ depends on study:

$m(\boldsymbol{x},s)=x_{1}/2+\sum_{j=2}^{4}x_{j}$

and

$\tau(\boldsymbol{x},s)=I(s\in\\{1,2,3,4\\})*g(x_{1})g(x_{2})+I(s\in\\{5,6,7,8\\})*x_{1}*I(x_{1}>0)+I(s\in\\{9,10\\})*0$

where $g(x)$ is as previously defined.

For each simulation setup, we generated 1,000 simulated datasets. Necessary
packages included causalToolbox for the S-learner and X-learner (Künzel et
al., 2019), grf for the causal forest(Athey et al., 2019), rpart for the
ensemble tree (Therneau et al., 2015), ranger for the ensemble forest (Wright
and Ziegler, 2015), glmnet for the ensemble lasso (Friedman et al., 2017), and
lme4 for the mixed effects meta-analysis (Bates, 2010). Some functions were
based off of those in the ifedtree package (Tan et al., 2021) but were adapted
to the setting in which data could be shared across trials. For each method
and each iteration, performance of the different approaches was assessed based
on the mean squared error (MSE) between the true individual CATE estimates and
the estimated individual CATE estimates, and these MSEs were ultimately
averaged across the 1,000 repetitions.

##  5 Simulation Results

The following tables and figures display the performance results across 1,000
iterations of each parameter combination/scenario. Table A2 in the Appendix
shows the average and standard deviation of the MSEs across all iterations of
a given scenario and approach. Figure 1 shows these average MSEs for every
approach for each scenario, broken down by the standard deviations of the
trial main and interaction effects. In the piecewise linear and non-linear
CATE scenarios, as the trial coefficients (both main and interaction effects)
increase in variability, the MSE increases, meaning the methods estimate
individual CATEs more poorly. This is consistent with the idea that when trial
membership is involved in the treatment effect function, the CATEs vary across
trials and therefore are harder to estimate when data is pooled across
studies. Notably, this increase in MSE happens much more quickly for the
complete pooling approaches.

Figure 1: Average MSE for each parameter combination across all approaches.*
*The y-axis was cut off for ease of visualization; note that for the High-High SD combination, some methods were therefore missing from the graph because the MSEs were very high (complete pooling and meta-analysis). The SD of study main and study interaction coefficient pairs are as follows: Low-Low: 0.5, 0; Med-Low: 1, 0; Med-Med: 1, 0.5; Med-High: 1, 1; High-High: 3, 1. Facet labels refer to the simulation scenarios: 1a, 1b, and 2.

In the piecewise linear scenario (1a), the most consistently effective
approaches in terms of MSE are when the causal forest is used as the single-
study method and when the aggregation approach is either pooling with trial
indicator, ensemble forest, or ensemble lasso. The X-learner also performs
relatively well in terms of MSE. Meta-analysis performs well, which is
expected because the model was set up to mostly match the true functional form
of the CATE in this scenario. For the non-linear scenario (1b), the ensemble
lasso and meta-analysis perform notably worse, which makes sense due to the
complexity of the functional form of the CATE, as it includes the product of
two expit functions, and the lasso and meta-analysis assume a parametric
linear relationship between covariates and outcome. The ensemble forest and
pooling with trial indicator again estimate the CATEs well, with all single-
study methods performing similarly. While the S-learner was not very effective
with the piecewise linear CATE (1a), it was more effective with the non-linear
CATE (1b).

Figure 1 also displays the results for the variable CATE scenario (2). Here,
the causal forest is clearly performing the best of the three single-study
methods, while the S-learner is not performing as well. The most effective
aggregation methods are again pooling with trial indicator and ensemble
forests, along with ensemble trees. The least successful aggregation approach
is the ensemble lasso, followed by meta-analysis and complete pooling.

Figure 2 displays the average MSE across all data generation setups (i.e.,
piecewise linear, non-linear, and variable CATE) and all iterations.
Consistent with our scenario-specific findings, the complete pooling
approaches are ineffective at estimating the individual CATEs compared to the
pooling with study indicator, ensemble tree, ensemble forest, ensemble lasso,
and meta-analysis. Within those more effective aggregation approaches, the
three single-study options perform somewhat similarly but with the S-learner
doing the worst in almost all cases. The overall best approaches, in terms of
average MSE, are the causal forest with pooling with trial indicator, the
causal forest with an ensemble forest, and the X-learner with an ensemble
forest.

Figure 2: Average MSE across all scenarios and iterations.

To more formally examine these results, we regressed the average MSE across
iterations on the methods and parameter combinations, just within the
piecewise linear and non-linear CATE scenarios and excluding meta-analysis.
Specifically, the regression is such that:

$\displaystyle MSE=$
$\displaystyle\beta_{0}+\beta_{1}*singlestudy+\beta_{2}*aggregation+\beta_{3}*singlestudy*aggregation$
$\displaystyle+\beta_{4}*main_{sd}+\beta_{5}*interaction_{sd}+\beta_{6}*scenario+\epsilon.$

From this regression, there were no significant differences in performance
across single-study methods, but all aggregation methods performed
significantly better than complete pooling. The ensemble forest had the best
average MSE for the S-learner and X-learner, and pooling with trial indicator
had the best average MSE for the causal forest.

Finally, we also performed 500 more iterations of each scenario and parameter
combination with the same methods previously described, but with honest causal
forests instead of traditional “adaptive” causal forests. The resulting
average MSEs are presented in the Appendix; we found similar results to the
original 1,000 repetitions with adaptive causal forests, but the honest causal
forests had slightly higher MSEs on average, indicating worse estimation
accuracy than the adaptive causal forests. The honest causal forests had
higher average MSE compared to the X-learner with each of the aggregation
methods except for pooling with trial indicator (Figure A6).

##  6 Application to Real Dataset

After the simulations demonstrated differences across methods in several data
generation setups, we applied the various methods to an existing dataset
containing multiple randomized controlled trials that compared the same two
medications.

###  6.1 Treatments for Major Depressive Disorder

The applied dataset used in the current paper consists of four randomized
controlled trials (Mahableshwarkar et al., 2013, 2015; Boulenger et al., 2014;
Baldwin et al., 2012), each of which included three treatments: duloxetine,
vortioxetine, and placebo, where duloxetine and vortioxetine are both
treatments for major depressive disorder (MDD). At the time of the trials,
duloxetine had been more commonly used to treat MDD so was primarily included
in the trials as an active reference, while vortioxetine was a newer treatment
not yet marketed Schatzberg et al. (2014). Each of the four trials compared at
least two different dosages of vortioxetine and therefore had more
participants taking vortioxetine as opposed to duloxetine or placebo. For the
purposes of the current application, we removed placebo participants and
lumped all dosages of vortioxetine together to investigate the potential
differences between the efficacy of the active medications (duloxetine and
vortioxetine), as well as identify features that might be moderating this
difference.

Participants in each of the four trials shared similar eligibility criteria.
All four trials required patients to be between the ages of 18 to 75, to have
a Major Depressive Episode (MDE) as a primary diagnosis according to the DSM-
IV-TR criteria over at least three months, and to have a Montgomery-Asberg
Depression Rating Scale (MADRS) Montgomery and Åsberg (1979) score of at least
22 (one trial) or 26 (three trials) at both screening and
baseline(Mahableshwarkar et al., 2013, 2015; Boulenger et al., 2014; Baldwin
et al., 2012). A primary outcome in the trials is the change in MADRS
(Montgomery-Asberg Depression Rating Scale) Montgomery and Åsberg (1979) score
from baseline to the last observed follow-up in the study. Participants were
meant to stay in the study for 8 weeks, at which point their final MADRS score
was collected. For those who did not remain in the trial for 8 weeks, a last
observation carried forward imputation approach was used for simplicity. This
imputation approach is not the best way to account for missing data and many
other options exist (Little et al., 2012), but it is used here for simplicity
because this example is primarily illustrative. Predictors/effect modifiers
used in the models were age, sex (female or male), smoking status (ever smoked
or never smoked), weight, baseline MADRS score, baseline HAM-A (Hamilton
Anxiety Rating) score Hamilton (1959), comorbidity indicators (if ever had
diabetes mellitus, hypothyroidism, anxiety), and medication indicators (if
they are concomittantly taking an antidepressant, antipsychotic, thyroid
medication). Since the outcome is the difference in MADRS score (MADRS at
follow up minus MADRS at baseline), a more negative outcome indicates a better
result.

After removing any individuals with missing treatment assignment or missing
outcomes and removing individuals who were assigned to the placebo group,
sample sizes were 575, 436, 418, and 418 for each of the trials. Further
descriptive information about the samples in the four RCTs is reported in the
Appendix (Table A3). Little missing covariate data was present in the sample;
however, conditional mean imputation was performed for missing values of
weight (n=1) and baseline HAM-A score (n=2).

Following data preparation, we used each of the aforementioned method
combinations (i.e., causal forest, S-learner, and X-learner as single-study
methods paired with complete pooling, pooling with trial indicator, ensemble
tree, ensemble forest, and ensemble lasso) to estimate the CATEs for every
individual across the four trials. We then compared the CATE estimates across
methods to see their concordance levels. Notably, it is not possible to
compare the method performances with the truth, as the true CATEs are unknown
in this real dataset.

###  6.2 Results

All methods broadly led to the conclusion of a positive average CATE. This
indicates that vortioxetine is estimated to have less of a beneficial effect
on the MADRS score on average. In each of the four RCTs, both treatments were
associated with a reduction in depressive symptom severity over time (shown
through a reduction in MADRS score), but this reduction was smaller for the
vortioxetine group than the duloxetine group. Table 1 contains the mean and
standard deviation of the CATEs according to each method. Broadly, the
S-learner approaches estimated lower CATEs on average than the other
approaches, and there is some consistency between the aggregation approaches
within each single-study method (S-learner, X-learner, and causal forest).
There were especially high levels of similarity in the average CATE estimates
across the causal forest methods, shown in the last column of Table 1. The
variability of the CATE estimates differs depending on the approach as well;
causal forest approaches had higher standard deviations than approaches that
used the S-learner and X-learner. Complete pooling also yielded the highest
standard deviations for CATE estimates out of all of the aggregation
approaches. As a comparison point, we used a multiple linear regression to
estimate an average treatment effect of 2.39 (SE = 0.5), which is similar to
the averages of the CATEs according to the X-learner and causal forest
approaches.

| S-Learner | X-Learner | Causal Forest
---|---|---|---
Complete Pooling | 1.38 (1.6) | 2.57 (1.4) | 2.37 (2.8)
Pooling with Trial Indicator | 0.91 (1.3) | 2.52 (1.3) | 2.37 (2.7)
Ensemble Tree | 0.89 (1.3) | 2.35 (1.5) | 2.23 (2.5)
Ensemble Forest | 0.89 (1.1) | 2.36 (1.4) | 2.30 (2.2)
Ensemble Lasso | 0.89 (1.2) | 2.32 (1.4) | 2.23 (2.1)
Table 1: Mean (SD) of CATEs from all individuals in sample according to
different single-study and aggregation method combinations.*
*The CATEs are individual-level estimates that indicate the difference in the estimated effect of vortioxetine versus duloxetine on the difference in MADRS score for a given patient. A positive CATE indicates that vortioxetine is estimated to have a smaller reduction of the MADRS score.

We then focused in on the resulting model from the causal forest with pooling
with trial indicator, since that approach performed the best on average in the
simulations. The CATE estimates according to this approach with their
associated 95% confidence intervals are displayed in Figure 3. These estimates
support that the majority of individuals have a positive CATE estimate, but
they also display very high levels of uncertainty, with all confidence
intervals including zero.

Figure 3: Point estimates and 95% confidence intervals for CATEs according to
causal forest with pooling with trial indicator.

To learn more about the moderation within the CATE model, we can explore
variable importance measures. In the causal forest as applied through the grf
package Athey et al. (2019), variable importance is calculated as a weighted
sum of the number of times the variable was used in a split at each level of
the forest. Figure 4 displays the variable importance measures according to
the grf package (Athey et al., 2019) for all covariates, first in separate
causal forests for each study (4(a)), and second according to the causal
forest with pooling with trial indicator (4(b)). From Figure 4(a), there are a
few variables that are consistently identified as effect moderators across
studies (age, weight, baseline MADRS score, and baseline HAM-A score), and
there are several that are not found to be major moderators (the comorbidity
and medication indicators). However, notably there are some differences
according to the separate models, indicating that the treatment effect
functions are slightly different within each study. Figure 4(b) then displays
the resulting importance measures from one aggregation model fit to all
studies. Here, we can see that the same four variables (age, weight, baseline
MADRS, and baseline HAM-A) are involved in a high proportion of the splits in
the causal forest, and study membership is involved in some splits as well.
The fact that these study indicators are not more highly involved in the
partitioning of the treatment effect is a good sign, though, that there is not
a very high level of heterogeneity in CATE estimates across studies.

(a) Variable importance for study-specific causal forest models.

(b) Variable importance for causal forest with pooling with trial indicator

Figure 4: Variable importance measures (a) within studies, and (b) according
to the causal forest with pooling with trial indicator

The variable importance plots do not demonstrate the direction of the
moderating effect, however. We briefly investigate these directional effects
through an interpretation tree (Figure 5) and through exploratory plots such
as Figure A7. This interpretation tree was formed by fitting a regression
tree, where the CATE estimates according to the causal forest with pooling
with trial indicator were the outcomes, and the features (predictors) were
every covariate in the original CATE model. The tree confirms what was shown
in Figure 4 – that age, weight, baseline MADRS, and baseline HAM-A score are
the strongest effect moderators. Study membership does not show up in this
interpretation tree, supporting that there is low heterogeneity across trials.
This is a helpful visual to see the direction of the relationships aggregated
across trials, but it is exploratory and should not be interpreted in great
detail. Another similar approach for investigating the CATE function in terms
of individual moderators is to fit the best linear projection of the CATE
estimates using a function in the grf package Athey et al. (2019); the
resulting coefficients from this regression using doubly-robust estimates of
the CATE are reported in Table A4.

Figure 5: Interpretation tree for causal forest with pooling with trial
indicator.*
*Circled numbers represent the average CATE estimate for individuals in that leaf.

Broadly, these interpretations of the CATE function derived by the causal
forest with pooling with trial indicator do not display high levels of
heterogeneity or highlight clear groups of individuals who might differ in
their treatment effect. Although the interpretation tree shows a partition of
the covariate space into groups that have different average CATEs, the other
plots and best linear projection do not show high levels of heterogeneity,
except for slightly higher average CATE estimates for older individuals. This
data application shows how to effectively apply the methods compared in
simulations to a real dataset and assess potential moderation. The methods all
agree broadly on the direction of the average treatment effect but have some
differences in the individual CATE estimates and which variables generate
effect heterogeneity.

##  7 Discussion

In this paper, we compared methods to estimate the conditional average
treatment effect in a single trial and methods to extend the single-trial
approaches to multiple trials. In the absence of notable cross-trial
heterogeneity of treatment effects, the methods examined all performed well,
but when trial membership was involved in the treatment effect function, some
methods performed worse than others. Specifically, and not surprisingly,
methods that ignore trial membership (complete pooling) do not effectively
estimate the CATE when there is cross-trial heterogeneity. On the other hand,
some methods performed well no matter the level of heterogeneity: pooling with
trial indicator and ensemble forests had consistently low mean squared error
despite increasing the variability of the trial membership coefficients in the
treatment effect. This was especially true when the single-study method used
was the causal forest (Figure 1).

When considering the three single-study approaches, the most consistently
favorable method in the simulations was the causal forest, followed by the
X-learner. The S-learner performed well in certain scenarios, such as scenario
1b, where the treatment effect function involved a bounded, non-linear expit
function. The performances of the S-learner and X-learner in our simulations
and applied example were consistent with results found previously (Künzel et
al., 2019), in that the S-learner seemed to be somewhat biased towards 0 in
the applied example (Table 1) and performed worse in the simulations when the
treatment effect function was complicated (variable CATE scenario and the
piecewise linear and non-linear CATE with high variability). The X-learner
performed well in the simulations with complex CATEs and with structural forms
of the CATE, again consistent with previous work (Künzel et al., 2019). The
causal forest performed well across all scenarios. These simulation results
and the results from the applied data example of MDD medications demonstrate
that it is important to carefully select the single-study method for a given
question, as each of the three options can provide different estimates. A good
starting point would be to consider expert knowledge of how heterogeneous
across studies and complicated the outcomes or treatment effect might be.
These results also indicate the need for more diagnostics in future work to
help researchers determine which approach to choose.

The simulations also incorporated some comparisons between the non-parametric
and parametric approaches. Specifically, the usage of a lasso regression as an
ensemble showed how a parametric ensemble could perform compared to the
ensemble tree and forest. The lasso performed very well when the treatment
effect function was piecewise linear (scenario 1a) but quickly suffered in
performance when the function was more non-linear (scenarios 1b and 2).
Furthermore, the inclusion of a mixed effects meta-analysis demonstrated a
common parametric technique used in the multiple-study setting. This model was
set up to perform well when the CATE function was piecewise linear (scenario
1a), but it yielded high MSE in the non-linear and complex scenarios that it
was not correctly parametrized for (scenario 1b and 2). These comparisons
demonstrate that non-parametric machine learning approaches are very
beneficial when the treatment effect function is complicated and non-linear,
as the non-parametric methods do not require correct specification of any
parameters. Although interpretability becomes more of a challenge, the non-
parametric methods allow for flexible relationships and hopefully high levels
of accuracy in CATE estimation.

In this work, we did not explore an exhaustive list of potential single-study
and aggregation methods. We attempted to select single-study methods that were
common, user-friendly, and shown to be effective or potentially effective in
previous literature. However, as this is an ever-growing field, future work
could include other single-study methods (Wendling et al., 2018; Powers et
al., 2018) to see how they compare to the ones used in this study. For
example, it would be interesting to investigate the performance of the
X-learner with a different base learner, such as Bayesian additive regression
trees (BART) Chipman et al. (2010); Künzel et al. (2019). The simulations
included three scenarios and some sub-scenarios that should be possible in
real data, but there of course are many other possible data generation setups
that could be implemented and could show slightly different results.

It is important to note that with the exception of the complete pooling and
meta-analysis approaches, the resulting CATE estimates are trial-specific.
Unless trial was not picked up in the aggregation methods, the majority of the
methods discussed will produce trial-specific estimates of the CATE. This
allows for improved accuracy of estimates but might be less helpful in real
world applications. We are interested in continuing to identify ways in which
researchers could aggregate across trials to develop estimates that are
accurate but not trial-specific – this could be crucial for use of the
resulting methods and models in practice, on data not coming from the specific
trials used in the model formulation. However, the trial-specific estimates
can still be useful; for example, if trials were done in separate hospitals,
CATEs of future patients could be predicted using the hospital that they are
being treated in, and the model that estimates their treatment effect should
be more accurate after taking into consideration the data from the other
hospitals.

In the MDD trials, duloxetine was included as a reference medication because
it was already marketed at the time of the trials, and patients were excluded
from the study if they had previously not responded to duloxetine. On the
other hand, vortioxetine was not yet marketed and was the more experimental
medication; therefore, some bias could arise due to participants being
excluded if they had previously not responded to duloxetine. Acknowledging
this, we were able to estimate treatment effects according to each method
combination, and we used variable importance and interpretation trees to
investigate which variables might be important moderators of the treatment
effect. Variable importance is a limited measure and can often be biased
towards continuous variables with more possible split points Strobl et al.
(2007), so we encourage caution when interpreting those results. This example
dataset shows how to combine multiple RCTs to get an improved assessment of
treatment effect heterogeneity and better determine which treatment would be
best suited to a given individual, based on their features and their estimated
CATE. Notably, the four trials used in this dataset were run by the same
organizations and had very similar protocols; this helps ensure that we can
confidently combine datasets but also might limit the potential heterogeneity
across trials that might exist in other applications. We also did not see high
levels of heterogeneity in the treatment effects based on other covariates in
these trials. A general idea is that studies need to be four times larger to
identify effect moderators compared to an average treatment effect Enderlein
(1988), and this study included precisely four trials. Therefore, our findings
would become more robust and we could more confidently assess heterogeneity or
lack thereof with the inclusion of more studies.

There are many openings for future work, some of which have been previously
mentioned. Broadly, it is important to further refine these methods and
identify which are most helpful in specific data scenarios. It will also be
helpful to determine when it is appropriate to develop universal CATE
estimates, versus when the CATE estimates should be trial-specific. This paper
demonstrated several approaches that take data from multiple studies and
estimate heterogeneous treatment effects, using flexible models that allow for
complex relationships - which is often the case in the real world.

##  8 Acknowledgments

The study was funded by the Patient-Centered Outcomes Research Institute
(PCORI) through PCORI Award ME-2020C3-21145 (PI: Stuart) and the National
Institute of Mental Health (NIMH) through Award R01MH126856 (PI: Stuart).
Disclaimer: Opinions and information in this content are those of the study
authors and do not necessarily represent the views of PCORI or NIMH.
Accordingly, PCORI and NIMH cannot make any guarantees with respect to the
accuracy or reliability of the information and data.

Furthermore, this paper is based on research using data from data
contributors, Takeda and Lundbeck, that has been made available through Vivli,
Inc. Vivli has not contributed to or approved, and is not in any way
responsible for, the contents of this publication.

## References

  * Enderlein [1988] G. Enderlein.  Fleiss, J. L.: The Design and Analysis of Clinical Experiments. Wiley, New York – Chichester – Brislane – Toronto – Singapore 1986, 432 S., £38.35.  _Biometrical Journal_ , 30(3):304–304, 1988.  ISSN 1521-4036.  doi:10.1002/bimj.4710300308.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.4710300308.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/bimj.4710300308.
  * Debray et al. [2015] Thomas P. A. Debray, Karel G. M. Moons, Gert Valkenhoef, Orestis Efthimiou, Noemi Hummel, Rolf H. H. Groenwold, Johannes B. Reitsma, and on behalf of the GetReal methods review group.  Get real in individual participant data (IPD) meta-analysis: a review of the methodology.  _Research Synthesis Methods_ , 6(4):293–309, December 2015.  ISSN 1759-2879, 1759-2887.  doi:10.1002/jrsm.1160.  URL https://onlinelibrary.wiley.com/doi/10.1002/jrsm.1160.
  * Seo et al. [2021] Michael Seo, Ian R. White, Toshi A. Furukawa, Hissei Imai, Marco Valgimigli, Matthias Egger, Marcel Zwahlen, and Orestis Efthimiou.  Comparing methods for estimating patient-specific treatment effects in individual patient data meta-analysis.  _Statistics in Medicine_ , 40(6):1553–1573, March 2021.  ISSN 0277-6715, 1097-0258.  doi:10.1002/sim.8859.  URL https://onlinelibrary.wiley.com/doi/10.1002/sim.8859.
  * Künzel et al. [2019] Sören R Künzel, Jasjeet S Sekhon, Peter J Bickel, and Bin Yu.  Metalearners for estimating heterogeneous treatment effects using machine learning.  _Proceedings of the national academy of sciences_ , 116(10):4156–4165, 2019.
  * Athey et al. [2019] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  _The Annals of Statistics_ , 47(2), April 2019.  ISSN 0090-5364.  doi:10.1214/18-AOS1709.  URL https://projecteuclid.org/journals/annals-of-statistics/volume-47/issue-2/Generalized-random-forests/10.1214/18-AOS1709.full.
  * Green and Kern [2012] D. P. Green and H. L. Kern.  Modeling Heterogeneous Treatment Effects in Survey Experiments with Bayesian Additive Regression Trees.  _Public Opinion Quarterly_ , 76(3):491–511, September 2012.  ISSN 0033-362X, 1537-5331.  doi:10.1093/poq/nfs036.  URL https://academic.oup.com/poq/article-lookup/doi/10.1093/poq/nfs036.
  * Kennedy [2020] Edward H. Kennedy.  Optimal doubly robust estimation of heterogeneous causal effects.  _arXiv:2004.14497 [math, stat]_ , June 2020.  URL http://arxiv.org/abs/2004.14497.  arXiv: 2004.14497.
  * Nie and Wager [2021] X Nie and S Wager.  Quasi-oracle estimation of heterogeneous treatment effects.  _Biometrika_ , 108(2):299–319, May 2021.  ISSN 0006-3444, 1464-3510.  doi:10.1093/biomet/asaa076.  URL https://academic.oup.com/biomet/article/108/2/299/5911092.
  * Dandl et al. [2022] Susanne Dandl, Torsten Hothorn, Heidi Seibold, Erik Sverdrup, Stefan Wager, and Achim Zeileis.  What makes forest-based heterogeneous treatment effect estimators work?  _arXiv preprint arXiv:2206.10323_ , 2022.
  * Yang et al. [2020] Shu Yang, Donglin Zeng, and Xiaofei Wang.  Elastic Integrative Analysis of Randomized Trial and Real-World Data for Treatment Heterogeneity Estimation.  _arXiv:2005.10579 [stat]_ , September 2020.  URL http://arxiv.org/abs/2005.10579.  arXiv: 2005.10579.
  * Yang et al. [2022] Shu Yang, Donglin Zeng, and Xiaofei Wang.  Improved Inference for Heterogeneous Treatment Effects Using Real-World Data Subject to Hidden Confounding.  _arXiv:2007.12922 [stat]_ , January 2022.  URL http://arxiv.org/abs/2007.12922.  arXiv: 2007.12922.
  * Kallus et al. [2018] Nathan Kallus, Aahlad Manas Puli, and Uri Shalit.  Removing Hidden Confounding by Experimental Grounding.  _arXiv:1810.11646 [cs, stat]_ , October 2018.  URL http://arxiv.org/abs/1810.11646.  arXiv: 1810.11646.
  * Rosenman et al. [2020] Evan Rosenman, Guillaume Basse, Art Owen, and Michael Baiocchi.  Combining Observational and Experimental Datasets Using Shrinkage Estimators.  _arXiv:2002.06708 [math, stat]_ , May 2020.  URL http://arxiv.org/abs/2002.06708.  arXiv: 2002.06708.
  * Brantner et al. [2023] Carly Lupton Brantner, Ting-Hsuan Chang, Trang Quynh Nguyen, Hwanhee Hong, Leon Di Stefano, and Elizabeth A Stuart.  Methods for integrating trials and non-experimental data to examine treatment effect heterogeneity.  _arXiv preprint arXiv:2302.13428_ , 2023.
  * Tan et al. [2021] Xiaoqing Tan, Chung-Chou H. Chang, and Lu Tang.  A Tree-based Federated Learning Approach for Personalized Treatment Effect Estimation from Heterogeneous Data Sources.  _arXiv:2103.06261 [cs, stat]_ , April 2021.  URL http://arxiv.org/abs/2103.06261.  arXiv: 2103.06261.
  * Rubin [1974] Donald B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688–701, 1974.  ISSN 1939-2176.  doi:10.1037/h0037350.  Place: US Publisher: American Psychological Association.
  * Wager and Athey [2018] Stefan Wager and Susan Athey.  Estimation and Inference of Heterogeneous Treatment Effects using Random Forests.  _Journal of the American Statistical Association_ , 113(523):1228–1242, July 2018.  ISSN 0162-1459, 1537-274X.  doi:10.1080/01621459.2017.1319839.  URL https://www.tandfonline.com/doi/full/10.1080/01621459.2017.1319839.
  * Therneau et al. [2015] Terry Therneau, Beth Atkinson, Brian Ripley, and Maintainer Brian Ripley.  Package ‘rpart’.  _Available online: cran. ma. ic. ac. uk/web/packages/rpart/rpart. pdf (accessed on 20 April 2016)_ , 2015.
  * Wright and Ziegler [2015] Marvin N Wright and Andreas Ziegler.  ranger: A fast implementation of random forests for high dimensional data in c++ and r.  _arXiv preprint arXiv:1508.04409_ , 2015.
  * Friedman et al. [2017] Jerome Friedman, Trevor Hastie, Noah Simon, Rob Tibshirani, Maintainer Trevor Hastie, and Depends Matrix.  Package ‘glmnet.’.  _Journal of statistical software_ , 33(1):1–22, 2017.
  * Bates [2010] Douglas M Bates.  lme4: Mixed-effects modeling with r, 2010.
  * Mahableshwarkar et al. [2013] Atul R. Mahableshwarkar, Paula L. Jacobsen, and Yinzhong Chen.  A randomized, double-blind trial of 2.5 mg and 5 mg vortioxetine (Lu AA21004) versus placebo for 8 weeks in adults with major depressive disorder.  _Current Medical Research and Opinion_ , 29(3):217–226, March 2013.  ISSN 0300-7995, 1473-4877.  doi:10.1185/03007995.2012.761600.  URL http://www.tandfonline.com/doi/full/10.1185/03007995.2012.761600.
  * Mahableshwarkar et al. [2015] Atul R. Mahableshwarkar, Paula L. Jacobsen, Yinzhong Chen, Michael Serenko, and Madhukar H. Trivedi.  A randomized, double-blind, duloxetine-referenced study comparing efficacy and tolerability of 2 fixed doses of vortioxetine in the acute treatment of adults with MDD.  _Psychopharmacology_ , 232(12):2061–2070, June 2015.  ISSN 0033-3158, 1432-2072.  doi:10.1007/s00213-014-3839-0.  URL http://link.springer.com/10.1007/s00213-014-3839-0.
  * Boulenger et al. [2014] Jean-Philippe Boulenger, Henrik Loft, and Christina Kurre Olsen.  Efficacy and safety of vortioxetine (Lu AA21004), 15 and 20 mg/day: a randomized, double-blind, placebo-controlled, duloxetine-referenced study in the acute treatment of adult patients with major depressive disorder.  _International Clinical Psychopharmacology_ , 29(3):138–149, May 2014.  ISSN 0268-1315.  doi:10.1097/YIC.0000000000000018.  URL http://journals.lww.com/00004850-201405000-00002.
  * Baldwin et al. [2012] David S Baldwin, Henrik Loft, and Marianne Dragheim.  A randomised, double-blind, placebo controlled, duloxetine-referenced, fixed-dose study of three dosages of lu aa21004 in acute treatment of major depressive disorder (mdd).  _European Neuropsychopharmacology_ , 22(7):482–491, 2012.
  * Schatzberg et al. [2014] Alan F Schatzberg, Pierre Blier, Larry Culpepper, Rakesh Jain, George I Papakostas, and Michael E Thase.  An overview of vortioxetine.  _The Journal of Clinical Psychiatry_ , 75(12):13677, 2014.
  * Montgomery and Åsberg [1979] Stuart A Montgomery and MARIE Åsberg.  A new depression scale designed to be sensitive to change.  _The British journal of psychiatry_ , 134(4):382–389, 1979.
  * Little et al. [2012] Roderick J Little, Ralph D’Agostino, Michael L Cohen, Kay Dickersin, Scott S Emerson, John T Farrar, Constantine Frangakis, Joseph W Hogan, Geert Molenberghs, Susan A Murphy, et al.  The prevention and treatment of missing data in clinical trials.  _New England Journal of Medicine_ , 367(14):1355–1360, 2012.
  * Hamilton [1959] MAX Hamilton.  The assessment of anxiety states by rating.  _British journal of medical psychology_ , 1959.
  * Wendling et al. [2018] T. Wendling, K. Jung, A. Callahan, A. Schuler, N. H. Shah, and B. Gallego.  Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases.  _Statistics in Medicine_ , 37(23):3309–3324, October 2018.  ISSN 02776715.  doi:10.1002/sim.7820.  URL https://onlinelibrary.wiley.com/doi/10.1002/sim.7820.
  * Powers et al. [2018] Scott Powers, Junyang Qian, Kenneth Jung, Alejandro Schuler, Nigam H. Shah, Trevor Hastie, and Robert Tibshirani.  Some methods for heterogeneous treatment effect estimation in high dimensions: Some methods for heterogeneous treatment effect estimation in high dimensions.  _Statistics in Medicine_ , 37(11):1767–1787, May 2018.  ISSN 02776715.  doi:10.1002/sim.7623.  URL https://onlinelibrary.wiley.com/doi/10.1002/sim.7623.
  * Chipman et al. [2010] Hugh A Chipman, Edward I George, and Robert E McCulloch.  Bart: Bayesian additive regression trees.  _The Annals of Applied Statistics_ , 4(1):266–298, 2010.
  * Strobl et al. [2007] Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis, and Torsten Hothorn.  Bias in random forest variable importance measures: Illustrations, sources and a solution.  _BMC bioinformatics_ , 8(1):1–21, 2007.

| Piecewise Linear (1a) | Non-Linear (1b) | Variable (2)
---|---|---|---
| Low-Low | Med-Low | Med-Med | Med-High | High-High | Low-Low | Med-Low | Med-Med | Med-High | High-High |
S - Pool | 0.43 (0.11) | 1.07 (0.42) | 1.3 (0.43) | 1.94 (0.58) | 8.75 (3.56) | 0.2 (0.09) | 0.77 (0.36) | 0.98 (0.4) | 1.57 (0.53) | 7.9 (3.33) | 0.42 (0.02)
X - Pool | 0.25 (0.11) | 0.92 (0.45) | 1.16 (0.45) | 1.82 (0.6) | 9.06 (3.76) | 0.22 (0.1) | 0.86 (0.4) | 1.1 (0.45) | 1.75 (0.58) | 8.83 (3.72) | 0.29 (0.01)
CF - Pool | 0.26 (0.12) | 0.99 (0.48) | 1.25 (0.48) | 1.98 (0.66) | 9.8 (4.04) | 0.23 (0.1) | 0.89 (0.42) | 1.15 (0.46) | 1.84 (0.61) | 9.17 (3.82) | 0.31 (0.01)
S - Indicator | 0.41 (0.09) | 0.6 (0.16) | 0.77 (0.18) | 1.27 (0.39) | 0.91 (0.6) | 0.02 (0.02) | 0.05 (0.07) | 0.08 (0.08) | 0.15 (0.14) | 0.49 (0.7) | 0.43 (0.02)
X - Indicator | 0.09 (0.05) | 0.22 (0.16) | 0.3 (0.18) | 0.53 (0.29) | 1.83 (1.58) | 0.06 (0.04) | 0.18 (0.16) | 0.25 (0.17) | 0.44 (0.27) | 1.78 (1.45) | 0.18 (0.01)
CF - Indicator | 0.03 (0) | 0.03 (0.01) | 0.05 (0.01) | 0.09 (0.03) | 0.11 (0.03) | 0.02 (0) | 0.02 (0) | 0.04 (0.01) | 0.07 (0.02) | 0.07 (0.02) | 0.03 (0.01)
S - Tree | 0.5 (0.05) | 0.52 (0.06) | 0.57 (0.1) | 0.74 (0.16) | 0.95 (0.28) | 0.21 (0.05) | 0.19 (0.04) | 0.34 (0.08) | 0.6 (0.18) | 1.25 (0.41) | 0.57 (0.03)
X - Tree | 0.11 (0.02) | 0.15 (0.03) | 0.21 (0.05) | 0.34 (0.11) | 0.78 (0.28) | 0.16 (0.02) | 0.18 (0.02) | 0.31 (0.07) | 0.53 (0.14) | 1.29 (0.4) | 0.16 (0.02)
CF - Tree | 0.06 (0.01) | 0.1 (0.03) | 0.17 (0.05) | 0.3 (0.1) | 0.74 (0.27) | 0.13 (0.02) | 0.15 (0.02) | 0.28 (0.06) | 0.49 (0.13) | 1.18 (0.4) | 0.06 (0.01)
S - Forest | 0.48 (0.04) | 0.46 (0.05) | 0.48 (0.08) | 0.58 (0.1) | 0.41 (0.09) | 0.13 (0.05) | 0.09 (0.04) | 0.12 (0.05) | 0.23 (0.12) | 0.14 (0.07) | 0.55 (0.03)
X - Forest | 0.08 (0.01) | 0.08 (0.01) | 0.1 (0.01) | 0.15 (0.04) | 0.15 (0.04) | 0.07 (0.01) | 0.07 (0.01) | 0.09 (0.01) | 0.16 (0.04) | 0.16 (0.04) | 0.17 (0.02)
CF - Forest | 0.04 (0) | 0.04 (0) | 0.06 (0.01) | 0.1 (0.03) | 0.11 (0.03) | 0.04 (0) | 0.04 (0) | 0.06 (0.01) | 0.11 (0.04) | 0.11 (0.03) | 0.07 (0.01)
S - Lasso | 0.5 (0.05) | 0.48 (0.06) | 0.5 (0.09) | 0.56 (0.1) | 0.34 (0.08) | 0.78 (0.03) | 0.77 (0.03) | 0.78 (0.03) | 0.84 (0.08) | 0.8 (0.05) | 0.64 (0.03)
X - Lasso | 0.13 (0.01) | 0.13 (0.01) | 0.13 (0.01) | 0.14 (0.01) | 0.15 (0.02) | 0.76 (0.02) | 0.76 (0.02) | 0.77 (0.02) | 0.79 (0.03) | 0.79 (0.02) | 0.38 (0.02)
CF - Lasso | 0.1 (0) | 0.1 (0) | 0.11 (0.01) | 0.12 (0.01) | 0.12 (0.01) | 0.75 (0.02) | 0.75 (0.02) | 0.76 (0.02) | 0.77 (0.02) | 0.77 (0.02) | 0.35 (0.01)
Meta-Analysis | 0.09 (0.03) | 0.1 (0.09) | 0.1 (0.1) | 0.1 (0.11) | 0.22 (1.11) | 0.93 (0.03) | 0.94 (0.11) | 0.93 (0.09) | 0.94 (0.13) | 1.01 (0.91) | 0.41 (0.02)

Table 2: Average MSE for each method and data generation setup.*
*The lowest average MSE for each scenario (column) is bolded. Figure 6: Average MSE across all scenarios and iterations using honest causal forests. Table 3: Descriptive statistics of participants of four randomized controlled trials, broken down by treatment group. Figure 7: CATE estimates by age of individual according to causal forest with pooling with trial indicator.*
*Note that uncertainty of the CATE estimates is not reflected in this plot. | Estimate | Standard Error | P-Value
---|---|---|---
(Intercept) | -6.32 | 5.06 | 0.21
Age | 0.09 | 0.04 | 0.03*
Female | 0.45 | 1.07 | 0.67
Smoker | -1.47 | 1.10 | 0.18
Weight | -0.01 | 0.03 | 0.72
Baseline MADRS | 0.09 | 0.13 | 0.49
Baseline HAM-A | 0.08 | 0.09 | 0.38
Has Diabetes Mellitus | -3.97 | 3.36 | 0.24
Has Hypothyroidism | -1.81 | 3.56 | 0.61
Has Anxiety | 3.58 | 3.63 | 0.32
Takes Antidepressant | 1.55 | 1.32 | 0.24
Takes Antipsychotic | -0.21 | 1.93 | 0.91
Takes Thyroid Medication | 2.23 | 4.22 | 0.60
Study NCT00635219 | -1.09 | 1.63 | 0.50
Study NCT01140906 | 2.71 | 1.62 | 0.09
Study NCT00672620 | 2.93 | 1.58 | 0.06
Table 4: Results of best linear projection of the CATE according to the causal
forest with pooling with trial indicator.
*Indicates a p-value less than 0.05.