# Heterogeneous Treatment Effect Estimation for Observational Data using
Model-based Forests

Susanne Dandl
LMU München  0000-0003-4324-4163 MCML
&Andeas Bender
LMU München  0000-0001-5628-8611 MCML
&Torsten Hothorn
Universität Zürich  0000-0001-8301-0471<EMAIL_ADDRESS>
Dandl, Bender, Hothorn Forest-based HTE Estimation The estimation of
heterogeneous treatment effects (HTEs) has attracted considerable interest in
many disciplines, most prominently in medicine and economics. Contemporary
research has so far primarily focused on continuous and binary responses where
HTEs are traditionally estimated by a linear model, which allows the
estimation of constant or heterogeneous effects even under certain model
misspecifications. More complex models for survival, count, or ordinal
outcomes require stricter assumptions to reliably estimate the treatment
effect. Most importantly, the noncollapsibility issue necessitates the joint
estimation of treatment and prognostic effects. Model-based forests allow
simultaneous estimation of covariate-dependent treatment and prognostic
effects, but only for randomized trials. In this paper, we propose
modifications to model-based forests to address the confounding issue in
observational data. In particular, we evaluate an orthogonalization strategy
originally proposed by Robinson (1988, Econometrica) in the context of model-
based forests targeting HTE estimation in generalized linear models and
transformation models. We found that this strategy reduces confounding effects
in a simulated study with various outcome distributions. We demonstrate the
practical aspects of HTE estimation for survival and ordinal outcomes by an
assessment of the potentially heterogeneous effect of Riluzole on the progress
of Amyotrophic Lateral Sclerosis. Heterogeneous treatment effects,
personalized medicine, random forest, observational data, censored survival
data, generalized linear model, transformation model Susanne Dandl, Andreas
Bender
Institut für Statistik, LMU München, Germany
Munich Center for Machine Learning (MCML), Germany
Torsten Hothorn
Institut für Epidemiologie, Biostatistik und Prävention, Universität Zürich
Hirschengraben 84, CH-8001 Zürich, Switzerland
E-mail:

##  1 Introduction

Over the past years, there has been emerging interest in methods to estimate
heterogeneous treatment effects (HTEs) in various application fields. In
healthcare, HTE estimation can be understood as a core principle driving
personalized medicine. As opposed to average treatment effects, which assume a
constant effect of a treatment on an outcome for the whole population, HTEs
account for the heterogeneity in the effect for subgroups or individuals based
on their characteristics. Most research on HTE estimation has mainly focused
on continuous and binary response variables. These methods have typically
built upon Rubin’s potential outcomes framework, a statistical approach to
formulating and inferring causal effects in various designs (Rubin, 1974,
2005).

Traditionally, statistical models were used to estimate the treatment effect,
but machine learning methods have been more and more adapted for these tasks
over the past decade. Machine learning models rely on weaker assumptions and
can automatically learn complex relationships such as higher order interaction
effects, resulting in greater predictive performance in a variety of
applications. In the case of continuous or binary responses, prominent methods
to estimate HTEs are based on random forests (Foster _et al._ , 2011; Lu _et
al._ , 2018; Athey _et al._ , 2019; Powers _et al._ , 2018; Su _et al._ ,
2018; Li _et al._ , 2022), Bayesian additive regression trees (BART) (Hill,
2011; Hu _et al._ , 2020), or neural networks (Shalit _et al._ , 2017; Curth
_et al._ , 2021; Chapfuwa _et al._ , 2021). Künzel _et al._ (2019) proposed
general frameworks – T-learners, S-learners, U-learners, and X-learners – that
base treatment effect estimates on arbitrary machine learning models.
Chernozhukov _et al._ (2018) coined the term double/debiased machine learning
models, which uses machine learning models for nuisance parameter estimations.
The approach still relies on parametric models for estimating treatment
effects, but Nie and Wager (2021) derived so-called R-learners that allow for
arbitrary (nonparametric or semiparametric) models.

Beyond continuous or binary responses, research on machine learning methods
for HTE estimation have primarily focused on (right-censored) survival data.
Methods have been proposed based on Bayesian additive regression trees (BART)
(Henderson _et al._ , 2018), random forest-type methods (Cui _et al._ , 2022;
Tabib and Larocque, 2020), or deep learning approaches (Curth _et al._ , 2021;
Chapfuwa _et al._ , 2021). Theoretically, any machine learning model for
survival analysis – such as random survival forests (Ishwaran _et al._ , 2008)
or a Cox regression-based deep neural network (deepSurv) (Katzman _et al._ ,
2018) – can estimate HTEs (Hu _et al._ , 2021). These models can estimate
survival or hazard functions in both treatment groups separately; HTEs are
then defined as the difference in derived properties of the two functions,
e.g., as differences in the median survival time. However, Hu _et al._ (2021)
found that methods specifically designed for HTE estimation, like the adapted
BART (Henderson _et al._ , 2018), produce more reliable estimates.

In general, for a continuous or binary outcome $Y$ conditional on treatment
$w$ and covariates $x$, the conditional average treatment effect
$\tau(\text{\boldmath$x$})$ (CATE) can be estimated from the model
$\mathbb{E}(Y\mid
W=w,\text{\boldmath$X$}=\text{\boldmath$x$})=\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w$
even if the model is misspecified, e.g., when the prognostic effect
$\mu(\text{\boldmath$x$})$ cannot be fully estimated due to missing covariate
information. Beyond mean regression, stricter assumptions are necessary both
for randomized and for observational studies to estimate HTEs. For example,
under a true Cox model with survivor function
$\exp(-\exp(h(t)+\mu(\text{\boldmath$x$})+\tau w))$ with log-cumulative
baseline hazard $h(t)$ at time $t$ and log-hazard ratio $\tau$, the prognostic
effect $\mu(\text{\boldmath$x$})$ must be specified correctly, even in a
randomized trial. Estimated marginal log-hazard ratios $\hat{\tau}$ – i.e.,
when the model is fitted under the constraint $\mu(\text{\boldmath$x$})\equiv
0$ – are shrunken towards zero if this constraint is unrealistic (Aalen _et
al._ , 2015). Naturally, this problem carries over to heterogeneous log-hazard
ratios $\tau(\text{\boldmath$x$})$.

Consequently, HTE estimation in more complex models requires the simultaneous
estimation of both the prognostic part $\mu(\text{\boldmath$x$})$ and the
predictive HTE $\tau(\text{\boldmath$x$})$. Model-based forests have been
demonstrated to allow estimation of $\mu(\text{\boldmath$x$})$ and
$\tau(\text{\boldmath$x$})$ in randomized trials (Seibold _et al._ , 2016,
2018; Korepanova _et al._ , 2020; Buri and Hothorn, 2020; Fokkema _et al._ ,
2018; Hothorn and Zeileis, 2021b). In a nutshell, model-based forests combine
the parametric modeling framework with random forests to estimate individual
treatment effects (Seibold _et al._ , 2018). By using generalized linear
models and transformation models, model-based forests can be adapted for
survival data (Seibold _et al._ , 2016, 2018; Korepanova _et al._ , 2020),
ordinal data (Buri and Hothorn, 2020), or clustered data (Fokkema _et al._ ,
2018). A unique feature of model-based forests is the simultaneous estimation
of both treatment and prognostic effects in the same forest model.

In observational studies the treatment group assignment is not under control
of the researcher and confounding effects could bias the estimation of HTEs.
In this work, we propose and evaluate novel variants of model-based forests
for HTE estimation in observational studies. Adaptions of Robinson’s
orthogonalization strategy for generalized linear models and transformation
models are discussed and implemented. We review key components of model-based
forests for HTE estimation in randomized trials in Section 2. In Section 3, we
start introducing the orthogonalization approach by Robinson (1988), which is
instrumental for achieving robustness to confounding effects in the non-
randomized situation. We motivate previous developments using linear
models(Dandl _et al._ , 2022) and leverage adaptations to more complex models
discussed by Gao and Hastie (2022) to define novel model-based forest variants
suitable for HTE in the observational setting. These variants’ performances
are empirically assessed in a simulation study with a range of outcome
distributions in Section 4. Finally, in Section 5 presenting a re-analysis of
the patient-specific effect of Riluzole in patients with Amyotrophic Lateral
Sclerosis (ALS), practical aspects of model estimation and interpretation are
discussed.

##  2 Review of model-based forests for randomized trials

We are interested in estimating HTEs based on i.i.d. observations
$(y,\text{\boldmath$x$},w)$, where $y$, $x$ and $w$ are realizations of the
outcome $Y$, covariates $\text{\boldmath$X$}\in\mathcal{X}$, and control vs.
treatment indicator $W\in\\{0,1\\}$. $Y(0)$ and $Y(1)$ denote the potential
outcomes under the two treatment conditions $W\in\\{0,1\\}$. Throughout this
paper, we assume that $X$ includes all relevant variables to explain
heterogeneity both in the treatment effect and the outcome $Y$, and that the
base model underlying model-based forests is correctly specified.

We review model-based forests for HTE estimation based on randomized trials as
introduced by Seibold _et al._ (2018) and Korepanova _et al._ (2020). Within
this section, we only consider settings where the treatment assignment is
randomized and, therefore, follows a binomial model
$W\mid\text{\boldmath$X$}=\text{\boldmath$x$}\sim\operatorname{B}(1,\pi(\text{\boldmath$x$}))$
with constant propensities $\pi(\text{\boldmath$x$})\equiv\pi$. We omit
discussion of the abstract framework underlying model-based forests and
instead discuss the important linear, generalized linear (Seibold _et al._ ,
2018), and transformation models (Korepanova _et al._ , 2020) in detail.

###  2.1 Linear model

For a continuous outcome $Y\in\mathbb{R}$ with symmetric error distribution, a
model-based forest might be defined based on the model

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)$
$\displaystyle=$
$\displaystyle\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w+\phi Z$ (1)

where the residuals are given by the error term $\phi Z$ with
$\mathbb{E}(Z|\text{\boldmath$X$},W)=0$ and standard deviation $\phi>0$ (Dandl
_et al._ , 2022). We are mainly interested in estimating
$\tau(\text{\boldmath$x$})$, the treatment effect that depends on predictive
variables in $x$. With model-based forests, however, we also obtain an
estimated value for the prognostic effect $\mu(\text{\boldmath$x$})$, which
depends on prognostic variables in $x$. A variable might be predictive and
prognostic at the same time. We refer to these situations as “overlays”.

Because we assume in this section that $\pi(\text{\boldmath$x$})\equiv\pi$
applies, $W\perp\\!\\!\\!\perp\text{\boldmath$X$}$ holds. Consequently,
$\tau(\text{\boldmath$x$})$ can be interpreted as a CATE

$\tau(\text{\boldmath$x$})=\text{CATE}(\text{\boldmath$x$})=\mathbb{E}(Y(1)-Y(0)\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
(2)

on the absolute scale. To estimate
$(\mu(\text{\boldmath$x$}),\tau(\text{\boldmath$x$}))^{\top}$ the $L_{2}$ loss

$\displaystyle\ell(\mu(\text{\boldmath$x$}),\tau(\text{\boldmath$x$}))=\nicefrac{{1}}{{2}}\left(Y-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w\right)^{2}$
(3)

is minimized w.r.t. $\mu$ and $\tau$ using an ensemble of trees. Inspired by
recursive partitioning techniques (Hothorn _et al._ , 2006; Zeileis _et al._ ,
2008), split variable and split point selection are separated. The split
variable is the variable that has the lowest $p$-value for the bivariate
permutation tests for the $H_{0}$-hypothesis that $\mu$ and $\tau$ are
constant and independent of any split variable. The cut-point is the point of
the chosen split variable at which the score functions

$\displaystyle s(\hat{\mu},\hat{\tau}):=(Y-\hat{\mu}-\hat{\tau}w)(1,w)^{\top}$

in the two resultant subgroups differ the most; details are available in
Appendix 2 of Seibold _et al._ (2018).

Once $B\in\mathbb{N}$ trees were fitted to subsamples of the training data,
predictions for the treatment effect for a new observation $x$ are obtained
via local maximum likelihood aggregation (Hothorn _et al._ , 2004;
Meinshausen, 2006; Lin and Jeon, 2006; Athey _et al._ , 2019; Hothorn and
Zeileis, 2021b). First, for the $i$-th training sample, the frequency
$\alpha_{i}$ with which it falls in the same leaf as $x$ over all $B$ trees is
measured. The obtained weighting vector $(\alpha_{1},...,\alpha_{n})$ is used
as an input for minimizing

$(\hat{\mu}(\text{\boldmath$x$}),\hat{\tau}(\text{\boldmath$x$}))^{\top}=\operatorname*{{arg\,min}}_{\mu,\tau}\sum_{i=1}^{n}\alpha_{i}(\text{\boldmath$x$})\ell_{i}(\mu,\tau)$
(4)

where $\ell_{i}$ denotes the loss for the $i$-th sample. Model-based forests
easily allow adaptions if HTEs for an outcome variable $Y$ that is not well
represented by equation (1) should be estimated. In this case, model-based
forests can build on generalized linear models or transformation models in the
recursive partitioning framework (Zeileis _et al._ , 2008). As detailed in the
following sections, the loss function $\ell$ in equation (3) changes from the
squared error to the negative (partial) log-likelihood of some appropriate
model.

###  2.2 Generalized linear models

When the conditional outcome distribution is better described through a
generalized linear model

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim\text{ExpFam}(\theta(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w),\phi)$

with parameter $\theta$ depending on the additive function
$\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w$, the conditional mean

$\displaystyle
g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w))=\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w=:\eta_{w}(\text{\boldmath$x$})$
(5)

is linear on the scale of a link function $g$. Thus, the interpretation of
$\tau(\text{\boldmath$x$})$ as CATE (2) generally no longer holds. Instead,
the predictive effect is understood as the difference in natural parameters
(DINA (Gao and Hastie, 2022))

$\displaystyle\tau(\text{\boldmath$x$})=\text{DINA}(\text{\boldmath$x$})=\eta_{1}(\text{\boldmath$x$})-\eta_{0}(\text{\boldmath$x$}).$
(6)

In contrast to the linear model case, HTEs $\tau(\text{\boldmath$x$})$ are now
defined on relative scales, such as odds ratios in binary logistic regression
models or multiplicative mean effects in a Poisson or Gaussian model with a
log-link. The negative log-likelihood contribution of some observation
$(Y,\text{\boldmath$x$},w)$ is

$\ell(\mu,\tau,\phi)=-\log(f(Y\mid\theta(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w),\phi))$

with $f$ as the conditional density of an exponential family distribution

$f(Y\mid\theta(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w),\phi).$

Model-based trees and forests (Zeileis _et al._ , 2008; Seibold _et al._ ,
2016, 2018) jointly estimate the prognostic effect $\mu(\text{\boldmath$x$})$
and the predictive effect $\tau(\text{\boldmath$x$})$. The procedure
simultaneously minimizes the negative log-likelihood with respect to
$\mu(\text{\boldmath$x$})$ and $\tau(\text{\boldmath$x$})$. In each node of
the model-based forest, $\mu$, $\tau$, and potentially $\phi$ are estimated by
minimizing

$\displaystyle\ell(\mu,\tau,\phi)=-\log(f(Y\mid\theta(\mu+\tau w),\phi))$ (7)

and regressing the bivariate gradient

$\displaystyle\left.\frac{\partial\ell(\mu,\tau,\phi)}{\partial(\mu,\tau)}\right|_{\hat{\mu},\hat{\tau},\hat{\phi}}$

on $x$. This means that one is not explicitly looking for changes in the scale
parameter $\phi$, but this could be implemented by looking at the three-
variate gradient

$\displaystyle\left.\frac{\partial\ell(\mu,\tau,\phi)}{\partial(\mu,\tau,\phi)}\right|_{\hat{\mu},\hat{\tau},\hat{\phi}}$

for example, in a heteroscedastic normal linear model

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w+\phi(\text{\boldmath$x$})Z.$

After the tree fitting phase, a HTE is estimated with equation (4) with
$\ell(\mu,\tau,\phi)$ of equation (7) as the corresponding loss function.

Thus, model-based forests can be directly applied to estimate HTEs on relative
scales for binary outcomes (binary logistic or probit regression, for
example), counts (Poisson or quasi-Poisson regression), or continuous outcomes
where a multiplicative effect is of interest (normal model with log-link).

###  2.3 Transformation models

More complex responses like ordered categorical or time-to-event outcomes are
not covered by generalized linear models but can be analysed using
transformation models; corresponding model-based forests for survival analysis
have been introduced by Korepanova _et al._ (2020). For some at least ordered
outcome $Y$, we write the conditional distribution function as

$\displaystyle\mathbb{P}(Y\leq
y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=F(h(y)-\underbrace{(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w)}_{=:\eta_{w}(\text{\boldmath$x$})}).$
(8)

The transformation function $h$ is monotone non-decreasing and the inverse
link function $F$ governs the interpretability of $\tau$ as log-odds ratios
($F=\operatorname{logit}^{-1}$), log-hazard ratios
($F=\operatorname{cloglog}^{-1}$), log-reverse time hazard ratios
($F=\operatorname{loglog}^{-1}$), or shift effects ($F=\Phi$, the cumulative
distribution function of the standard normal). The shift term
$\eta_{w}(\text{\boldmath$x$})$ differs between the two treatment groups
$w\in\\{0,1\\}$. The distribution functions of the potential outcomes are
$F(h(y)-\mu(\text{\boldmath$x$}))$ for $Y(0)$ and
$F(h(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$}))$ for $Y(1)$. The
negative log-likelihood of a discrete or interval-censored observation
$(\underaccent{\bar}{\ry},\bar{y}]$ (where $\underaccent{\bar}{\ry}$ is the
lower interval bound, $\bar{y}$ is the upper) is

$\displaystyle\ell_{\text{Trafo}}(h,\mu,\tau)$ $\displaystyle=$
$\displaystyle-\log(\mathbb{P}(\underaccent{\bar}{\ry}<Y\leq\bar{y}\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w))$
$\displaystyle=$
$\displaystyle-\log(F(h(\bar{y})-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w)-F(h(\underaccent{\bar}{\ry})-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w)).$

For a continuous datum $y\in\mathbb{R}$, we obtain

$\displaystyle\ell_{\text{Trafo}}(h,\mu,\tau)=-\\{\log(F^{\prime}(h(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w))+\log(h^{\prime}(y))\\};$

details are given in Hothorn _et al._ (2018). Transformation forests apply the
model-based recursive partitioning principle and estimate $\tau$ in each node
along with the transformation function $h$ (a “nuisance” parameter) by
minimising $\ell_{\text{Trafo}}(h,\mu\equiv 0,\tau)$ (Hothorn and Zeileis,
2021b). Because $h$ contains an intercept term, the parameter $\mu$ is not
identified. We thus estimate the model under the constraint $\mu\equiv 0$.
Variable and cut-points are selected using the bivariate gradient

$\displaystyle\left.\frac{\partial\ell_{\text{Trafo}}(h,\mu\equiv
0,\tau)}{\partial(\mu,\tau)}\right|_{\mu=0,\hat{\tau}}$

This model family includes proportional odds logistic regression (for ordered
categorical, count or continuous outcomes), Box-Cox type models, Cox
proportional hazards model, Weibull proportional hazards models for discrete
and continuous outcomes, reverse time proportional hazards models relying on
Lehmann alternatives, and many more (Hothorn _et al._ , 2018). Forests for
ordinal outcomes were evaluated by Buri and Hothorn (2020), and a general
approach to “transformation forests” is described in Hothorn and Zeileis
(2021b).

Application of the ideas underlying model-based forests allows HTEs to be
estimated for such outcomes under all types of random censoring and truncation
(Korepanova _et al._ , 2020). For example, for Weibull distributed outcomes
under right censoring, $h(y)=\nu_{1}+\nu_{2}\log(y)$ is chosen for the
conditional distribution function in equation (8) (Hothorn _et al._ , 2018).

In this case, we define $Y$ as the event time, $C$ as the censoring time and
$T=\min(Y,C)$ as the observed time. For identification of
$\tau(\text{\boldmath$x$})$ under potential censoring, the following
assumption must hold (Cui _et al._ , 2022):

######  Assumption 1 (Ignorable censoring).

Censoring time $C$ is independent of survival time $Y$ conditional on
treatment indicator $W$ and covariates $X$

$(Y(0),Y(1))\perp\\!\\!\\!\perp
C\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w.$

An important special case represents the Cox proportional hazards model, where
the profile likelihood over the baseline hazard function defines the partial
log-likelihood $\ell_{\text{PL}}(\mu,\tau)$ with $\mu\equiv 0$. The scores
with respect to the constant $\mu\equiv 0$ are known as martingale residuals.
Model-based forests for such models, and extensions to time-varying prognostic
and predictive effects, are discussed in Korepanova _et al._ (2020).

###  2.4 Noncollapsibility

As mentioned in the introduction, one problem with the Cox model is that
misspecifications of prognostic effects $\mu(\text{\boldmath$x$})$ lead to
biased estimates such that the estimated hazard ratios cannot be interpreted
causally. This issue arises from the noncollapsiblity of the Cox model, the
notion of which is characterized by the fact that in these models, the mean of
the conditional effect estimates defined over covariates $X$ does not coincide
with the marginal effect over $X$. Because the noncollapsiblity of the Cox
model arises from its nonadditivity of the hazard function, models such as the
Weibull model do not suffer from this issue because they satisfy the
additivity condition. Consequently, misspecifications of prognostic effects do
not affect treatment effect estimates (Aalen _et al._ , 2015).

The noncollapsibility issue is not limited to the Cox model but also affects
members of the exponential family without identity or linear link functions.
Without adjustments, effect estimates can only be interpreted causally if
there is no treatment effect ($\tau\equiv 0$) or there are no prognostic
covariates (Daniel _et al._ , 2021).

If this is not the case, specific methods are needed; ignoring the estimation
of $\mu(\text{\boldmath$x$})$ at all and only focusing on
$\tau(\text{\boldmath$x$})$ does not solve the problem. Conditioning on
available prognostic variables is a common solution and is already applied by
model-based forests, because they estimate both the prognostic effect
$\mu(\text{\boldmath$x$})$ and $\tau(\text{\boldmath$x$})$. The ensemble of
trees used to estimate these effects provides a high degree of flexibility and
might therefore retain some of the potential complexity in the underlying
$\mu(\text{\boldmath$x$})$ to mitigate misspecification. Whether conditioning
resolves the non-collapsibility issue depends heavily on the assumption that
all prognostic variables are known which is often not the case in the real
world (Aalen _et al._ , 2015).

For members of the exponential family and the Cox model, Gao and Hastie (2022)
derived a method to account for noncollapsibility in the context of
observational data with confounding effects. While we consider the
noncollapsibility issue beyond the scope of this work, we briefly review the
work of Gao and Hastie and discuss its applicability to model-based forests in
Section A of the Supplementary Material.

##  3 Model-based forests for observational studies

In the previous section, we described model-based forests in the randomized
setting under the assumption that $\pi(\text{\boldmath$x$})=\pi$. In
observational studies in which the treatment group assignment is not under the
control of the researcher, the propensity score (and therefore, the
probability of being in the treatment group) often depends on covariates $x$

$\pi(\text{\boldmath$x$}):=\mathbb{P}(W=1\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\mathbb{E}(W\mid\text{\boldmath$X$}=\text{\boldmath$x$}).$
(9)

In this case, confounding effects could bias the estimation of treatment
effects $\tau(\text{\boldmath$x$})$, and stricter assumptions are necessary in
order to interpret $\tau(\text{\boldmath$x$})$ causally (Rosenbaum and Rubin,
1983).

######  Assumption 2 (Ignorability/Unconfoundedness).

The treatment assignment is independent of the potential outcomes conditional
on covariates $x$

$(Y(0),Y(1))\perp\\!\\!\\!\perp W\mid\text{\boldmath$X$}=\text{\boldmath$x$}.$

######  Assumption 3 (Positivity).

The propensity score $\pi(\text{\boldmath$x$})$ must be bounded away from 0
and 1

$0<\pi(\text{\boldmath$x$})=\mathbb{P}(W=1\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\mathbb{E}(W\mid\text{\boldmath$X$}=\text{\boldmath$x$})<1.$

Assumption 2 could be violated by an unmeasured confounder, while Assumption 3
could be violated if all observations in a certain group (defined via $x$) are
in the treatment group.

Dandl _et al._ (2022) showed for mean regression models that model-based
forests are not robust to confounding effects and need further adaptions to
estimate causal effects in case of observational data. One strategy for
dealing with confounding effects is the orthogonalization strategy originally
introduced by Robinson (1988), which has received considerable attention in
recent years (Chernozhukov _et al._ , 2018; Athey _et al._ , 2019; Nie and
Wager, 2021). The reformulation of the linear model

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$ $\displaystyle=$
$\displaystyle\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})W+\phi Z$ (10)

to

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$ $\displaystyle=$
$\displaystyle
m(\text{\boldmath$x$})-m(\text{\boldmath$x$})+\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})W+\phi
Z$ (11) $\displaystyle=$ $\displaystyle
m(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(W-\pi(\text{\boldmath$x$}))+\phi
Z$

given the conditional mean function

$m(\text{\boldmath$x$}):=\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})=\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})\pi(\text{\boldmath$x$}),$
(12)

motivates this approach (Dandl _et al._ , 2022).

Overall, the orthogonalization strategy consists of two steps: First, nuisance
parameters
$m(\text{\boldmath$x$})=\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
and
$\pi(\text{\boldmath$x$})=\mathbb{P}(W=1\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
are estimated. Originally, Robinson (1988) used kernel estimators, but any
machine learning method could be employed (Chernozhukov _et al._ , 2018; Nie
and Wager, 2021). Regressing $Y-\hat{m}(\text{\boldmath$x$})$ on
$W-\hat{\pi}(\text{\boldmath$x$})$ then yields unbiased estimates for
$\tau(\text{\boldmath$x$})$. Subtracting $\hat{m}(\text{\boldmath$x$})$ and
$\hat{\pi}(\text{\boldmath$x$})$ from $Y$ and $W$, respectively, partially
eliminates the association between $X$ and $Y$ and between $X$ and $W$,
respectively. The orthogonalization strategy has the distinct advantage over
other methods against confounding – such as inverse propensity weighting and
matching – that it is stable for extreme propensity scores and forgoes
stratification (Gao and Hastie, 2022).

Robinson (1988) and Chernozhukov _et al._ (2018) use parametric models to
estimate treatment effects based on residualized $W$ and $Y$, but these models
could be replaced by non-parametric or local parametric models (Nie and Wager,
2021; Wager and Athey, 2018) – such as model-based forests. For mean
regression, Dandl _et al._ (2022) adapted the orthogonalization strategy to
model-based forests. Their approach closely follows causal forests, which were
the first to combine the orthogonalization strategy with tree-based estimators
for $\tau(\text{\boldmath$x$})$.

Gao and Hastie (2022) proposed extensions of Robinson’s strategy to members of
the exponential family and the Cox model, where DINA (6) is of interest. Gao
and Hastie (2022) assume
$\tau(\text{\boldmath$x$})=\text{\boldmath$x$}^{\top}\boldsymbol{\beta}$ and
use parametric models to estimate $\tau(\text{\boldmath$x$})$, but they
conclude that non-parametric or local parametric models could be applied
instead. We review model-based forests in combination with linear models for
observational data in the next section and summarize the idea by Gao and
Hastie (2022) in Section 3.2. On this basis, we assess how the
orthogonalization strategy could be employed in model-based forests beyond
mean regression with generalized linear models and transformation models as
base models.

###  3.1 Review of Dandl et al. (2022)

As noted above, Athey _et al._ (2019) were the first to combine the
orthogonalization strategy of Robinson with tree-based estimators to estimate
$\tau(\text{\boldmath$x$})$. First, the marginal model
$m(\text{\boldmath$x$})=\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
and propensity score
$\pi(\text{\boldmath$x$})=\mathbb{E}(W\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
are estimated by regression forests. Afterwards, causal forests estimate
individual treatment effects $\tau(\text{\boldmath$x$})$ in the model

$(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=\hat{m}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))+\phi
Z$ (13)

using the “locally centered” outcomes $Y-\hat{m}(\text{\boldmath$x$})$ and
treatment indicators $W-\hat{\pi}(\text{\boldmath$x$})$.

Equation (13) shows that causal forests and model-based forests share common
foundations for mean regression. The main difference is that the splitting
scheme of model-based forests allows splitting according to heterogeneity in
both treatment and prognostic effects, whereas causal forests only split with
respect to heterogeneity in treatment effects (in equation (11),
$\mu(\text{\boldmath$x$})$ cancels out).

Dandl _et al._ (2022) identified which elements of both approaches lead to
improved performance in randomized trials and observational studies by
defining and evaluating blended versions of model-based forests and causal
forests:

  * (1)

$\text{mob}(\hat{W},\hat{Y})$, which applies model-based forests to the model

$\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=\hat{m}(\text{\boldmath$x$})+\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$})),$

i.e. after centering the treatment indicator $w$ and the outcome $Y$. Both
parameters $\tilde{\mu}$ and $\tau$ are estimated simultaneously.

  * (2)

$\text{mob}(\hat{W})$, which applies model-based forests to the model

$\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$})),$

i.e. after only centering the treatment indicator $w$ but not outcome $Y$.
Both $\mu$ and $\tau$ are estimated.

  * (3)

cfmob, a method that applies model-based forests to the model

$\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)=\hat{m}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$})),$

i.e. after only centering the treatment indicator $w$ and splitting only
according to $\hat{\tau}$. That is, only the parameters $\tau$ are estimated
in this variant.

Their blended approaches competed with the original implementations of
(uncentered) model-based forests and causal forests in an extensive simulation
study. In case of confounding, the authors identified local centering of
treatment indicator $w$ and simultaneous estimation of both predictive and
prognostic effects of the treatment indication ($\text{mob}(\hat{W})$) as the
key driver for good performance. Additionally, centering $Y$
($\text{mob}(\hat{W},\hat{Y})$) is recommended, since it further improved
performances in some cases. Splitting only according to $\hat{\tau}$ but not
$\hat{\mu}$ (cfmob) resulted in lower performance. Even for settings with
confounding, the performance of cfmob was inferior to that of uncentered
model-based forests.

###  3.2 Review of Gao and Hastie (2022)

Robinson (1988) derived the orthogonalization strategy only for semi-
parametric additive models with $Y\in\mathbb{R}$. Gao and Hastie (2022)
extended the idea to a broader class of distributions including the
exponential family and Cox’ model.

Local centering of the treatment indicator works analogously to mean
regression. First, propensity scores
$\pi(\text{\boldmath$x$})=\mathbb{P}(W\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
are estimated. The effects of the covariates $X$ on the treatment assignment
are then regressed out by subtracting $\hat{\pi}(\text{\boldmath$x$})$ from
$W$.

Orthogonalization of $Y$ is not straightforward due to the link function that
relates the linear predictor $\eta_{w}(\text{\boldmath$x$})$ in equation (5)
to the outcome $Y$. To understand how Gao and Hastie derived
$m(\text{\boldmath$x$})$ to center $Y$, we consider equation (10) as a model
of the exponential family with identity link function $g$. Now we can rewrite
equation (12) to

$\displaystyle g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}))$
$\displaystyle=$
$\displaystyle\mathbb{E}_{W}(g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)))$
$\displaystyle=$
$\displaystyle\pi(\text{\boldmath$x$})\underbrace{(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$}))}_{=\eta_{1}(\text{\boldmath$x$})}+(1-\pi(\text{\boldmath$x$}))\underbrace{\mu(\text{\boldmath$x$})}_{=\eta_{0}(\text{\boldmath$x$})}$
$\displaystyle=$
$\displaystyle\mu(\text{\boldmath$x$})+\pi(\text{\boldmath$x$})\tau(\text{\boldmath$x$})=m(\text{\boldmath$x$}).$

Similarly, we derive
$g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$}))$ for all other
distributions of the exponential family by

$m(\text{\boldmath$x$})=\pi(\text{\boldmath$x$})\eta_{1}(\text{\boldmath$x$})+(1-\pi(\text{\boldmath$x$}))\eta_{0}(\text{\boldmath$x$}).$
(14)

We can regard the estimated $m(\text{\boldmath$x$})$ as an offset in the
linear predictor

$\hat{m}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(W-\hat{\pi}(\text{\boldmath$x$})).$

Note that equation (3.2) states that (only) for the Gaussian distribution we
can directly estimate
$m(\text{\boldmath$x$})=\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
without estimating $\eta_{0}(\text{\boldmath$x$})$ and
$\eta_{1}(\text{\boldmath$x$})$. We can also derive
$\hat{m}(\text{\boldmath$x$})$ for transformation models based on the
definition of $\eta_{0}$ and $\eta_{1}$ in equation (8). As mentioned in
Section 2.4, compared to the difference in conditional means, the difference
in natural parameters additionally suffers from the noncollapsibility issue
(Greenland _et al._ , 1999). Gao and Hastie (2022) also extend the Robinson
strategy to tackle not only the confounding but also the noncollapsibility
issue for members of the exponential family (without a linear or log link
function, otherwise confounding is not an issue) and the Cox model. While the
noncollapsibility issue is beyond the scope of this work, we briefly summarize
and discuss the work of Gao and Hastie in Section A of the Supplementary
Material.

###  3.3 Novel model-based forests for observational data

As stated above, our main goal is to assess how the orthogonalization strategy
proposed for continuous outcomes could be extended to models beyond mean
regression, specifically generalized linear models and transformation models.
Based on Dandl _et al._ (2022) and Gao and Hastie (2022) we propose two
different versions of model-based forests, which should be more robust against
confounding. Following Dandl _et al._ (2022), we formulate research questions
for these versions, which we aim to answer empirically in Section 4. An
overview of all proposed versions is given in Table 1.

The first version of model-based forests directly applies Robinson’s
orthogonalization strategy: First, we estimate propensities
$\pi(\text{\boldmath$x$})$ as well as $\eta_{0}(\text{\boldmath$x$})$ and
$\eta_{1}(\text{\boldmath$x$})$ to derive $\hat{m}(\text{\boldmath$x$})$.
Then, we update the linear predictor of equation (5) by centering $W$ by
$\hat{\pi}(\text{\boldmath$x$})$ and by adding the offset
$\hat{m}(\text{\boldmath$x$})$. For generalized linear models, we obtain

$g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w))=\hat{m}(\text{\boldmath$x$})+\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))$

and for the conditional distribution function of equation (8) in case of
transformation models

$F[h(y)-\\{\hat{m}(\text{\boldmath$x$})+\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$})\\}].$

Based on the updated models, both prognostic and predictive effects
$\tilde{\mu}(\text{\boldmath$x$})$ and $\tau(\text{\boldmath$x$})$ are
simultaneously estimated by model-based forests.

In the simulation study and practical example in Sections 4 and 5, we use
regression forests to estimate $\pi(\text{\boldmath$x$})$ and gradient
boosting machines (with tailored loss functions) to estimate $\eta_{0}$ and
$\eta_{1}$. In the following, we denote this version of model-based forests as
Robinson in recognition of Robinson (1988) while model-based forests without
centering $W$ and without offset $\hat{m}(\text{\boldmath$x$})$ are called
Naive.

##### RQ 1

To what extent does centering $W$ by $\hat{\pi}(\text{\boldmath$x$})$ and
including $\hat{m}(\text{\boldmath$x$})$ as an offset affect the performance
of model-based forests in the presence of confounding?

Similar to Dandl et al. – who saw an improvement in performance when only
centering $W$ (compared to the naive model-based forests) – we define an
approach called Robinson${}_{\hat{W}}$ that applies model-based forests to
models with linear predictors

$\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$})).$

##### RQ 2

Do centered treatment indicator model-based forests perform better than
uncentered model-based forests in the presence of confounding?

##### RQ 3

Are model-based forests with centered treatment indicators relevantly
outperformed by model-based forests with $\hat{m}(\text{\boldmath$x$})$ as an
additional offset in the presence of confounding?

Table 1: Overview of proposed model-based forest versions. Method | Linear Predictor | Definitions
---|---|---
Naive | $\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})\,\,w$ |
Robinson${}_{\hat{W}}$ | $\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))$ | $\pi(\text{\boldmath$x$})=\mathbb{P}(W=1|\text{\boldmath$X$}=\text{\boldmath$x$})$
Robinson | $\hat{m}(\text{\boldmath$x$})+\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))$ | $m(\text{\boldmath$x$})=\pi(\text{\boldmath$x$})\eta_{1}(\text{\boldmath$x$})-(1-\pi(\text{\boldmath$x$}))\eta_{0}(\text{\boldmath$x$})$

##  4 Empirical evaluation

We evaluated the performance of our proposed model-based forest versions
(Table 1) in a simulation study. The study includes different outcome types,
different predictive and prognostic effects, and a varying number of
observations and covariates. Model-based forests were fitted with the
model4you R add-on package (Seibold _et al._ , 2019). Similar to Dandl _et
al._ (2022), we base our study settings on the four setups (A, B, C and D) of
Nie and Wager (2021). In addition, in Section B of the Supplementary Material,
we show the results for simulation settings first proposed by Wager and Athey
(2018) and later reused by Athey _et al._ (2019).

###  4.1 Data generating process

Given $P=\\{10,20\\}$, for Setup A, we sampled $\text{\boldmath$X$}\sim
U([0,1]^{P})$. For all other setups, we used $\text{\boldmath$X$}\sim
N(0,\mathds{1}_{P\times P})$. The treatment indicator was binomially
distributed with
$W\mid\text{\boldmath$X$}=\text{\boldmath$x$}\sim\operatorname{B}(1,\pi(\text{\boldmath$x$}))$.
The propensity function $\pi(\text{\boldmath$x$})$ differed for the four
considered setups:

$\displaystyle\pi(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\pi_{A}(x_{1},x_{2})=\max\\{0.1,\min\\{\sin(\pi
x_{1}x_{2}),1-0.1\\}\\}\\\ \pi_{B}\equiv 0.5\\\
\pi_{C}(x_{2},x_{3})=1/(1+\exp(x_{2}+x_{3}))\\\
\pi_{D}(x_{1},x_{2})=1/(1+\exp(-x_{1})+\exp(-x_{2})).\end{array}\right.$

$\pi(\text{\boldmath$x$})\equiv 0.5$ in Setup B implies a randomized study.
The treatment effect function $\tau(\cdot)$ and the prognostic effect function
$\mu(\cdot)$ also differed between the setups

$\displaystyle\tau(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\tau_{A}(x_{1},x_{2})=(x_{1}+x_{2})/2\\\
\tau_{B}(x_{1},x_{2})=x_{1}+\log(1+\exp(x_{2}))\\\ \tau_{C}\equiv 1\\\
\tau_{D}(x_{1},x_{2},x_{3},x_{4},x_{5})=\max\\{x_{1}+x_{2}+x_{3},0\\}-\max\\{x_{4}+x_{5},0\\}.\end{array}\right.$
$\displaystyle\mu(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\mu_{A}(x_{1},x_{2},x_{3},x_{4},x_{5})=\sin(\pi
x_{1}x_{2})+2(x_{3}-0.5)^{2}+x_{4}+0.5x_{5}\\\
\mu_{B}(x_{1},x_{2},x_{3})=\max\\{x_{1}+x_{2},x_{3},0\\}+\max\\{x_{4}+x_{5},0\\}\\\
\mu_{C}(x_{1},x_{2},x_{3})=2\log(1+\exp(x_{1}+x_{2}+x_{3}))\\\
\mu_{D}(x_{1},x_{2},x_{3},x_{4},x_{5})=(\max\\{x_{1}+x_{2}+x_{3},0\\}+\max\\{x_{4}+x_{5},0\\})/2.\end{array}\right.$

Setup A has extensive confounding that must be eliminated before estimating an
easily predictable treatment effect function $\tau(\text{\boldmath$x$})$.
Setup B needs no confounding adjustment for reliable estimation of $\tau$.
Although Setup C contains strong confounding, the propensity score function is
easier to estimate than the prognostic effect, while the treatment effect is
constant. In Setup D, the treatment and control arms are unrelated, and
therefore, learning the conditional expected outcomes of both arms jointly is
not beneficial (Nie and Wager, 2021; Dandl _et al._ , 2022).

We studied four different simulation models

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{N}(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-0.5),1)$
(18a) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{B}(1,\text{expit}(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-0.5)))$
(18b) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{M}\text{with}\log(O(y_{k}\mid\text{\boldmath$x$},w))=\vartheta_{k}-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})(w-0.5)$
(18c) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{W}\text{with}\log(H(y\mid\text{\boldmath$x$},w))=2\log(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})(w-0.5)$
(18d)

Model (18a) is a normal linear regression model, model (18b) is a binary
logistic regression model, model (18c) is a 4-nomial model with log-odds
function
$\vartheta_{k}-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})(w-0.5)$ with
threshold parameters $\vartheta_{k}=\text{logit}(k/4)$ for $k=1,2,3$, and
model (18d) is a Weibull model with log-cumulative hazard function
$2\log(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})(w-0.5)$. We added
$50\,\%$ random right-censoring to the Weibull-generated data. Additionally,
we applied a Cox proportional hazards model to the Weibull data to determine
if the performance of model-based forests degrades when the forests do not
take the true underlying model as their base model.

Due to $w-0.5$ in all scenarios, half of the (negative) predictive effect
$\tau(\text{\boldmath$x$})$ was added to the prognostic effect. We refer to
the implied scenario – where one variable which is both prognostic (impact in
$\mu(\text{\boldmath$x$})$) and predictive (impact in
$\tau(\text{\boldmath$x$})$) exists – as overlay. Apart from Setup C in which
the treatment effect is constant and independent of any covariate, overlay was
present for all scenarios.

Like Dandl _et al._ (2022), we compared all study settings and outcome types
for a varying number of samples $N\in\\{800,1600\\}$ and dimensions
$P\in\\{10,20\\}$. All model-based forests were grown with the same
hyperparameter options specified in Section 7. We used random forests as
implemented in the grf package to estimate $\pi(\text{\boldmath$x$})$ for
centering $W$ (Tibshirani _et al._ , 2021). To estimate
$\eta_{0}(\text{\boldmath$x$})$ and $\eta_{1}(\text{\boldmath$x$})$ to derive
$\hat{m}(\text{\boldmath$x$})$, we relied on different tree-based estimators
depending on the outcome type. For normally distributed outcomes (models
(18a)), we used grf regression forests (Tibshirani _et al._ , 2021). For all
other outcomes, we relied on gradient boosting machines (with adapted loss
functions) as implemented in mboost and gbm (Hothorn _et al._ , 2021b;
Greenwell _et al._ , 2020). The employed distribution varied depending on the
outcome type.

In accordance with Dandl _et al._ (2022), we evaluated the models with respect
to the mean squared error
$\mathbb{E}_{\text{\boldmath$X$}}\\{(\hat{\tau}(\text{\boldmath$X$})-\tau(\text{\boldmath$X$}))^{2}\\}$
on a test sample of size $1000$. The results are shown in Figure 1 and were
statistically analyzed by means of a normal linear mixed model with a log-
link. The model explained the estimated mean squared error for
$\hat{\tau}(\text{\boldmath$x$})$ by a four-way interaction of the data
generating process, sample size $N$, dimension $P$, and random forest variant.
We estimated the mean squared error ratios between different model-based
forest versions according to the two research questions stated in Section 3.3.
The corresponding tables are given in Tables 2 to 4.

Figure 1: Model-based forest results for the empirical study (Section 4), Cox means a Cox model applied to the Weibull data. For the Weibull and Cox model, treatment effects $\tau(\text{\boldmath$x$})$ are estimated as conditional log hazard ratios. Direct comparison of model-based forests without centering (Naive), model-based forests with local centering according to Robinson (1988) of $Y$ and $W$ (originally proposed) (Robinson) or only of $W$ (Robinson${}_{\widehat{W}}$).  Table 2: Results of RQ 1 for the experimental setups in Section 4. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of the Naive model-based forests, and cells printed in italics indicate an inferior reference.  |  |  |  |  Mean squared error ratio for RQ 1: Robinson vs. Naive
---|---|---|---|---
DGP | N | P |  | Normal | Binomial | Multinomial | Weibull | Cox
Setup A | 800 | 10 |  | 0.465 (0.421, 0.512) | 1.173 (1.045, 1.316) | 0.690 (0.629, 0.758) | 0.672 (0.609, 0.742) | 0.712 (0.650, 0.781)
|  | 20 |  | 0.396 (0.359, 0.438) | 1.161 (1.014, 1.330) | 0.600 (0.540, 0.666) | 0.605 (0.547, 0.669) | 0.654 (0.596, 0.718)
| 1600 | 10 |  | 0.414 (0.362, 0.474) | 1.042 (0.892, 1.216) | 0.582 (0.512, 0.662) | 0.580 (0.508, 0.663) | 0.589 (0.519, 0.669)
|  | 20 |  | 0.341 (0.295, 0.395) | 0.898 (0.751, 1.075) | 0.503 (0.428, 0.591) | 0.471 (0.405, 0.548) | 0.495 (0.430, 0.570)
Setup B | 800 | 10 |  | 0.643 (0.607, 0.681) | 1.021 (0.929, 1.121) | 0.868 (0.830, 0.907) | 0.692 (0.653, 0.733) | 0.703 (0.668, 0.739)
|  | 20 |  | 0.658 (0.625, 0.693) | 0.977 (0.894, 1.067) | 0.906 (0.870, 0.943) | 0.716 (0.681, 0.754) | 0.731 (0.699, 0.763)
| 1600 | 10 |  | 0.603 (0.557, 0.653) | 0.981 (0.873, 1.101) | 0.852 (0.803, 0.903) | 0.657 (0.607, 0.710) | 0.656 (0.612, 0.702)
|  | 20 |  | 0.588 (0.544, 0.636) | 0.912 (0.811, 1.026) | 0.869 (0.822, 0.917) | 0.648 (0.603, 0.697) | 0.653 (0.614, 0.695)
Setup C | 800 | 10 |  | 0.153 (0.144, 0.163) | 0.474 (0.432, 0.520) | 0.250 (0.233, 0.268) | 0.174 (0.160, 0.189) | 0.176 (0.162, 0.191)
|  | 20 |  | 0.156 (0.147, 0.166) | 0.359 (0.322, 0.400) | 0.219 (0.204, 0.235) | 0.157 (0.142, 0.172) | 0.161 (0.146, 0.177)
| 1600 | 10 |  | 0.154 (0.139, 0.171) | 0.361 (0.316, 0.412) | 0.260 (0.238, 0.284) | 0.181 (0.160, 0.206) | 0.187 (0.165, 0.212)
|  | 20 |  | 0.157 (0.142, 0.173) | 0.300 (0.256, 0.351) | 0.215 (0.195, 0.238) | 0.147 (0.129, 0.169) | 0.152 (0.133, 0.174)
Setup D | 800 | 10 |  | 0.818 (0.801, 0.835) | 1.109 (1.037, 1.187) | 0.996 (0.968, 1.026) | 1.036 (1.008, 1.065) | 1.085 (1.057, 1.113)
|  | 20 |  | 0.851 (0.835, 0.867) | 1.126 (1.058, 1.199) | 1.054 (1.028, 1.082) | 1.055 (1.029, 1.081) | 1.099 (1.075, 1.124)
| 1600 | 10 |  | 0.783 (0.762, 0.805) | 1.075 (0.985, 1.175) | 0.994 (0.957, 1.032) | 0.968 (0.934, 1.004) | 1.029 (0.995, 1.063)
|  | 20 |  | 0.803 (0.783, 0.824) | 1.131 (1.046, 1.223) | 1.016 (0.983, 1.051) | 1.021 (0.989, 1.053) | 1.076 (1.046, 1.108)
Table 3: Results of RQ 2 for the experimental setups in Section 4. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of the Naive model-based forests, and cells printed in italics indicate an inferior reference.  |  |  |  |  Mean squared error ratio for RQ 2: Robinson${}_{\hat{W}}$ vs. Naive
---|---|---|---|---
DGP | N | P |  | Normal | Binomial | Multinomial | Weibull | Cox
Setup A | 800 | 10 |  | 1.029 (0.910, 1.164) | 0.820 (0.729, 0.922) | 1.259 (1.142, 1.388) | 0.924 (0.820, 1.042) | 0.844 (0.752, 0.947)
|  | 20 |  | 1.060 (0.933, 1.204) | 0.784 (0.679, 0.905) | 1.282 (1.144, 1.437) | 0.935 (0.825, 1.060) | 0.835 (0.740, 0.942)
| 1600 | 10 |  | 1.126 (0.953, 1.330) | 0.915 (0.781, 1.072) | 1.370 (1.194, 1.571) | 1.067 (0.911, 1.250) | 1.015 (0.870, 1.184)
|  | 20 |  | 1.163 (0.970, 1.395) | 0.887 (0.726, 1.084) | 1.302 (1.086, 1.561) | 1.063 (0.881, 1.283) | 0.994 (0.831, 1.188)
Setup B | 800 | 10 |  | 1.555 (1.468, 1.647) | 0.980 (0.892, 1.077) | 1.152 (1.102, 1.205) | 1.445 (1.363, 1.531) | 1.423 (1.353, 1.496)
|  | 20 |  | 1.520 (1.444, 1.600) | 1.024 (0.938, 1.119) | 1.104 (1.060, 1.150) | 1.396 (1.327, 1.469) | 1.368 (1.309, 1.430)
| 1600 | 10 |  | 1.658 (1.530, 1.796) | 1.019 (0.907, 1.144) | 1.174 (1.107, 1.245) | 1.524 (1.409, 1.648) | 1.525 (1.424, 1.634)
|  | 20 |  | 1.700 (1.574, 1.837) | 1.097 (0.975, 1.233) | 1.151 (1.090, 1.216) | 1.542 (1.435, 1.657) | 1.532 (1.440, 1.629)
Setup C | 800 | 10 |  | 1.871 (1.743, 2.009) | 1.377 (1.243, 1.526) | 1.577 (1.456, 1.708) | 2.331 (2.128, 2.553) | 2.388 (2.182, 2.614)
|  | 20 |  | 2.081 (1.944, 2.226) | 1.294 (1.138, 1.470) | 1.718 (1.588, 1.859) | 2.565 (2.318, 2.839) | 2.611 (2.363, 2.886)
| 1600 | 10 |  | 1.774 (1.573, 2.001) | 2.619 (2.288, 2.999) | 1.759 (1.594, 1.942) | 2.198 (1.920, 2.517) | 2.141 (1.874, 2.446)
|  | 20 |  | 1.817 (1.629, 2.026) | 1.800 (1.512, 2.144) | 1.675 (1.494, 1.877) | 2.541 (2.203, 2.932) | 2.566 (2.228, 2.956)
Setup D | 800 | 10 |  | 1.136 (1.113, 1.161) | 0.910 (0.851, 0.974) | 0.992 (0.964, 1.021) | 0.916 (0.890, 0.942) | 0.883 (0.860, 0.906)
|  | 20 |  | 1.098 (1.077, 1.120) | 0.898 (0.844, 0.956) | 0.942 (0.918, 0.966) | 0.909 (0.886, 0.932) | 0.881 (0.861, 0.901)
| 1600 | 10 |  | 1.147 (1.114, 1.180) | 0.950 (0.871, 1.037) | 0.994 (0.958, 1.032) | 0.965 (0.929, 1.001) | 0.923 (0.892, 0.954)
|  | 20 |  | 1.126 (1.097, 1.157) | 0.890 (0.823, 0.961) | 0.972 (0.940, 1.005) | 0.922 (0.893, 0.952) | 0.888 (0.862, 0.914)
Table 4: Results of RQ 3 for the experimental setups in Section 4. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of Robinson${}_{\hat{W}}$, and cells printed in italics indicate an inferior reference.  |  |  |  |  Mean squared error ratio for RQ 3: Robinson vs. Robinson${}_{\hat{W}}$
---|---|---|---|---
DGP | N | P |  | Normal | Binomial | Multinomial | Weibull | Cox
Setup A | 800 | 10 |  | 0.972 (0.859, 1.099) | 1.220 (1.085, 1.373) | 0.794 (0.720, 0.876) | 1.082 (0.959, 1.220) | 1.185 (1.056, 1.329)
|  | 20 |  | 0.944 (0.831, 1.072) | 1.276 (1.105, 1.472) | 0.780 (0.696, 0.874) | 1.070 (0.944, 1.212) | 1.197 (1.061, 1.351)
| 1600 | 10 |  | 0.888 (0.752, 1.049) | 1.093 (0.933, 1.281) | 0.730 (0.637, 0.838) | 0.937 (0.800, 1.098) | 0.985 (0.844, 1.149)
|  | 20 |  | 0.860 (0.717, 1.030) | 1.127 (0.922, 1.378) | 0.768 (0.641, 0.921) | 0.941 (0.780, 1.135) | 1.006 (0.841, 1.203)
Setup B | 800 | 10 |  | 0.643 (0.607, 0.681) | 1.020 (0.929, 1.121) | 0.868 (0.830, 0.907) | 0.692 (0.653, 0.733) | 0.703 (0.669, 0.739)
|  | 20 |  | 0.658 (0.625, 0.692) | 0.976 (0.894, 1.067) | 0.906 (0.869, 0.943) | 0.716 (0.681, 0.754) | 0.731 (0.699, 0.764)
| 1600 | 10 |  | 0.603 (0.557, 0.654) | 0.981 (0.874, 1.102) | 0.852 (0.803, 0.903) | 0.656 (0.607, 0.710) | 0.656 (0.612, 0.702)
|  | 20 |  | 0.588 (0.544, 0.635) | 0.912 (0.811, 1.026) | 0.869 (0.822, 0.917) | 0.649 (0.603, 0.697) | 0.653 (0.614, 0.695)
Setup C | 800 | 10 |  | 0.534 (0.498, 0.574) | 0.726 (0.655, 0.804) | 0.634 (0.586, 0.687) | 0.429 (0.392, 0.470) | 0.419 (0.383, 0.458)
|  | 20 |  | 0.481 (0.449, 0.514) | 0.773 (0.680, 0.878) | 0.582 (0.538, 0.630) | 0.390 (0.352, 0.431) | 0.383 (0.346, 0.423)
| 1600 | 10 |  | 0.564 (0.500, 0.636) | 0.382 (0.333, 0.437) | 0.569 (0.515, 0.628) | 0.455 (0.397, 0.521) | 0.467 (0.409, 0.534)
|  | 20 |  | 0.550 (0.494, 0.614) | 0.555 (0.467, 0.661) | 0.597 (0.533, 0.669) | 0.393 (0.341, 0.454) | 0.390 (0.338, 0.449)
Setup D | 800 | 10 |  | 0.880 (0.861, 0.899) | 1.099 (1.027, 1.175) | 1.008 (0.979, 1.037) | 1.092 (1.061, 1.123) | 1.133 (1.104, 1.163)
|  | 20 |  | 0.911 (0.893, 0.929) | 1.113 (1.046, 1.185) | 1.062 (1.035, 1.089) | 1.101 (1.073, 1.128) | 1.136 (1.110, 1.162)
| 1600 | 10 |  | 0.872 (0.848, 0.898) | 1.052 (0.964, 1.148) | 1.006 (0.969, 1.044) | 1.037 (0.999, 1.076) | 1.084 (1.048, 1.121)
|  | 20 |  | 0.888 (0.865, 0.912) | 1.124 (1.040, 1.215) | 1.029 (0.995, 1.064) | 1.085 (1.050, 1.120) | 1.126 (1.094, 1.160)

###  4.2 Results

The results for the normal distribution coincide with the results obtained by
Dandl _et al._ (2022) summarized in Section 3.1. To some degree, they also
hold for the other distributions. The boxplots are not directly comparable
between different data generating processes because of different signal-to-
noise ratios. In general, a more informative outcome (binary $<$ ordered $<$
right-censored $<$ exact normal), more data (higher $N$), and less noise
(lower $P$) leads to better results. Using a Cox model compared to a Weibull
model (last two rows of Figure 1) did not lead to a major decrease in
performance, although knowledge of the true functional form of the
transformation function did not enter the Cox modeling process.

For Setup A, model-based forests without centering (Naive) were unable to cope
with complex confounding, but solely centering of the treatment indicator
(Robinson${}_{\hat{W}}$) was valuable. Additionally adding
$\hat{m}(\text{\boldmath$x$})$ as an offset (Robinson) did not further improve
the results for the normal, binomial, and Weibull distributions, but an
improvement was observed for the multinomial distribution.

For Setup B, the Robinson forests performed slightly better in disentangling
the more complicated prognostic and predictive effects compared to Naive and
Robinson${}_{\hat{W}}$ model-based forests. An exception is the binomial
model: without overlay, Robinson${}_{\hat{W}}$ forests performed similarly to
Robinson forests.

In Setup C, over all distributions, uncentered model-based forests (Naive)
failed to overcome the strong confounding effect and therefore did not provide
accurate estimates for the treatment effect. The performance was fundamentally
improved by centering the treatment indicator (Robinson${}_{\hat{W}}$) and was
further improved by additionally adding $\hat{m}(\text{\boldmath$x$})$ as an
offset (Robinson).

In Setup D – with unrelated treatment and control arms – all methods had a
higher mean squared error than in the other setups, as jointly modeling the
expected conditional outcomes for both arms has no benefit. Apart from the
normal distributions, Robinson forests were inferior to the
Robinson${}_{\hat{W}}$ and Naive model-based forests.

The empirical evidence of our simulation study can be summarized as follows:
If confounding was present, model-based forests performed better when
centering $W$ by $\hat{\pi}(\text{\boldmath$x$})$ (Robinson${}_{\hat{W}}$)
compared to not centering $W$ (Naive). Adding $\hat{m}(\text{\boldmath$x$})$
as an offset (Robinson) further improved the performance – especially in cases
with very strong confounding.

##  5 Effect of Riluzole on progression of ALS

Amyotrophic lateral sclerosis (ALS) is a progressive nervous system disease
causing loss of muscle control. The status of the disease as well as the rate
of progression is commonly evaluated by the ALS functional rating scale
(ALSFRS) (Brooks _et al._ , 1996; Cedarbaum _et al._ , 1999). Here, physical
abilities such as speaking, handwriting, and walking are assessed and rated on
a scale from 0 (inability) to 4 (normal ability). In 1995, the FDA approved
the first drug to manage and slow progression of ALS, named Riluzole. The
largest database for study results on the effect of Riluzole offers the Pooled
Resource Open-Access Clinical Trials (PROACT) database – initiated by the non-
profit organization Prize4Life (http://www.prize4life.org). The data comes
from different randomized and observational studies not disclosed in the data.
Thus, the assumption of random treatment assignment is quite hard to justify
in an analysis. Patient characteristics and treatment group sizes might vary
greatly between the centers, which affect both the probability of receiving
treatment as well as the outcome. To account for these potential confounding
effects, we compared the treatment effects estimated by the naive model-based
forests to the ones estimated with local centering by Robinson. As in Section
4, we use random forests to estimate the propensity scores to center $W$ and
gradient boosting machines (with adapted loss functions) to estimate the
values of the linear predictors $\eta_{0}(\text{\boldmath$x$})$ and
$\eta_{1}(\text{\boldmath$x$})$ to center $Y$. Model-based forests, random
forests, and gradient boosting machines rely on the hyperparameter values
stated in Section 7. As for Seibold _et al._ (2018) and Korepanova _et al._
(2020), 16 phase II and phase III randomized trials and one observational
study from the PROACT database serve as a training dataset. We analyze the
effect of Riluzole with respect to two outcome variables: survival time and
the handwriting ability score approximately six months after treatment – an
item of the ALSFRS. We omitted observations with missing outcome values. As
splitting variables, Seibold _et al._ (2018) used demographic, medical
history, and family history data, which were informative in the sense that not
more than half of their values were missing.

Figure 2: Kaplan-Meier curves of survival probability for both treatment arms.

###  5.1 Survival Time

The dataset for the survival time contains 3306 observations and 18
covariates. Of the 3306 observations, 2199 received Riluzole. Because very few
patients had event times that exceed those of the others by a factor of two,
we artificially censored five observations with (censoring or event) times of
more than 750 days. The Kaplan-Meier estimates of survival probabilities for
both treatment arms of the preprocessed dataset are shown in Figure 2.
Overall, the estimated survival curves are very close to each other, and the
treated group has only a slight survival advantage compared to the untreated
group. As a base model, we use a Cox proportional hazards model. We compared
treatment effects from two approaches: the naive uncentered model-based
forests (Naive) and the model-based forest with Robinson’s orthogonalization
(Robinson).

####  5.1.1 Personalized models

For the naive model-based forests, the underlying Cox proportional hazards
base model for the survival outcome $T$ was, on the hazard scale,

$\lambda(t)=\lambda_{0}(t)\exp(\mu+\tau w)$

Because $\lambda_{0}(t)$ contains an intercept term, $\mu$ is not identified
(and was constraint to $\mu\equiv 0$). The treatment effect $\tau$ is the log-
hazard ratio of the treated versus untreated patients and our aim is to
replace a constant marginal effect $\tau$ with a heterogenuous (and thus
conditional) log-hazard ratio $\tau(\text{\boldmath$x$})$ and, simultaneously,
to estimate prognostic effects $\mu(\text{\boldmath$x$})$.

For Robinson’s strategy, we first centered the treatment indicator $W$ by
estimating the propensity scores
$\pi(\text{\boldmath$x$})=\mathbb{P}(W\mid\text{\boldmath$X$}=\text{\boldmath$x$})$
using a regression forest. Figure 3 compares the distributions of estimated
propensity scores (left) and of the estimated centered treatment
$W-\hat{\pi}(\text{\boldmath$x$})$ (Robinson’s strategy, right), both obtained
from regression forests. We can already see a decent overlap of propensity
scores in the two treatment arms without centering, but the overlap increases
if the strategy by Robinson was applied.

Figure 3: Distribution of estimated propensities
$\hat{\pi}(\text{\boldmath$x$})$ (left) and estimated propensities of the
centered treatment indicators (right, Robinson’s strategy) as estimated by
regression forests for the two treatment groups.

In addition to centering $W$, Robinson’s strategy requires the estimation of
$m(\text{\boldmath$x$})$ to use as an offset (see Section 3). As in Section 4,
we used gradient boosting machines (with the negative log partial likelihood
of the Cox proportional hazards model as a loss) to estimate the natural
parameters $\eta_{0}(\text{\boldmath$x$})$ and $\eta_{1}(\text{\boldmath$x$})$
for the control and treatment group, respectively (Friedman, 2001). The offset
$m(\text{\boldmath$x$})$ for each observation is equal to the sum of natural
parameter estimates weighted by $\hat{\pi}(\text{\boldmath$x$})$ (see equation
(14)). The final base model for model-based forests using Robinson’s
orthogonalization is

$\lambda_{R}(t)=\lambda_{0}(t)\exp(\mu+\tau(w-\hat{\pi}(\text{\boldmath$x$}))+\hat{m}(\text{\boldmath$x$})).$

Figure 4: Kernel density estimates of the personalized treatment estimates for
the naive model-based forest (Naive) and for the model-based forest with
Robinson orthogonalization (Robinson).

####  5.1.2 Model-based forests

The corresponding base models serve as an input for the model-based forests to
estimate personalized effects of Riluzole. Figure 4 compares the kernel
density estimates of $\tau(\text{\boldmath$x$})$ for each forest version
(Naive and Robinson). The naive approach reveals that on average the treatment
reduced the hazard compared to no treatment, whereas the model-based forest
with centering according to Robinson obtained weaker effects of Riluzole with
more mass centered around 0.

A meta-analysis of previous studies by Andrews _et al._ (2020), also yielded a
mixed picture: only eight of the 15 studies meeting their inclusion criteria
showed a statistically significant increase of median survival time due to
Riluzole.

Over all strategies, for both approaches there were some patients for which
Riluzole was estimated to increase the hazard. The dependency plots in Figures
S. 4 and S. 5 in the Supplementary Material provide indications of the
characteristics of the group of harmed individuals. For example, both the
naive and centering approach agree that for patients with atrophy or
fasciculation, Riluzole intake would increase the hazard. The estimated
effects differed most between the uncentered forest (Naive) and the
orthogonalized forest (Robinson) for the covariate sex (Figure S. 4 (c)), the
covariate of whether patients swallow, and for the covariate specifying
whether cases in the same generation exist (Figure S. 5 (f) and (i)).

For the variables time onset treatment, age, height and weakness the
dependency plots (Figure S. 5 (a), (d), (e) and Figure S. 6 (g)) of the Naive
forest agree with the ones of Seibold _et al._ (2018): for middle-aged people
with a longer time between disease onset and start of treatment, lower height,
and no weakness, the treatment appears to be more beneficial. By considering
confounding effects due to orthogonalization (Robinson forests), these effects
diminished. For Korepanova _et al._ (2020) the effect of Riluzole was also
rather weak and showed low heterogeneity across covariates.

###  5.2 Handwriting Ability Score

The dataset for the handwriting ability score – an ordinal outcome with five
categories – contains 2538 observations and 58 covariates. Besides the
covariate age, all covariates had missing values (but less than 50 % of the
values were missing per variable enforced by the preprocessing step stated at
the beginning of this section). Of the 2538 observations, 1754 received
Riluzole, and 784 did not. Figure 5 displays the frequency of the ability
scores for both treatment groups. Most of the patients have an ability score
of 3 or 4 (normal ability); only a few have ability scores less than 2. Note
that the plot shows the conditional proportions given the treatment indicator.
We chose a proportional odds logistic regression model as a base model for the
model-based forests – once without further adaptions (Naive), and once
parameterized with centered $W$ and with an offset (Robinson).

In addition to the handwriting ability score after six months, the ability
score values at treatment start are also available. In the following, we
denote $Y_{6}$ as the handwriting score after six months and $Y_{0}$ as the
handwriting score at the beginning of the treatment period. To account for the
ability level at treatment start, $Y_{0}$ served as an additional splitting
variable for both model-based forests (Naive and Robinson) and was included in
$\mathbf{X}$. The alluvial plot in Figure 6 breaks down the change in each
ability class over six months. Overall, for most patients, the handwriting
ability remained constant over the six months or worsened slightly. Rarely,
patients experienced a progression to both extremes (0 to 4, or 4 to 0). These
results hold regardless of whether patients received Riluzole or not.

Figure 5: Relative frequency distribution plot of the handwriting ability
score ($Y_{6}$) (left) and of changes of the handwriting ability score over
six months ($Y_{6}-Y_{0}$) (right) for both treatment arms. Frequencies were
calculated relative to the treatment indicator.  Figure 6: Alluvial plot of
the progression of the handwriting ability score over six months for both
treatment arms.

####  5.2.1 Personalized models

The proportional odds logistic regression model for the naive model-based
forests is defined as (Agresti, 2002; Venables and Ripley, 2002)

$\operatorname{logit}(\mathbb{P}(Y_{6}\leq
k|\text{\boldmath$X$}=\text{\boldmath$x$},W=w,Y_{0}=y_{0}))=\vartheta_{k}(\text{\boldmath$x$},y_{0})-\tau(\text{\boldmath$x$},y_{0})w$

with $k\in\\{0,...,3\\}$ as the ordinal ability score classes. The parameters
$\vartheta_{k}$ are increasing thresholds, depending on covariates $x$ and the
initial score $y_{0}$. Due to the proportional odds assumption, the treatment
effect $\tau(\text{\boldmath$x$},y_{0})$ is the same for all scores $k$.
Negative $\tau(\text{\boldmath$x$},y_{0})$ indicate a negative effect of
Riluzole, as treated patients are expected to have a higher odds of low
writing ability scores compared to untreated patients.

As for the survival forest, we used regression forests to estimate propensity
scores $\pi(\text{\boldmath$x$},y_{0})$ and a gradient boosting machine (with
adapted loss functions for the proportional odds model) to estimate the
natural parameters $\eta_{0}(\text{\boldmath$x$},y_{0})$ and
$\eta_{1}(\text{\boldmath$x$},y_{0})$. The personalized model for the model-
based forest with Robinson orthogonalization was specified as

$\operatorname{logit}(\mathbb{P}(Y_{6}\leq
k|\text{\boldmath$X$}=\text{\boldmath$x$},W=w,Y_{0}=y_{0}))=\vartheta_{k}(\text{\boldmath$x$},y_{0})-[\hat{m}(\text{\boldmath$x$},y_{0})+\tau(\text{\boldmath$x$},y_{0})\\{w-\hat{\pi}(\text{\boldmath$x$},y_{0})\\}]$

with $\hat{m}(\text{\boldmath$x$},y_{0})$ as defined in equation (14).

Figure 7 compares the estimated treatment indicators with $W$ as the outcome
in the random forest without centering (left), with
$(W-\hat{\pi}(\text{\boldmath$x$},y_{0}))$ as the outcome in the random forest
(right). Before centering, there is a lack of overlap of the propensity
scores; the distribution of $\hat{\pi}$ for the control group is bimodal, and
the distribution for the treatment group is heavily left-skewed. After
centering, the distributions of the estimated
$W-\hat{\pi}(\text{\boldmath$x$},y_{0})$ for the treatment groups move closer
together and have a similar unimodal shape. However, there is still a lack of
overlap of the groups, which indicates that important covariates to explain
the remaining heterogeneity in the two treatment groups seem to be missing.

Figure 7: Estimates returned by the regression forest (rf) for
orthogonalization of the treatment indicator: left for original $W$ as an
outcome in the rf such that it estimates propensity scores
$\pi(\text{\boldmath$x$},y_{0})$; right for the centered treatment indicator
$W-\hat{\pi}(\text{\boldmath$x$},y_{0})$ as an outcome in the rf.

####  5.2.2 Model-based forests

The proportional odds logistic regression models served as a base model for
the (Naive and Robinson) model-based forests to derive personalized treatment
effects. Figure 8 displays the kernel density estimates of
$\tau(\text{\boldmath$x$},y_{0})$ for each forest version (Naive and
Robinson).

Figure 8: Kernel density estimates of the personalized treatment estimates for
the naive model-based forest (Naive) vs. the forest with Robinson
orthogonalized (Robinson).

Both random forests estimate on average a negative effect of Riluzole. Naive
model-based forest estimated on average a log-odds of $\bar{\tau}=-0.08$,
which indicates that treated patients have a 0.08 points higher log-odds for
low writing scores than untreated patients. The distribution of
$\hat{\tau}(\text{\boldmath$x$},y_{0})$ for the model-based forest relying on
the Robinson orthogonalization is slightly shifted to the left
($\bar{\tau}=-0.10$). For a larger subgroup of patients, the naive approach
estimates a negative effect of Riluzole
($-1\leq\tau(\text{\boldmath$x$},y_{0})\leq-0.5$), meaning that patients
receiving treatment with Riluzole have higher odds of low writing scores than
untreated patients. According to the dependency plots (Figures S. 6 to S. 11
in the Supplementary Material), this subgroup could be identified as having
the low initial ability scores (left side of Figure S. 6 (a)). For all other
splitting variables, the distributions of estimated treatment effects are very
similar.

##  6 Discussion and outlook

HTE estimation is a challenging problem, especially for observational studies
and even more when the outcome cannot be modeled by a linear model. In this
work, we investigated several versions of model-based forests for the
estimation of potentially complex HTEs $\tau(\text{\boldmath$x$})$ based on
observational data with various outcome types based on the orthogonalization
strategy by Robinson (Robinson, 1988). These investigations suggest the
following workflow for model-based forests: (1) estimate propensities
$\pi(\text{\boldmath$x$})$ using some machine learning procedure (binary
random forests are a good default), (2) center the treatment indicator
$w-\hat{\pi}(\text{\boldmath$x$})$ for each observation, (3) setup an
appropriate model for the outcome conditioning on the centered treatment and –
if possible – add an offset for centering $Y$, (4) use model-based forests to
estimate predictive and prognostic effects $\tau(\text{\boldmath$x$})$ and
$\mu(\text{\boldmath$x$})$ simultaneously. Notably,
$\tau(\text{\boldmath$x$})$ is the CATE only in specific models, especially a
linear or log-linear model. We demonstrate these steps by estimating the
individual effects of Riluzole for ALS patients using survival times and
ordinal ability scores as outcomes.

Our work still leaves open questions for example how model-based forests
perform for survival data for which the censoring procedure is not randomized
but depends on $X$, or how ($k$-fold) cross-fitting influences the
performance, where only one part of the data is used to estimate nuisance
parameters and the other part to estimate $\tau(\text{\boldmath$x$})$
(Chernozhukov _et al._ , 2018). We leave investigations to these questions to
future research.

Last but not least, we want to emphasize that all approaches for estimating
HTEs – including those presented in this work – rely on strong and typically
untestable assumptions. For example, for models beyond mean regression,
$\hat{\tau}(\text{\boldmath$x$})$ cannot be expected to be robust against
missing covariates or other violations of model assumptions due to non-
collapsibility. Consequently, results from these approaches in practical
applications should be evaluated with the utmost caution, reservation, and
humility.

##  7 Computational details

For all computations, we used R version 4.1.1 (R Core Team, 2022), with the
following add-on packages: model4you (Seibold _et al._ , 2021), trtf (Hothorn,
2021), partykit (Hothorn and Zeileis, 2021a), grf (Tibshirani _et al._ ,
2021), mboost (Hothorn _et al._ , 2021b), and gbm (Greenwell _et al._ , 2020).

Model-based forests were always grown with M $=500$ trees (model4you::pmforest
default) with a minimum node size of node $=14$, number of chosen variables
per split mtry $=P$, and subsampling. These settings were also used by Dandl
_et al._ (2022). Transformation forests implemented in the trtf package fitted
the Weibull transformation forests of Section 4 (Hothorn, 2021; Hothorn and
Zeileis, 2021b).

Propensity scores $\pi(\text{\boldmath$x$})$ were estimated with grf (honest)
regression forests with 125 trees, a minimum node size of 5, and subsampling.
Natural parameters $\eta_{0}(\text{\boldmath$x$})$ and
$\eta_{1}(\text{\boldmath$x$})$ and probability of not being censored were
estimated with gradient boosting machines implemented in the mboost or gbm
packages. The used maximum tree depth was 2 (default of mboost::blackboost),
and a loss function that differed depending on the outcome type was also
employed (Hothorn _et al._ , 2021b; Greenwell _et al._ , 2020).

Ratios and confidence intervals presented in Table 2 were calculated using
generalized linear mixed models of the glmmTMB package (Magnusson _et al._ ,
2021). Post-hoc inference relied on the multcomp package (Hothorn _et al._ ,
2021a).

All study settings are available in a dedicated R package called htesim (Dandl
and Hothorn, 2021). It is published on Github:
https://github.com/dandls/htesim.

## Funding

TH received funding from the Swiss National Science Foundation, Grant No.
200021_184603, Horizon 2020 Research and Innovation Programme of the European
Union under grant agreement number 681094, and is supported by the Swiss State
Secretariat for Education, Research and Innovation (SERI) under contract
number 15.0137.

## References

  * Aalen _et al._ (2015) Aalen OO, Cook RJ, Røysland K (2015).  “Does Cox Analysis of a Randomized Survival Study Yield a Causal Treatment Effect?”  _Lifetime Data Analysis_ , 21(4), 579–593.  10.1007/s10985-015-9335-y.
  * Agresti (2002) Agresti A (2002).  _Categorical Data Analysis_.  John Wiley & Sons, Inc.  10.1002/0471249688.
  * Andrews _et al._ (2020) Andrews JA, Jackson CE, Heiman-Patterson TD, Bettica P, Brooks BR, Pioro EP (2020).  “Real-world Evidence of Riluzole Effectiveness in Treating Amyotrophic Lateral Sclerosis.”  _Amyotrophic Lateral Sclerosis and Frontotemporal Degeneration_ , 21(7–8), 509–518.  10.1080/21678421.2020.1771734.
  * Athey _et al._ (2019) Athey S, Tibshirani J, Wager S (2019).  “Generalized Random Forests.”  _The Annals of Statistics_ , 47(2), 1148–1178.  10.1214/18-aos1709.
  * Brooks _et al._ (1996) Brooks B, Sanjak M, Ringel S, England J, Brinkmann J (1996).  “The Amyotrophic Lateral Sclerosis Functional Rating Scale: Assessment of Activities of Daily Living in Patients With Amyotrophic Lateral Sclerosis.”  _Archives of Neurology_ , 53(2), 141–147.  10.1001/archneur.1996.00550020045014.
  * Buri and Hothorn (2020) Buri M, Hothorn T (2020).  “Model-based Random Forests for Ordinal Regression.”  _International Journal of Biostatistics_ , 16(2), 20190063\.  10.1515/ijb-2019-0063.
  * Cedarbaum _et al._ (1999) Cedarbaum JM, Stambler N, Malta E, Fuller C, Hilt D, Thurmond B, Nakanishi A (1999).  “The ALSFRS-R: A Revised ALS Functional Rating Scale that Incorporates Assessments of Respiratory Function.”  _Journal of the Neurological Sciences_ , 169(1), 13–21.  https://doi.org/10.1016/S0022-510X(99)00210-5.
  * Chapfuwa _et al._ (2021) Chapfuwa P, Assaad S, Zeng S, Pencina MJ, Carin L, Henao R (2021).  “Enabling Counterfactual Survival Analysis with Balanced Representations.”  CHIL ’21, pp. 133–145. Association for Computing Machinery, New York, NY, USA.  ISBN 9781450383592.  10.1145/3450439.3451875.
  * Chernozhukov _et al._ (2018) Chernozhukov V, Chetverikov D, Demirer M, Duflo E, Hansen C, Newey W, Robins J (2018).  “Double/Debiased Machine Learning for Treatment and Structural Parameters.”  _The Econometrics Journal_ , 21(1), C1–C68.  10.1111/ectj.12097.
  * Cui _et al._ (2022) Cui Y, Kosorok MR, Sverdrup E, Wager S, Ruoqing (2022).  “Estimating Heterogeneous Treatment Effects with Right-Censored Data via Causal Survival Forests.”  _Technical report_ , arXiv 2001.09887 v3.  URL https://arxiv.org/abs/2001.09887.
  * Curth _et al._ (2021) Curth A, Lee C, van der Schaar M (2021).  “SurvITE: Learning Heterogeneous Treatment Effects from Time-to-Event Data.”  _Technical report_ , arXiv 2110.14001.  URL https://arxiv.org/abs/2110.14001.
  * Dandl and Hothorn (2021) Dandl S, Hothorn T (2021).  _htesim: Conducting Extensive Simulation Studies of Heterogeneous Treatment Effect Estimation_.  R package version 0.0.0.9000, URL https://github.com/susanne-207/htesim.
  * Dandl _et al._ (2022) Dandl S, Hothorn T, Seibold H, Sverdrup E, Wager S, Zeileis A (2022).  “What Makes Forest-Based Heterogeneous Treatment Effect Estimators Work?”  _Technical report_ , arXiv 2206.10323.  URL https://arxiv.org/abs/2206.10323.
  * Daniel _et al._ (2021) Daniel R, Zhang J, Farewell D (2021).  “Making Apples from Oranges: Comparing Noncollapsible Effect Estimators and Their Standard Errors after Adjustment for Different Covariate Sets.”  _Biometrical Journal_ , 63(3), 528–557.  https://doi.org/10.1002/bimj.201900297.
  * Fokkema _et al._ (2018) Fokkema M, Smits N, Zeileis A, Hothorn T, Kelderman H (2018).  “Detecting Treatment-Subgroup Interactions in Clustered Data with Generalized Linear Mixed-Effects Model Trees.”  _Behavior Research Methods_ , 50(6), 2016–2034.  10.3758/s13428-017-0971-x.
  * Foster _et al._ (2011) Foster JC, Taylor J, Ruberg S (2011).  “Subgroup Identification from Randomized Clinical Trial Data.”  _Statistics in Medicine_ , 30(24), 2867–2880.  10.1002/sim.4322.
  * Friedman (2001) Friedman JH (2001).  “Greedy Function Approximation: A Gradient Boosting Machine.”  _The Annals of Statistics_ , 29(5), 1189–1232.  10.1214/aos/1013203451.
  * Gao and Hastie (2022) Gao Z, Hastie T (2022).  “Estimating Heterogeneous Treatment Effects for General Responses.”  _Technical report_ , arXiv 2103.04277 v4.  URL https://arxiv.org/abs/2103.04277.
  * Greenland (1996) Greenland S (1996).  “Absence of Confounding Does Not Correspond to Collapsibility of the Rate Ratio or Rate Difference.”  _Epidemiology_ , 7, 498–501.
  * Greenland _et al._ (1999) Greenland S, Pearl J, Robins JM (1999).  “Confounding and Collapsibility in Causal Inference.”  _Statistical Science_ , 14(1), 29–46.  10.1214/ss/1009211805.
  * Greenwell _et al._ (2020) Greenwell B, Boehmke B, Cunningham J, Developers G (2020).  _gbm: Generalized Boosted Regression Models_.  R package version 2.1.8, URL https://CRAN.R-project.org/package=gbm.
  * Henderson _et al._ (2018) Henderson NC, Louis TA, Rosner GL, Varadhan R (2018).  “Individualized Treatment Effects with Censored Data via Fully Nonparametric Bayesian Accelerated Failure Time Models.”  _Biostatistics_ , 21(1), 50–68.  ISSN 1465-4644.  10.1093/biostatistics/kxy028.
  * Hill (2011) Hill JL (2011).  “Bayesian Nonparametric Modeling for Causal Inference.”  _Journal of Computational and Graphical Statistics_ , 20(1), 217–240.  10.1198/jcgs.2010.08162.
  * Hothorn (2021) Hothorn T (2021).  _trtf: Transformation Trees and Forests_.  R package version 0.3-8, URL https://CRAN.R-project.org/package=trtf.
  * Hothorn _et al._ (2021a) Hothorn T, Bretz F, Westfall P (2021a).  _multcomp: Simultaneous Inference in General Parametric Models_.  R package version 1.4-17, URL https://CRAN.R-project.org/package=multcomp.
  * Hothorn _et al._ (2021b) Hothorn T, Bühlmann P, Kneib T, Schmid M, Hofner B (2021b).  _mboost: Model-Based Boosting_.  R package version 2.9-5, URL https://CRAN.R-project.org/package=mboost.
  * Hothorn _et al._ (2006) Hothorn T, Hornik K, Zeileis A (2006).  “Unbiased Recursive Partitioning: A Conditional Inference Framework.”  _Journal of Computational and Graphical Statistics_ , 15(3), 651–674.  10.1198/106186006x133933.
  * Hothorn _et al._ (2004) Hothorn T, Lausen B, Benner A, Radespiel-Tröger M (2004).  “Bagging Survival Trees.”  _Statistics in Medicine_ , 23(1), 77–91.  10.1002/sim.1593.
  * Hothorn _et al._ (2018) Hothorn T, Möst L, Bühlmann P (2018).  “Most Likely Transformations.”  _Scandinavian Journal of Statistics_ , 45(1), 110–134.  10.1111/sjos.12291.
  * Hothorn and Zeileis (2021a) Hothorn T, Zeileis A (2021a).  _partykit: A Toolkit for Recursive Partytioning_.  R package version 1.2-15, URL https://R-Forge.R-project.org/projects/partykit/.
  * Hothorn and Zeileis (2021b) Hothorn T, Zeileis A (2021b).  “Predictive Distribution Modelling Using Transformation Forests.”  _Journal of Computational and Graphical Statistics_ , 14, 144–148.  10.1080/10618600.2021.1872581.
  * Hu _et al._ (2020) Hu L, Gu C, Lopez MJ, Ji J, Wisnivesky J (2020).  “Estimation of Causal Effects of Multiple Treatments in Observational Studies with a Binary Outcome.”  _Statistical Methods in Medical Research_ , 29, 3218–3234.  10.1177/0962280220921909.
  * Hu _et al._ (2021) Hu L, Ji J, Li F (2021).  “Estimating Heterogeneous Survival Treatment Effect in Observational Data using Machine Learning.”  _Statistics in Medicine_ , 40, 4691–4713.  10.1002/sim.9090.
  * Ishwaran _et al._ (2008) Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS (2008).  “Random Survival Forests.”  _The Annals of Applied Statistics_ , 2(3), 841–860.  10.1214/08-aoas169.
  * Katzman _et al._ (2018) Katzman J, Shaham U, Cloninger A, Bates J, Jiang T, Kluger Y (2018).  “DeepSurv: Personalized Treatment Recommender System using a Cox Proportional Hazards Deep Neural Network.”  _BMC Medical Research Methodology_ , 18, 24.  10.1186/s12874-018-0482-1.
  * Korepanova _et al._ (2020) Korepanova N, Seibold H, Steffen V, Hothorn T (2020).  “Survival Forests under Test: Impact of the Proportional Hazards Assumption on Prognostic and Predictive Forests for ALS Survival.”  _Statistical Methods in Medical Research_ , 29(5), 1403–1419.  10.1177/0962280219862586.
  * Künzel _et al._ (2019) Künzel SR, Sekhon JS, Bickel PJ, Yu B (2019).  “Metalearners for Estimating Heterogeneous Treatment Effects Using Machine Learning.”  _Proceedings of the National Academy of Sciences of the United States of America_ , 116(10), 4156–4165.  10.1073/pnas.1804597116.
  * Li _et al._ (2022) Li L, Levine RA, Fan J (2022).  “Causal Effect Random Forest of Interaction Trees for Learning Individualized Treatment Regimes with Multiple Treatments in Observational Studies.”  _Stat_ , 11(1), e457.  10.1002/sta4.457.
  * Lin and Jeon (2006) Lin Y, Jeon Y (2006).  “Random Forests and Adaptive Nearest Neighbors.”  _Journal of the American Statistical Association_ , 101(474), 578–590.  10.1198/016214505000001230.
  * Lu _et al._ (2018) Lu M, Sadiq S, Feaster DJ, Ishwaran H (2018).  “Estimating Individual Treatment Effect in Observational Data Using Random Forest Methods.”  _Journal of Computational and Graphical Statistics_ , 27(1), 209–219.  10.1080/10618600.2017.1356325.
  * Magnusson _et al._ (2021) Magnusson A, Skaug H, Nielsen A, Berg C, Kristensen K, Maechler M, van Bentham K, Bolker B, Brooks M (2021).  _glmmTMB: Generalized Linear Mixed Models using Template Model Builder_.  R package version 1.1.2.3, URL https://CRAN.R-project.org/package=glmmTMB.
  * Meinshausen (2006) Meinshausen N (2006).  “Quantile Regression Forests.”  _Journal of Machine Learning Research_ , 7, 983–999.  10.1007/s10994-014-5452-1.
  * Nie and Wager (2021) Nie X, Wager S (2021).  “Quasi-Oracle Estimation of Heterogeneous Treatment Effects.”  _Biometrika_ , 108, 299–319.  10.1093/biomet/asaa076.
  * Powers _et al._ (2018) Powers S, Qian J, Jung K, Schuler A, Shah NH, Hastie T, Tibshirani R (2018).  “Some Methods for Heterogeneous Treatment Effect Estimation in High Dimensions.”  _Statistics in Medicine_ , 37(11), 1767–1787.  10.1002/sim.7623.
  * R Core Team (2022) R Core Team (2022).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL https://www.R-project.org/.
  * Robinson (1988) Robinson PM (1988).  “Root-N-Consistent Semiparametric Regression.”  _Econometrica_ , 56(4), 931–954.  10.2307/1912705.
  * Rosenbaum and Rubin (1983) Rosenbaum PR, Rubin DB (1983).  “The Central Role of the Propensity Score in Observational Studies for Causal Effects.”  _Biometrika_ , 70(1), 41–55.
  * Rubin (1974) Rubin DB (1974).  “Estimating causal effects of treatments in randomized and nonrandomized studies.”  _Journal of Educational Psychology_ , 66, 688–701.
  * Rubin (2005) Rubin DB (2005).  “Causal Inference Using Potential Outcomes: Design, Modeling, Decisions.”  _Journal of the American Statistical Association_ , 100(469), 322–331.
  * Seibold _et al._ (2016) Seibold H, Zeileis A, Hothorn T (2016).  “Model-Based Recursive Partitioning for Subgroup Analyses.”  _International Journal of Biostatistics_ , 12(1), 45–63.  10.1515/ijb-2015-0032.
  * Seibold _et al._ (2018) Seibold H, Zeileis A, Hothorn T (2018).  “Individual Treatment Effect Prediction for Amyotrophic Lateral Sclerosis Patients.”  _Statistical Methods in Medical Research_ , 27(10), 3104–3125.  10.1177/0962280217693034.
  * Seibold _et al._ (2019) Seibold H, Zeileis A, Hothorn T (2019).  “model4you: An R Package for Personalised Treatment Effect Estimation.”  _Journal of Open Research Software_ , 7(17), 1–6.  10.5334/jors.219.
  * Seibold _et al._ (2021) Seibold H, Zeileis A, Hothorn T (2021).  _model4you: Stratified and Personalised Models Based on Model-Based Trees and Forests_.  R package version 0.9-7, URL https://R-Forge.R-project.org/projects/partykit/.
  * Shalit _et al._ (2017) Shalit U, Johansson FD, Sontag D (2017).  “Estimating Individual Treatment Effect: Generalization Bounds and Algorithms.”  _Technical report_ , arXiv 1606.03976.  URL https://arxiv.org/abs/1606.03976.
  * Su _et al._ (2018) Su X, Peña AT, Liu L, Levine RA (2018).  “Random Forests of Interaction Trees for Estimating Individualized Treatment Effects in Randomized Trials.”  _Statistics in Medicine_ , 37(17), 2547–2560.  10.1002/sim.7660.
  * Tabib and Larocque (2020) Tabib S, Larocque D (2020).  “Non-parametric Individual Treatment Effect Estimation for Survival Data with Random Forests.”  _Bioinformatics_ , 36, 629–636.  10.1093/bioinformatics/btz602.
  * Tibshirani _et al._ (2021) Tibshirani J, Athey S, Sverdrup E, Wager S (2021).  _grf: Generalized Random Forests_.  R package version 2.0.2, URL https://github.com/grf-labs/grf.
  * Venables and Ripley (2002) Venables WN, Ripley BD (2002).  _Modern Applied Statistics with S_.  Springer New York.  10.1007/978-0-387-21706-2.
  * Wager and Athey (2018) Wager S, Athey S (2018).  “Estimation and Inference of Heterogeneous Treatment Effects Using Random Forests.”  _Journal of the American Statistical Association_ , 113(523), 1228–1242.  10.1080/01621459.2017.1319839.
  * Zeileis _et al._ (2008) Zeileis A, Hothorn T, Hornik K (2008).  “Model-Based Recursive Partitioning.”  _Journal of Computational and Graphical Statistics_ , 17(2), 492–514.  10.1198/106186008x319331.

##  Appendix A Noncollapsibility

As mentioned in Section 2.4, for members of the exponential family without an
identity or linear link function the marginal and conditional treatment
effects are not collapsible. This means that the mean of the conditional
treatment effects given a covariate are not equal to the marginal treatment
effect estimate over the same covariate (Greenland _et al._ , 1999). This
happens if the covariate conditioned on is associated with the outcome of
interest. Caution is necessary on multiple stages of the estimation process of
$\tau(\text{\boldmath$x$})$ as soon as we condition on other covariates, for
example, because these covariates are assumed to be sufficient to control for
confounding (Daniel _et al._ , 2021).

In case of Robinson’s orthogonalization, misspecification of
$m(\text{\boldmath$x$})$ translates into biased estimators for
$\tau(\text{\boldmath$x$})$, even under randomized treatments. This also
applies if one ignores the estimation of $\mu(\text{\boldmath$x$})$ at all and
only concentrates on $\tau(\text{\boldmath$x$})$. This is not the case for the
linear model (identity link function) since misspecifications are absorbed in
the additive error term and do not influence the estimation of
$\tau(\text{\boldmath$x$})$ (Gao and Hastie, 2022).

###  A.1 Review Gao and Hastie (2022)

Gao and Hastie (2022) extended the orthogonalization strategy of Robinson
(1988) to improve robustness to both confounding and noncollapsibility. The
authors propose

$a(\text{\boldmath$x$})=\frac{\pi(\text{\boldmath$x$})\frac{\partial\gamma(\eta_{1}(\text{\boldmath$x$}))}{\partial\eta}}{\pi(\text{\boldmath$x$})\frac{\partial\gamma(\eta_{1}(\text{\boldmath$x$}))}{\partial\eta}+(1-\pi(\text{\boldmath$x$}))\frac{\partial\gamma(\eta_{0}(\text{\boldmath$x$}))}{\partial\eta}}$
(19)

and

$\nu(\text{\boldmath$x$})=a(\text{\boldmath$x$})n_{1}(\text{\boldmath$x$})+(1-a(\text{\boldmath$x$}))n_{0}(\text{\boldmath$x$})$

instead of $\pi(\text{\boldmath$x$})$ (equation (9)) and
$m(\text{\boldmath$x$})$ (equation (14)), respectively, where $\gamma(\eta)$
denotes the inverse of the canonical link function. Its derivative is equal to
the variance function of the exponential family. Therefore,
$a(\text{\boldmath$x$})$ is larger if an observation is likely to be treated
(which also holds for Robinson’s orthogonalization) or if the response
variance is higher under treatment compared to no treatment. As a consequence
of the latter, the influence of spuriously influential natural parameter
values is reduced for more robustness to misspecifications (Gao and Hastie,
2022).

For Gaussian responses, $a(\text{\boldmath$x$})=\pi(\text{\boldmath$x$})$ and
$\nu(\text{\boldmath$x$})=m(\text{\boldmath$x$})$ holds, while for other
distributions the terms differ. For example, for Bernoulli distributed $Y$,
the closed form $a(\text{\boldmath$x$})$ is

$a(\text{\boldmath$x$})=\frac{\pi(\text{\boldmath$x$})}{\pi(\text{\boldmath$x$})+(1-\pi(\text{\boldmath$x$}))\frac{p_{0}(\text{\boldmath$x$})(1-p_{0}(\text{\boldmath$x$}))}{p_{1}(\text{\boldmath$x$})(1-p_{1}(\text{\boldmath$x$}))}}$
(20)

where
$p_{w}(\text{\boldmath$x$})=\mathbb{P}(Y=1|\text{\boldmath$X$}=\text{\boldmath$x$},W=w)$.

The noncollapsibility issue is not only present for distributions of the
exponential family. Also the Cox model suffers from noncollapsibility
(Greenland, 1996; Aalen _et al._ , 2015). This is in contrast to accelerated
failure time models (such as the Weibull proportional hazards model), which
can be rewritten as location-scale models and therefore are indeed collapsible
(Aalen _et al._ , 2015). For the Cox model, Gao and Hastie remark that with
knowledge of the baseline hazard function and without censoring, the
cumulative hazard function follows an exponential distribution. For the
exponential distribution, $a(\text{\boldmath$x$})$ and
$\nu(\text{\boldmath$x$})$ are equal to $\pi(\text{\boldmath$x$})$ and
$m(\text{\boldmath$x$})$ (Gao and Hastie, 2022).

In case of random censoring, the probability of not being censored under both
treatment arms needs to be considered for the estimation of
$a(\text{\boldmath$x$})$ and $\nu(\text{\boldmath$x$})$

$a(\text{\boldmath$x$})=\frac{\pi(\text{\boldmath$x$})\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=1)}{\pi(\text{\boldmath$x$})\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=1)+(1-\pi(\text{\boldmath$x$}))\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=0)}$ (21)

$\nu(\text{\boldmath$x$})=a(\text{\boldmath$x$})\eta_{1}(\text{\boldmath$x$})+(1-a(\text{\boldmath$x$}))\eta_{0}(\text{\boldmath$x$}).$
(22)

The nuisance parameter $a(\text{\boldmath$x$})$ is larger if an observation is
likely to be treated or likely to be not censored. Consequently, the influence
of likely to be not censored observations for the estimation of
$\tau(\text{\boldmath$x$})$ is increased. Above’s $a(\text{\boldmath$x$})$ and
$\nu(\text{\boldmath$x$})$ guarantee protection to misspecified nuisance
parameter if the baseline hazard is known. If it is unknown and the partial
likelihood is used – this is not guaranteed. Despite this lack of guarantee,
Gao and Hastie, 2022, obtained promising results in their simulation study
(Gao and Hastie, 2022).

###  A.2 Strategies against confounding and noncollapsibility

An interesting question is if replacing $\hat{\pi}(\text{\boldmath$x$})$ and
$\hat{m}(\text{\boldmath$x$})$ by $\hat{a}(\text{\boldmath$x$})$ and
$\hat{\nu}(\text{\boldmath$x$})$, respectively, also helps to additionally
tackle noncollapsibility when applying model-based forests. We can update the
linear predictor for model-based forests in case of generalized linear models
to

$g(\mathbb{E}(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w))=\hat{\nu}(\text{\boldmath$x$})+\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{a}(\text{\boldmath$x$})).$

Gao and Hastie additionally derived estimators for $a(\text{\boldmath$x$})$
and $\nu(\text{\boldmath$x$})$ for the Cox model which – compared to the
Weibull model – is not collapsible. For the Cox model, the natural parameter
of equation (8) could be updated to

$\eta_{w}(\text{\boldmath$x$})=\hat{\nu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{a}(\text{\boldmath$x$}))$

with $a(\text{\boldmath$x$})$ and $\nu(\text{\boldmath$x$})$ as defined in
equations (21) and (22).

We call this version of model-based forests in the following Gao approach.
Before we apply model-based forests, we need to estimate
$\pi(\text{\boldmath$x$})$, $\eta_{0}(\text{\boldmath$x$})$,
$\eta_{1}(\text{\boldmath$x$})$ as well as
$\frac{\partial\upsilon(\eta_{1}(\text{\boldmath$x$}))}{\partial\eta}$ for
exponential families and $\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=w)$ for Cox models. As in Section
3.3, we state some research questions that are empirically inspected in the
upcoming section.

##### RQ 4:

How do model-based forests centered according to Gao and Hastie (Gao) perform
compared to model-based forest with Robinson strategy (Robinson) for the
simulation settings of Section 4?

Similar to RQ 2, we could solely center $W$ by $a(\text{\boldmath$x$})$
without including an offset. We call this approach Gao${}_{\hat{W}}$ in the
following.

##### RQ 5:

How do model-based forest with solely centered $W$ by
$\hat{a}(\text{\boldmath$x$})$ (Gao${}_{\hat{W}}$) perform compared to model-
based forests with solely centered $W$ by $\hat{\pi}(\text{\boldmath$x$})$
Robinson${}_{\hat{W}}$ for the simulation study settings of Section 4?

Table S. 1: Updated overview of proposed model-based forest versions (Table 1)
for observational data.

Method | Linear Predictor | Definitions
---|---|---
Naive | $\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})\,\,w$ |
Robinson${}_{\hat{W}}$ | $\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))$ | $\pi(\text{\boldmath$x$})=\mathbb{P}(W=1|\text{\boldmath$X$}=\text{\boldmath$x$})$
Robinson | $\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{\pi}(\text{\boldmath$x$}))+\hat{m}(\text{\boldmath$x$})$ | $m(\text{\boldmath$x$})=\pi(\text{\boldmath$x$})\eta_{1}(\text{\boldmath$x$})-(1-\pi(\text{\boldmath$x$}))\eta_{0}(\text{\boldmath$x$})$
Gao${}_{\hat{W}}$ | $\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{a}(\text{\boldmath$x$}))$ | $a(\text{\boldmath$x$})=\frac{\pi(\text{\boldmath$x$})\frac{\partial\gamma(\eta_{1}(\text{\boldmath$x$}))}{\partial\eta}}{\pi(\text{\boldmath$x$})\frac{\partial\gamma(\eta_{1}(\text{\boldmath$x$}))}{\partial\eta}+(1-\pi(\text{\boldmath$x$}))\frac{\partial\gamma(\eta_{0}(\text{\boldmath$x$}))}{\partial\eta}}$
Gao | $\tilde{\mu}(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})(w-\hat{a}(\text{\boldmath$x$}))+\hat{\nu}(\text{\boldmath$x$})$ | $\nu(\text{\boldmath$x$})=a(\text{\boldmath$x$})n_{1}(\text{\boldmath$x$})+(1-a(\text{\boldmath$x$}))n_{0}(\text{\boldmath$x$})$

  * •

Note: for the Cox model
$a(\text{\boldmath$x$})=\frac{\pi(\text{\boldmath$x$})\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=1)}{\pi(\text{\boldmath$x$})\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=1)+(1-\pi(\text{\boldmath$x$}))\mathbb{P}(C\geq
Y|\text{\boldmath$X$}=\text{\boldmath$x$},W=0)}$ is used.

###  A.3 Data-generating process

To investigate the research questions of Section A.2., we compared the
performance of model-based forests with Gao’s strategy proposed in Section A.2
(Gao and Gao${}_{\hat{W}}$) to model-based forests with Robinson’s startegy
(Robinson and Robinson${}_{\hat{W}}$) for settings A, B, C, D described in
Section 4. Because we expect that the strategy of Gao is especially valuable
for settings with misspecified prognostic effect, e.g. because prognostic
covariates are missing, we additionally created Setup A’ from Setup A by
removing covariate $\text{\boldmath$X$}_{3}$ from the training data.
Therefore, the DGP of Setup A and Setup A’ are identical, the only difference
being that the training data did not contain $\text{\boldmath$X$}_{3}$
although $\text{\boldmath$X$}_{3}$ affects the prognostic effect.

Because the normal linear model and Weibull model are collapsible and Gao’s
strategy is equal to Robinson’s strategy (Sections 2.4 and A.1), we applied
our proposed approaches based on Gao and Hastie (2022) only to the binomial
model and the Cox model. Transformation models such as the proportional odds
model for multinomial data were not covered by the authors.

We used the same model-based forest parameter setup and evaluation scheme as
in Section 4.

Figure S. 1: Model-based forest results for the empirical study (Section 4), Cox means a Cox model applied to the Weibull data. For the Cox model, treatment effects $\tau(\text{\boldmath$x$})$ are estimated as conditional log hazard ratios. Direct comparison of model-based forests without centering (Naive), model-based forests with local centering according to Robinson (1988) or Gao and Hastie (2022) of $Y$ and $W$ (originally proposed) (Robinson, Gao) or only of $W$ (Robinson${}_{\widehat{W}}$, Gao${}_{\widehat{W}}$).  Table S. 2: Results of RQ 4 for the experimental setups in Section 4. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of Robinson forests, cells printed in italics indicate an inferior reference.  |  |  |  | Mean squared error ratio for RQ 4: Gao vs. Robinson
---|---|---|---|---
DGP | N | P |  | Binomial | Cox
Setup A | 800 | 10 |  | 1.258 (1.152, 1.373) | 1.203 (1.077, 1.344)
|  | 20 |  | 1.307 (1.180, 1.449) | 1.307 (1.170, 1.461)
| 1600 | 10 |  | 1.067 (0.933, 1.220) | 1.121 (0.947, 1.326)
|  | 20 |  | 1.183 (1.009, 1.388) | 1.155 (0.955, 1.398)
Setup A’ | 800 | 10 |  | 1.201 (1.105, 1.304) | 1.140 (1.011, 1.285)
|  | 20 |  | 1.354 (1.233, 1.488) | 1.272 (1.127, 1.435)
| 1600 | 10 |  | 1.047 (0.915, 1.200) | 1.055 (0.895, 1.243)
|  | 20 |  | 1.184 (1.014, 1.382) | 1.114 (0.911, 1.362)
Setup B | 800 | 10 |  | 1.042 (0.958, 1.134) | 0.984 (0.920, 1.052)
|  | 20 |  | 0.987 (0.909, 1.073) | 0.906 (0.853, 0.963)
| 1600 | 10 |  | 0.987 (0.885, 1.100) | 0.977 (0.889, 1.074)
|  | 20 |  | 0.926 (0.824, 1.042) | 0.922 (0.845, 1.006)
Setup C | 800 | 10 |  | 1.388 (1.263, 1.524) | 1.417 (1.261, 1.592)
|  | 20 |  | 1.616 (1.448, 1.804) | 1.401 (1.228, 1.598)
| 1600 | 10 |  | 1.276 (1.104, 1.476) | 1.360 (1.146, 1.615)
|  | 20 |  | 1.485 (1.255, 1.758) | 1.400 (1.163, 1.686)
Setup D | 800 | 10 |  | 0.996 (0.939, 1.057) | 0.916 (0.889, 0.943)
|  | 20 |  | 0.965 (0.913, 1.020) | 0.925 (0.902, 0.949)
| 1600 | 10 |  | 0.964 (0.890, 1.044) | 0.910 (0.875, 0.946)
|  | 20 |  | 0.948 (0.884, 1.015) | 0.907 (0.877, 0.938)
Table S. 3: Results of RQ 5 for the experimental setups in Section 4. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of Robinson${}_{\hat{W}}$ forests, cells printed in italics indicate an inferior reference.  |  |  |  |  Mean squared error ratio for RQ 5: Gao${}_{\hat{W}}$ vs. Robinson${}_{\hat{W}}$
---|---|---|---|---
DGP | N | P |  | Binomial | Cox
Setup A | 800 | 10 |  | 1.299 (1.168, 1.445) | 1.127 (0.986, 1.288)
|  | 20 |  | 1.425 (1.255, 1.618) | 1.190 (1.038, 1.366)
| 1600 | 10 |  | 1.162 (1.009, 1.339) | 1.110 (0.940, 1.310)
|  | 20 |  | 1.339 (1.128, 1.589) | 1.144 (0.944, 1.386)
Setup A’ | 800 | 10 |  | 1.261 (1.139, 1.397) | 1.096 (0.952, 1.263)
|  | 20 |  | 1.427 (1.264, 1.610) | 1.195 (1.033, 1.382)
| 1600 | 10 |  | 1.096 (0.950, 1.264) | 1.060 (0.896, 1.255)
|  | 20 |  | 1.305 (1.101, 1.548) | 1.114 (0.906, 1.370)
Setup B | 800 | 10 |  | 0.988 (0.904, 1.079) | 1.005 (0.959, 1.053)
|  | 20 |  | 0.959 (0.883, 1.042) | 1.037 (0.995, 1.081)
| 1600 | 10 |  | 0.947 (0.849, 1.056) | 0.968 (0.910, 1.031)
|  | 20 |  | 0.905 (0.811, 1.009) | 0.982 (0.929, 1.038)
Setup C | 800 | 10 |  | 1.228 (1.141, 1.323) | 1.636 (1.561, 1.715)
|  | 20 |  | 1.658 (1.524, 1.804) | 1.585 (1.510, 1.664)
| 1600 | 10 |  | 0.716 (0.660, 0.776) | 1.552 (1.437, 1.677)
|  | 20 |  | 1.272 (1.149, 1.408) | 1.588 (1.481, 1.702)
Setup D | 800 | 10 |  | 1.011 (0.948, 1.079) | 1.004 (0.973, 1.037)
|  | 20 |  | 0.981 (0.923, 1.042) | 0.987 (0.960, 1.016)
| 1600 | 10 |  | 0.969 (0.891, 1.054) | 1.027 (0.987, 1.069)
|  | 20 |  | 0.970 (0.899, 1.048) | 1.003 (0.968, 1.039)

###  A.4 Results

For Setup A, solely centering $W$ by $\hat{a}(\text{\boldmath$x$})$
(Gao${}_{\hat{W}}$) achieved better results than additionally adding the
offset $\hat{\nu}(\text{\boldmath$x$})$ (Gao). Model-based forests with
Robinson’s strategy (Robinson, Robinson${}_{\hat{W}}$) overall performed
better than model-based forests with Gao’s strategy (Gao, Gao${}_{\hat{W}}$).
Surpressing $X_{3}$ in the training dataset (Setup A’), did not deteriorate
the performance of all methods such that the ranking of methods was retained.

For Setup B, model-based forests centered by Gao and Robinson model-based
forests performed akin for binary outcomes. Also Robinson${}_{\hat{W}}$ and
Gao${}_{\hat{W}}$ model-based forests achieved similar performance.

In Setup C, Gao’s strategy for the Cox and logistic regression model overall
fare worse than Robinson’s strategy. In Setup D, Gao${}_{\hat{W}}$ forests
performed as good as Robinson${}_{\hat{W}}$ forests for the Cox and logistic
regression models. Notably, for the Cox model, Gao forests outperformed
Robinson forests.

Overall, the orthogonalization strategy of Gao for the exponential family –
that aims at addressing the noncollapsibility issue – did not perform as well
as expected. Our expectation was that the strategy would reduce the effect of
overfitting the marginal effect $\hat{m}(\text{\boldmath$x$})$ on the
treatment effect estimate. Overall, however, the estimation of additional
nuisance parameters tended to worsen the performance results on average – at
least for the binomial model. For the Cox model, Gao’s strategy, which
additionally takes the probability for not getting censored into account, did
not worsen performance. Further experiments are necessary in which the
censoring probability is not constant but depends on covariates $x$.

##  Appendix B Empirical evaluation based on Wager and Athey (2018)

We evaluated the performance of our proposed model-based forest versions also
with the study setting of Wager and Athey (2018), which were later reused by
Athey _et al._ (2019). Given uniformly distributed covariates
$\text{\boldmath$X$}\sim U([0,1]^{P})$ of dimensionality $P\in\\{10,20\\}$ and
a binomially distributed treatment indicator
$W\mid\text{\boldmath$X$}=\text{\boldmath$x$}\sim\operatorname{B}(1,\pi(\text{\boldmath$x$}))$,
the propensity function $\pi(\cdot)$ either did or did not depend on $x$

$\displaystyle\pi(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\pi\equiv
0.5\\\ \pi(x_{1})=\nicefrac{{1}}{{4}}(1+\beta_{2,4}(x_{1}))\\\
\pi(x_{3})=\nicefrac{{1}}{{4}}(1+\beta_{2,4}(x_{3}))\\\
\pi(x_{4})=\nicefrac{{1}}{{4}}(1+\beta_{2,4}(x_{4}))\end{array}\right.$

where $\beta_{2,4}$ is the $\beta$-density with shape $2$ and scale $4$. The
probability $\pi\equiv 0.5$ indicates no confounding and thus a randomized
trial. The treatment effect function $\tau(\cdot)$ was either $0$ (no
treatment effect) or depended on a smooth interaction function of $x_{1}$ and
$x_{2}$

$\displaystyle\tau(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\tau\equiv
0\\\
\tau(x_{1},x_{2})=\prod_{p=1,2}\left(1+\left(1+\exp\left(-20\left(x_{p}-\nicefrac{{1}}{{3}}\right)\right)\right)^{-1}\right).\end{array}\right.$

The prognostic effect function $\mu(\cdot)$ was either $0$ (no prognostic
effect) or linear in $x_{1}$ or $x_{3}$

$\displaystyle\mu(\text{\boldmath$x$})=\left\\{\begin{array}[]{l}\mu\equiv
0\\\ \mu(x_{1})=2x_{1}-1\\\ \mu(x_{3})=2x_{3}-1.\end{array}\right.$

We studied four different simulation models

$\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{N}(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w,1)$
(26a) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{B}(1,\text{expit}(\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w))$
(26b) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{M}\text{with}\log(O(y_{k}\mid\text{\boldmath$x$},w))=\vartheta_{k}-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w$
(26c) $\displaystyle(Y\mid\text{\boldmath$X$}=\text{\boldmath$x$},W=w)\sim$
$\displaystyle\operatorname{W}\text{with}\log(H(y\mid\text{\boldmath$x$},w))=2\log(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w$
(26d)

Model (26a) is a normal linear regression model, model (26b) a binary logistic
regression model, model (26c) is a 4-nomial model with log-odds function
$\vartheta_{k}-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w$ with
threshold parameters $\vartheta_{k}=\text{logit}(k/4)$ for $k=1,2,3$, and
model (26d) is a Weibull model with log-cumulative hazard function
$2\log(y)-\mu(\text{\boldmath$x$})-\tau(\text{\boldmath$x$})w$. We added
$50\,\%$ random right-censoring to the Weibull-generated data and also applied
a Cox proportional hazards model in addition to the Weibull model.

For the additive predictor
$\mu(\text{\boldmath$x$})+\tau(\text{\boldmath$x$})w$ we considered the 16
scenarios as specified in Table S. 4. Compared to Part A of this table, in
Part B half of the (negative) predictive effect is added to the prognostic
effect. We term the implied scenario where at least one variable exists which
is both prognostic (impact in $\mu(\text{\boldmath$x$})$) and predictive
(impact in $\tau(\text{\boldmath$x$})$) as overlay. $W(x_{1})$, $W(x_{3})$ and
$W(x_{4})$ depict that $W$ was drawn from a Bernoulli distribution with
$\pi(x_{1})$, $\pi(x_{3})$ or $\pi(x_{4})$, respectively.

Table S. 4: Experimental setup B. Confounding is present for non-constant propensities $\pi(\text{\boldmath$x$})$, an instrumental variable impacts $\pi(\text{\boldmath$x$})$ exclusively, heterogeneity of the treatment effect $\tau(\text{\boldmath$x$})$ is present when $\tau$ is non-constant, and overlay refers to variables being prognostic (impact in $\mu(\text{\boldmath$x$})$) and predictive (impact in $\tau(\text{\boldmath$x$})$) at the same time.  | Additive Predictor | Confounding | Instrument | Heterogeneity | Overlay
---|---|---|---|---|---
| $\mu(x_{3})+0\cdot W(x_{3})$ | yes | no | no | no
Part A | $\tau(x_{1},x_{2})W$ | no | no | yes | no
$\mu(x_{1})+\tau(x_{1},x_{2})W(x_{1})$ | yes | no | yes | yes
$\mu(x_{1})+\tau(x_{1},x_{2})W$ | no | no | yes | yes
| $\mu(x_{3})+\tau(x_{1},x_{2})W$ | no | no | yes | no
| $\mu(x_{3})+\tau(x_{1},x_{2})W(x_{3})$ | yes | no | yes | no
| $\tau(x_{1},x_{2})W(x_{3})$ | no | yes | yes | no
| $\mu(x_{3})+\tau(x_{1},x_{2})W(x_{4})$ | no | yes | yes | no
Part B | $\mu(x_{3})+0\cdot(W(x_{3})-0.5)$ | yes | no | no | no
$\tau(x_{1},x_{2})(W-0.5)$ | no | no | yes | yes
$\mu(x_{1})+\tau(x_{1},x_{2})(W(x_{1})-0.5)$ | yes | no | yes | yes
| $\mu(x_{1})+\tau(x_{1},x_{2})(W-0.5)$ | no | no | yes | yes
| $\mu(x_{3})+\tau(x_{1},x_{2})(W-0.5)$ | no | no | yes | yes
| $\mu(x_{3})+\tau(x_{1},x_{2})(W(x_{3})-0.5)$ | yes | no | yes | yes
| $\tau(x_{1},x_{2})(W(x_{3})-0.5)$ | no | yes | yes | yes
| $\mu(x_{3})+\tau(x_{1},x_{2})(W(x_{4})-0.5)$ | no | yes | yes | yes

In Part A of Table S. 4, the prognostic term and the predictive term are
separate and there is only overlay of prognostic and predictive effects when
both terms depend on $x_{1}$, i.e. $x_{1}$ is both prognostic and predictive
in this scenario. The treatment assignment probability may depend on $x_{1}$,
$x_{3}$, or $x_{4}$. In the third scenario, $x_{1}$ is a predictive confounder
(with impact on $\mu$, $\tau$, and $\pi$) and in the last two scenarios,
$x_{3}$ and $x_{4}$ can be understood as instruments with direct impact on
treatment assignment but without direct impact on the response. In Part B of
this table, half of the predictive effect is added to the prognostic effect,
so there is always overlay of both types of effects.

Again, we used random forests to estimate $\pi(\text{\boldmath$x$})$ and
gradient boosting machines to estimate $\eta_{0}(\text{\boldmath$x$})$ and
$\eta_{1}(\text{\boldmath$x$})$ as described in Section 4. We also applied the
same performance assessment (mean squared error evaluated on 1000 test
samples). The results are presented in Figures S. 2 and S. 3. The results for
the statistical analysis of RQ 1 to RQ 3 based on a normal linear mixed model
are presented in Table S. 5 to S. 7.

####  B.0.1 Results

Figure S. 2: Model-based forest results for Part A (Table S. 4), Cox means a
Cox model applied to the Weibull data. For the Weibull and Cox model,
treatment effects $\tau(\text{\boldmath$x$})$ are estimated as conditional log
hazard ratios. Direct comparison of model-based forests without centering
(Naive), model-based forests with local centering according to Robinson (1988)
of $Y$ and $W$ (Robinson) or only of $W$ (Robinson${}_{\widehat{W}}$).  Figure
S. 3: Model-based forest results for Part B (Table S. 4), Cox means a Cox
model applied to the Weibull data. For the Weibull and Cox model, treatment
effects $\tau(\text{\boldmath$x$})$ are estimated as conditional log hazard
ratios. Direct comparison of model-based forests without centering (Naive),
model-based forests with local centering according to Robinson (1988) of $Y$
and $W$ (Robinson) or only of $W$ (Robinson${}_{\widehat{W}}$).

Table S. 5: Results of RQ 1 for the experimental setups in Section B. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of the naive model-based forests, cells printed in italics indicate an inferior reference of naive model-based forests.  |  |  |  | Mean squared error ratio for RQ 1: Robinson vs. Naive
---|---|---|---|---
Part DGP | N | P |  | Normal | Binomial | Multinomial | Weibull | Cox
A $\mu(x_{3})+0\cdot W(x_{3})$ | 800 | 10 |  | 0.837 (0.762, 0.919) | 1.025 (0.713, 1.474) | 0.766 (0.660, 0.888) | 0.927 (0.769, 1.118) | 0.975 (0.759, 1.254)
|  | 20 |  | 0.808 (0.723, 0.903) | 1.123 (0.726, 1.735) | 0.761 (0.624, 0.928) | 0.931 (0.755, 1.149) | 0.975 (0.738, 1.287)
| 1600 | 10 |  | 0.795 (0.707, 0.894) | 0.912 (0.570, 1.462) | 0.742 (0.601, 0.915) | 0.821 (0.641, 1.053) | 0.848 (0.609, 1.180)
|  | 20 |  | 0.697 (0.601, 0.809) | 0.942 (0.526, 1.686) | 0.696 (0.533, 0.907) | 0.829 (0.613, 1.123) | 0.852 (0.569, 1.276)
$\tau(x_{1},x_{2})W$ | 800 | 10 |  | 1.106 (1.036, 1.180) | 3.295 (2.924, 3.714) | 0.976 (0.918, 1.037) | 1.222 (1.111, 1.344) | 1.220 (1.098, 1.356)
|  | 20 |  | 1.094 (1.023, 1.171) | 1.210 (1.044, 1.402) | 1.047 (0.986, 1.112) | 1.127 (1.019, 1.248) | 1.379 (1.234, 1.541)
| 1600 | 10 |  | 1.085 (1.005, 1.171) | 1.012 (0.878, 1.167) | 0.897 (0.827, 0.972) | 1.080 (0.938, 1.245) | 1.093 (0.946, 1.262)
|  | 20 |  | 1.128 (1.034, 1.231) | 1.021 (0.826, 1.261) | 1.004 (0.918, 1.098) | 1.108 (0.954, 1.287) | 1.140 (0.969, 1.342)
$\mu(x_{1})+\tau(x_{1},x_{2})W(x_{1})$ | 800 | 10 |  | 1.118 (1.054, 1.185) | 0.805 (0.706, 0.917) | 0.924 (0.885, 0.966) | 0.970 (0.895, 1.052) | 0.916 (0.833, 1.008)
|  | 20 |  | 0.988 (0.929, 1.051) | 0.943 (0.831, 1.070) | 1.032 (0.990, 1.074) | 0.923 (0.853, 0.998) | 0.912 (0.831, 1.002)
| 1600 | 10 |  | 1.078 (1.005, 1.156) | 0.188 (0.157, 0.225) | 0.872 (0.814, 0.935) | 0.875 (0.787, 0.973) | 0.805 (0.712, 0.910)
|  | 20 |  | 1.049 (0.970, 1.133) | 0.633 (0.532, 0.753) | 0.981 (0.919, 1.048) | 0.845 (0.747, 0.956) | 0.857 (0.742, 0.990)
$\mu(x_{1})+\tau(x_{1},x_{2})W$ | 800 | 10 |  | 1.076 (1.011, 1.146) | 0.465 (0.408, 0.531) | 0.964 (0.917, 1.013) | 1.013 (0.927, 1.106) | 1.059 (0.956, 1.173)
|  | 20 |  | 1.063 (1.000, 1.130) | 0.991 (0.862, 1.139) | 0.877 (0.831, 0.926) | 0.911 (0.832, 0.997) | 0.990 (0.889, 1.104)
| 1600 | 10 |  | 1.025 (0.947, 1.109) | 0.307 (0.259, 0.363) | 0.795 (0.738, 0.857) | 0.946 (0.834, 1.072) | 0.901 (0.788, 1.031)
|  | 20 |  | 1.027 (0.943, 1.119) | 0.918 (0.769, 1.096) | 0.880 (0.809, 0.958) | 0.933 (0.812, 1.072) | 0.863 (0.740, 1.007)
$\mu(x_{3})+\tau(x_{1},x_{2})W$ | 800 | 10 |  | 0.984 (0.931, 1.041) | 1.166 (1.032, 1.318) | 0.802 (0.760, 0.846) | 0.821 (0.753, 0.894) | 0.951 (0.860, 1.052)
|  | 20 |  | 0.875 (0.827, 0.926) | 0.897 (0.789, 1.020) | 0.800 (0.763, 0.839) | 0.690 (0.635, 0.751) | 0.724 (0.654, 0.801)
| 1600 | 10 |  | 0.951 (0.883, 1.023) | 0.946 (0.807, 1.110) | 0.776 (0.716, 0.840) | 0.820 (0.723, 0.930) | 0.880 (0.759, 1.019)
|  | 20 |  | 0.902 (0.834, 0.976) | 0.848 (0.693, 1.038) | 0.741 (0.685, 0.802) | 0.688 (0.598, 0.790) | 0.735 (0.622, 0.868)
$\mu(x_{3})+\tau(x_{1},x_{2})W(x_{3})$ | 800 | 10 |  | 0.865 (0.821, 0.911) | 1.161 (1.031, 1.308) | 0.778 (0.736, 0.823) | 0.820 (0.757, 0.888) | 0.844 (0.765, 0.931)
|  | 20 |  | 0.764 (0.726, 0.804) | 1.102 (0.974, 1.248) | 0.747 (0.711, 0.786) | 0.717 (0.663, 0.774) | 0.731 (0.666, 0.803)
| 1600 | 10 |  | 0.944 (0.882, 1.011) | 0.867 (0.744, 1.011) | 0.744 (0.690, 0.803) | 0.748 (0.666, 0.840) | 0.804 (0.698, 0.926)
|  | 20 |  | 0.918 (0.855, 0.986) | 0.982 (0.821, 1.174) | 0.697 (0.645, 0.754) | 0.632 (0.559, 0.715) | 0.701 (0.608, 0.809)
$\tau(x_{1},x_{2})W(x_{3})$ | 800 | 10 |  | 1.213 (1.142, 1.289) | 1.537 (1.372, 1.722) | 1.013 (0.954, 1.076) | 1.265 (1.154, 1.386) | 1.221 (1.096, 1.360)
|  | 20 |  | 1.148 (1.080, 1.220) | 1.349 (1.189, 1.530) | 1.179 (1.112, 1.250) | 1.335 (1.207, 1.477) | 1.468 (1.301, 1.657)
| 1600 | 10 |  | 1.227 (1.138, 1.323) | 1.178 (1.011, 1.372) | 0.916 (0.844, 0.993) | 1.150 (1.015, 1.304) | 0.965 (0.843, 1.105)
|  | 20 |  | 1.200 (1.106, 1.301) | 1.278 (1.052, 1.553) | 1.152 (1.051, 1.262) | 1.265 (1.094, 1.462) | 1.366 (1.167, 1.599)
$\mu(x_{3})+\tau(x_{1},x_{2})W(x_{4})$ | 800 | 10 |  | 1.017 (0.961, 1.076) | 1.191 (1.063, 1.334) | 0.877 (0.833, 0.924) | 0.913 (0.843, 0.990) | 0.942 (0.854, 1.040)
|  | 20 |  | 0.907 (0.857, 0.960) | 0.999 (0.880, 1.135) | 0.901 (0.860, 0.943) | 0.780 (0.720, 0.845) | 0.810 (0.736, 0.891)
| 1600 | 10 |  | 1.048 (0.978, 1.123) | 0.841 (0.722, 0.980) | 0.847 (0.788, 0.909) | 0.946 (0.841, 1.063) | 0.932 (0.811, 1.071)
|  | 20 |  | 0.967 (0.898, 1.042) | 0.943 (0.785, 1.133) | 0.847 (0.790, 0.910) | 0.752 (0.661, 0.857) | 0.782 (0.669, 0.913)
B $\mu(x_{3})+0\cdot(W(x_{3})-0.5)$ | 800 | 10 |  | 0.907 (0.820, 1.002) | 1.033 (0.925, 1.154) | 0.820 (0.723, 0.930) | 1.028 (0.867, 1.219) | 1.069 (0.850, 1.346)
|  | 20 |  | 0.818 (0.727, 0.921) | 1.047 (0.918, 1.194) | 0.739 (0.627, 0.872) | 0.980 (0.797, 1.205) | 1.025 (0.773, 1.359)
| 1600 | 10 |  | 0.783 (0.693, 0.885) | 0.896 (0.772, 1.041) | 0.746 (0.634, 0.879) | 0.827 (0.662, 1.032) | 0.869 (0.642, 1.177)
|  | 20 |  | 0.705 (0.597, 0.831) | 0.908 (0.755, 1.092) | 0.643 (0.503, 0.821) | 0.801 (0.610, 1.050) | 0.833 (0.577, 1.204)
$\tau(x_{1},x_{2})(W-0.5)$ | 800 | 10 |  | 1.013 (0.948, 1.082) | 1.112 (1.050, 1.177) | 0.973 (0.922, 1.026) | 1.222 (1.123, 1.330) | 1.414 (1.286, 1.555)
|  | 20 |  | 1.037 (0.971, 1.107) | 1.118 (1.047, 1.194) | 0.987 (0.928, 1.050) | 1.313 (1.189, 1.449) | 1.440 (1.291, 1.606)
| 1600 | 10 |  | 1.039 (0.956, 1.129) | 0.964 (0.890, 1.043) | 0.964 (0.896, 1.036) | 1.102 (0.974, 1.245) | 1.147 (1.004, 1.311)
|  | 20 |  | 1.033 (0.947, 1.127) | 0.971 (0.888, 1.062) | 0.983 (0.900, 1.074) | 1.140 (0.987, 1.317) | 1.240 (1.047, 1.468)
$\mu(x_{1})+\tau(x_{1},x_{2})(W(x_{1})-0.5)$ | 800 | 10 |  | 0.954 (0.901, 1.011) | 1.074 (1.016, 1.135) | 1.082 (1.028, 1.139) | 1.130 (1.048, 1.218) | 1.094 (0.999, 1.198)
|  | 20 |  | 0.931 (0.879, 0.986) | 1.131 (1.063, 1.203) | 1.140 (1.083, 1.201) | 1.032 (0.953, 1.118) | 1.069 (0.970, 1.179)
| 1600 | 10 |  | 1.017 (0.944, 1.096) | 0.994 (0.926, 1.068) | 0.987 (0.924, 1.054) | 0.947 (0.854, 1.051) | 0.944 (0.835, 1.067)
|  | 20 |  | 0.974 (0.900, 1.053) | 0.935 (0.858, 1.018) | 0.977 (0.904, 1.057) | 0.984 (0.870, 1.113) | 1.029 (0.892, 1.187)
$\mu(x_{1})+\tau(x_{1},x_{2})(W-0.5)$ | 800 | 10 |  | 1.024 (0.962, 1.089) | 0.969 (0.916, 1.024) | 0.874 (0.827, 0.923) | 1.152 (1.061, 1.250) | 1.179 (1.078, 1.289)
|  | 20 |  | 0.940 (0.884, 1.000) | 1.071 (1.004, 1.143) | 0.942 (0.889, 0.997) | 1.192 (1.091, 1.303) | 1.257 (1.138, 1.388)
| 1600 | 10 |  | 1.002 (0.925, 1.085) | 0.871 (0.809, 0.937) | 0.902 (0.838, 0.971) | 1.028 (0.917, 1.153) | 1.101 (0.963, 1.259)
|  | 20 |  | 0.976 (0.895, 1.064) | 0.871 (0.794, 0.956) | 0.867 (0.796, 0.944) | 1.001 (0.872, 1.149) | 1.039 (0.892, 1.211)
$\mu(x_{3})+\tau(x_{1},x_{2})(W-0.5)$ | 800 | 10 |  | 0.821 (0.780, 0.865) | 0.957 (0.905, 1.012) | 0.792 (0.752, 0.833) | 0.872 (0.805, 0.944) | 0.903 (0.820, 0.994)
|  | 20 |  | 0.705 (0.669, 0.744) | 0.845 (0.793, 0.900) | 0.728 (0.691, 0.766) | 0.730 (0.671, 0.793) | 0.775 (0.702, 0.856)
| 1600 | 10 |  | 0.850 (0.790, 0.915) | 0.840 (0.777, 0.909) | 0.787 (0.732, 0.847) | 0.863 (0.766, 0.974) | 0.935 (0.810, 1.081)
|  | 20 |  | 0.791 (0.734, 0.851) | 0.805 (0.737, 0.878) | 0.721 (0.668, 0.779) | 0.780 (0.682, 0.892) | 0.807 (0.687, 0.948)
$\mu(x_{3})+\tau(x_{1},x_{2})(W(x_{3})-0.5)$ | 800 | 10 |  | 0.810 (0.770, 0.853) | 1.084 (1.022, 1.149) | 0.927 (0.882, 0.973) | 0.979 (0.908, 1.056) | 0.927 (0.844, 1.017)
|  | 20 |  | 0.720 (0.684, 0.759) | 0.911 (0.857, 0.968) | 0.757 (0.721, 0.794) | 0.787 (0.733, 0.844) | 0.783 (0.719, 0.853)
| 1600 | 10 |  | 0.884 (0.825, 0.948) | 0.909 (0.844, 0.978) | 0.870 (0.814, 0.930) | 0.946 (0.855, 1.045) | 0.945 (0.831, 1.075)
|  | 20 |  | 0.829 (0.772, 0.890) | 0.847 (0.776, 0.925) | 0.755 (0.701, 0.813) | 0.830 (0.738, 0.934) | 0.832 (0.718, 0.965)
$\tau(x_{1},x_{2})(W(x_{3})-0.5)$ | 800 | 10 |  | 1.152 (1.085, 1.222) | 1.277 (1.206, 1.353) | 1.082 (1.023, 1.144) | 1.300 (1.200, 1.409) | 1.212 (1.101, 1.334)
|  | 20 |  | 1.107 (1.040, 1.179) | 1.228 (1.154, 1.306) | 1.127 (1.069, 1.190) | 1.470 (1.335, 1.618) | 1.498 (1.335, 1.680)
| 1600 | 10 |  | 1.132 (1.045, 1.226) | 1.136 (1.054, 1.225) | 1.040 (0.970, 1.116) | 1.228 (1.095, 1.377) | 1.123 (0.980, 1.286)
|  | 20 |  | 1.126 (1.032, 1.227) | 1.122 (1.026, 1.226) | 1.096 (1.006, 1.194) | 1.233 (1.077, 1.410) | 1.262 (1.078, 1.477)
$\mu(x_{3})+\tau(x_{1},x_{2})(W(x_{4})-0.5)$ | 800 | 10 |  | 0.864 (0.819, 0.911) | 1.021 (0.966, 1.080) | 0.823 (0.784, 0.864) | 0.907 (0.838, 0.982) | 0.891 (0.806, 0.986)
|  | 20 |  | 0.718 (0.682, 0.755) | 0.976 (0.921, 1.035) | 0.843 (0.806, 0.883) | 0.863 (0.805, 0.925) | 0.850 (0.780, 0.925)
| 1600 | 10 |  | 0.935 (0.869, 1.006) | 0.903 (0.841, 0.971) | 0.790 (0.737, 0.846) | 0.935 (0.837, 1.045) | 0.954 (0.833, 1.093)
|  | 20 |  | 0.861 (0.801, 0.925) | 0.897 (0.822, 0.978) | 0.762 (0.706, 0.821) | 0.899 (0.799, 1.011) | 0.911 (0.786, 1.055)

Table S. 6: Results of RQ 2 for the experimental setups in Section B. Comparison of mean squared errors for $\hat{\tau}(\text{\boldmath$x$})$ in the different scenarios. Estimates and simultaneous $95$ % confidence intervals were obtained from a normal linear mixed model with log-link. Cells printed in bold font correspond to a superior reference of the naive model-based forests, cells printed in italics indicate an inferior reference of naive model-based forests.  |  |  |  |  Mean squared error ratio for RQ 2: Robinson${}_{\hat{W}}$ vs. Naive
---|---|---|---|---
DGP | N | P |  | Normal | Binomial | Multinomial | Weibull | Cox
$\mu(x_{3})+0\cdot W(x_{3})$ | 800 | 10 |  | 1.182 (1.076, 1.298) | 1.012 (0.709, 1.446) | 1.284 (1.106, 1.491) | 1.109 (0.922, 1.335) | 1.063 (0.831, 1.361)
|  | 20 |  | 1.186 (1.059, 1.328) | 0.917 (0.598, 1.408) | 1.304 (1.069, 1.591) | 1.018 (0.821, 1.263) | 0.987 (0.744, 1.310)
| 1600 | 10 |  | 1.277 (1.137, 1.436) | 1.131 (0.711, 1.800) | 1.375 (1.116, 1.693) | 1.238 (0.967, 1.584) | 1.224 (0.883, 1.695)
|  | 20 |  | 1.392 (1.197, 1.618) | 1.084 (0.608, 1.930) | 1.430 (1.096, 1.866) | 1.203 (0.889, 1.629) | 1.200 (0.804, 1.790)
$\tau(x_{1},x_{2})W$ | 800 | 10 |  | 0.913 (0.856, 0.974) | 0.537 (0.499, 0.577) | 0.998 (0.939, 1.061) | 0.785 (0.712, 0.866) | 0.787 (0.706, 0.876)
|  | 20 |  | 0.913 (0.853, 0.976) | 0.808 (0.695, 0.938) | 0.954 (0.898, 1.013) | 0.906 (0.820, 1.002) | 0.766 (0.688, 0.854)
| 1600 | 10 |  | 0.920 (0.853, 0.993) | 0.822 (0.703, 0.961) | 1.060 (0.976, 1.151) | 0.879 (0.760, 1.017) | 0.816 (0.700, 0.952)
|  | 20 |  | 0.888 (0.814, 0.969) | 0.938 (0.756, 1.164) | 0.989 (0.904, 1.081) | 0.877 (0.753, 1.022) | 0.826 (0.697, 0.978)
$\mu(x_{1})+\tau(x_{1},x_{2})W(x_{1})$ | 800 | 10 |  | 0.948 (0.896, 1.004) | 1.000 (0.867, 1.153) | 0.959 (0.916, 1.005) | 0.898 (0.824, 0.979) | 0.986 (0.892, 1.089)
|  | 20 |  | 1.038 (0.977, 1.103) | 0.941 (0.823, 1.076) | 0.817 (0.781, 0.854) | 0.918 (0.843, 0.999) | 0.941 (0.851, 1.040)
| 1600 | 10 |  | 0.951 (0.887, 1.018) | 1.006 (0.785, 1.289) | 1.026 (0.954, 1.103) | 0.911 (0.810, 1.025) | 0.972 (0.848, 1.114)
|  | 20 |  | 1.001 (0.928, 1.080) | 0.959 (0.782, 1.175) | 0.924 (0.862, 0.991) | 1.020 (0.894, 1.164) | 1.148 (0.993, 1.327)
$\mu(x_{1})+\tau(x_{1},x_{2})W$ | 800 | 10 |  | 0.938 (0.881, 0.998) | 1.826 (1.593, 2.093) | 1.107 (1.055, 1.161) | 0.964 (0.881, 1.054) | 0.932 (0.841, 1.033)
|  | 20 |  | 0.945 (0.889, 1.004) | 0.984 (0.854, 1.133) | 1.096 (1.037, 1.158) | 1.108 (1.013, 1.213) | 1.052 (0.946, 1.170)
| 1600 | 10 |  | 0.978 (0.904, 1.058) | 2.057 (1.722, 2.456) | 1.196 (1.108, 1.291) | 1.036 (0.912, 1.177) | 1.051 (0.916, 1.205)
|  | 20 |  | 0.975 (0.895, 1.063) | 1.018 (0.848, 1.221) | 1.134 (1.042, 1.234) | 1.058 (0.919, 1.217) | 1.115 (0.955, 1.303)
$\mu(x_{3})+\tau(x_{1},x_{2})W$ | 800 | 10 |  | 1.018 (0.963, 1.076) | 0.823 (0.727, 0.933) | 1.257 (1.191, 1.326) | 1.231 (1.130, 1.340) | 1.107 (1.003, 1.221)
|  | 20 |  | 1.147 (1.084, 1.213) | 1.064 (0.933, 1.214) | 1.267 (1.208, 1.329) | 1.473 (1.356, 1.601) | 1.435 (1.298, 1.586)
| 1600 | 10 |  | 1.049 (0.975, 1.130) | 1.014 (0.863, 1.192) | 1.276 (1.179, 1.382) | 1.217 (1.072, 1.380) | 1.172 (1.013, 1.355)
|  | 20 |  | 1.108 (1.024, 1.199) | 1.138 (0.927, 1.397) | 1.358 (1.256, 1.469) | 1.479 (1.288, 1.698) | 1.434 (1.217, 1.690)
$\mu(x_{3})+\tau(x_{1},x_{2})W(x_{3})$ | 800 | 10 |  | 0.971 (0.918, 1.027) | 0.794 (0.702, 0.898) | 1.061 (0.998, 1.127) | 0.894 (0.815, 0.980) | 0.938 (0.841, 1.045)
|  | 20 |  | 0.980 (0.925, 1.038) | 0.777 (0.679, 0.889) | 0.946 (0.893, 1.003) | 0.877 (0.798, 0.964) | 0.904 (0.809, 1.011)
| 1600 | 10 |  | 0.989 (0.922, 1.061) | 0.819 (0.686, 0.977) | 1.127 (1.039, 1.222) | 0.961 (0.841, 1.099) | 0.996 (0.852, 1.164)
|  | 20 |  | 0.972 (0.902, 1.049) | 0.843 (0.692, 1.026) | 1.002 (0.916, 1.096) | 1.021 (0.884, 1.181) | 1.006 (0.854, 1.186)
$\tau(x_{1},x_{2})W(x_{3})$ | 800 | 10 |  | 1.043 (0.988, 1.100) | 0.570 (0.502, 0.646) | 1.074 (1.014, 1.138) | 1.027 (0.948, 1.112) | 1.055 (0.961, 1.159)
|  | 20 |  | 1.103 (1.045, 1.165) | 0.721 (0.634, 0.820) | 0.982 (0.930, 1.036) | 1.155 (1.065, 1.251) | 1.116 (1.018, 1.223)
| 1600 | 10 |  | 1.014 (0.948, 1.084) | 0.774 (0.660, 0.910) | 1.140 (1.053, 1.235) | 1.072 (0.958, 1.200) | 1.056 (0.925, 1.206)
|  | 20 |  | 1.021 (0.949, 1.098) | 0.883 (0.736, 1.058) | 1.151 (1.063, 1.246) | 1.158 (1.028, 1.305) | 1.028 (0.901, 1.173)
$\mu(x_{3})+\tau(x_{1},x_{2})W(x_{4})$ | 800 | 10 |  | 1.490 (1.420, 1.563) | 0.814 (0.725, 0.914) | 1.295 (1.233, 1.360) | 1.476 (1.374, 1.586) | 1.597 (1.466, 1.738)
|  | 20 |  | 1.706 (1.625, 1.792) | 1.019 (0.898, 1.156) | 1.272 (1.219, 1.329) | 1.701 (1.581, 1.831) | 1.810 (1.661, 1.972)
| 1600 | 10 |  | 1.268 (1.193, 1.348) | 1.026 (0.873, 1.205) | 1.381 (1.292, 1.477) | 1.604 (1.451, 1.774) | 1.782 (1.585, 2.003)