# Individual Treatment Effect Prediction
for ALS Patients

Heidi Seibold
Universität Zürich
&Achim Zeileis
Universität Innsbruck
&Torsten Hothorn
Universität Zürich<EMAIL_ADDRESS>
Seibold, Zeileis, Hothorn Individual Treatment Effect Prediction for ALS
Patients Individual Treatment Effect Prediction for ALS Patients A treatment
for a complicated disease may be helpful for some but not all patients, which
makes predicting the treatment effect for new patients important yet
challenging. Here we develop a method for predicting the treatment effect
based on patient characteristics and use it for predicting the effect of the
only drug (Riluzole) approved for treating Amyotrophic Lateral Sclerosis
(ALS). Our proposed method of model-based random forests detects similarities
in the treatment effect among patients and on this basis computes personalised
models for new patients. The entire procedure focuses on a base model, which
usually contains the treatment indicator as a single covariate and takes the
survival time or a health or treatment success measurement as primary outcome.
This base model is used both to grow the model-based trees within the forest,
in which the patient characteristics that interact with the treatment are
split variables, and to compute the personalised models, in which the
similarity measurements enter as weights. We applied the personalised models
using data from several clinical trials for ALS from the PRO-ACT database. Our
results indicate that some ALS patients benefit more from the drug Riluzole
than others. Our method allows shifting from stratified medicine to
personalised medicine and can also be used in assessing the treatment effect
for other diseases studied in a clinical trial. Personalised medicine, random
forest, treatment effect, model-based recursive partitioning Personalised
medicine, random forest, treatment effect, model-based recursive partitioning
Heidi Seibold
Institut für Epidemiologie, Biostatistik und Prävention, Abteilung
Biostatistik
Universität Zürich
Hirschengraben 84, CH-8001 Zürich, Switzerland
Achim Zeileis
Department of Statistics
Faculty of Economics and Statistics
University of Innsbruck
Universitätsstr. 15
A-6020 Innsbruck, Austria
Torsten Hothorn
Institut für Epidemiologie, Biostatistik und Prävention, Abteilung
Biostatistik
Universität Zürich
Hirschengraben 84, CH-8001 Zürich, Switzerland
E-mail:
URL: http://user.math.uzh.ch/hothorn/

##  1 Introduction

Amyotrophic lateral sclerosis (ALS) is a deadly disease that affects motor
neurons in the brain and spinal cord, i.e. the neurons responsible for
voluntary muscle control. Riluzole (Rilutek) is the only approved drug for
this disease to date. According to the European Medicines Agency (2012),
Riluzole prolongs the median survival of ALS patients, depending on the dose,
by a few months. Several side effects, such as sickness, weakness or increased
liver enzyme levels are mentioned (European Medicines Agency, 2012). Knowledge
how Riluzole works on the nervous system of ALS patients is limited. The PRO-
ACT database (Atassi _et al._ , 2014) is the largest database containing
clinical trial data of ALS patients available and was initiated to retrieve
more information on the disease. It contains data from 17 ALS studies
conducted between 1990 and 2010\. Using these data, we aimed at finding out
more about the effect of Riluzole on the health and survival of patients.

Before statistical analysis and $p$-values entered into medical progress,
doctors treated patients individually based on their experiences and knowledge
(Weisberg, 2015). Since the beginning of the “golden age of randomised
clinical trials”, however, medication became more and more standardised.
Nowadays, much knowledge about the effect of drugs has accumulated,
cornerstone drugs such as antibiotics have been used for decades and many
diseases can be treated successfully; however, providing new drugs for the
general public is more difficult. Diseases such as ALS are too complex to
treat all patients in the same way. Therefore, there is a need to return to
more individualised treatments, but this time with the use of statistical
concepts.

In the past years, there has been an immense effort towards personalised
medicine in the analysis of randomised controlled trials. The goal is to
identify predictive factors, i.e. factors that interact with the treatment
(Italiano, 2011), such as biomarkers, other treatments, and environmental
circumstances. In the following, we will refer to these factors as patient
characteristics. Prognostic factors, i.e. factors that directly affect the
patient’s outcome, are only of secondary interest, but should not be
neglected, because they not only change the general level of the outcome –
showing in the individual intercept – but might also be predictive and
prognostic (Seibold _et al._ , 2015). For drugs for which the biological mode
of action is unknown, predictive and prognostic factors should first be
identified in a data-driven way. New hypotheses can then be generated and new
trials can be planned based on these hypotheses. In this first step we are
asking whether a certain patient characteristic is relevant and not why.

Many new statistical methods in the field of stratified medicine, i.e.
subgroup analysis, have been developed. Subgroup analyses aim at finding
groups of patients that have differential treatment effects. Most of the
methods are based on recursive partitioning (trees) and/or interaction models
(Ciampi _et al._ , 1995; Kehl and Ulm, 2006; Dusseldorp and Van Mechelen,
2013; Loh _et al._ , 2015; Tian _et al._ , 2014; Foster _et al._ , 2015; Zhang
_et al._ , 2012). The tree methods for subgroup analyses have specialised
splitting procedures for partitioning the patients into groups with higher and
lower treatment effect. Interaction models evaluate the interaction between
the treatment and given patient characteristics. The idea behind methods of
subgroup analyses in general is to obtain a treatment effect
$\beta(\mathbf{z})$ that depends on the patient characteristics $\mathbf{z}$.
For example, the treatment effect could depend on the age of patients, in
which patients less than 40 years of age improve through the treatment,
patients between 40 and 60 do not improve and patients older than 60 years
improve, but less than the patients under 40 years:

$\displaystyle\beta(\mathbf{z})$
$\displaystyle=\left\\{\begin{array}[]{rl}1&\text{if }z_{\text{age}}<40\\\
0&\text{if }40\leq z_{\text{age}}<60\\\ 0.5&\text{if }60\leq z_{\text{age}}\\\
\end{array}\right.$ (4)

However, the assumption that the treatment effect is a step function may be
too restrictive, and $\beta(\mathbf{z})$ in reality may be a smooth
interaction function. In other words, personalised medicine is required
instead of stratified medicine. Because methods for subgroup analyses again
generalise the treatment effect for a group of patients, it can only be
considered as a step in the direction toward personalised medicine. We provide
a method that can estimate smooth treatment effect functions using model-based
random forests and weighted models. More importantly, this method provides an
estimate for the treatment effect of a future patient, thereby allowing a
decision to be made whether treatment of this patient is appropriate.

##  2 Methods

Seibold _et al._ (2015) introduced a means of conducting subgroup analysis
using model-based recursive partitioning. One first defines a model
$\mathcal{M}((Y,\mathbf{X}),\boldsymbol{\vartheta})$ with primary endpoint
$Y$, covariates $\mathbf{X}$ including the treatment indicator

$\displaystyle X_{A}$ $\displaystyle=\left\\{\begin{array}[]{l
l}1&\quad\text{if patient received the (new) treatment}\\\ 0&\quad\text{if
patient received no treatment (or standard of care),}\end{array}\right.$ (7)

and parameter vector $\boldsymbol{\vartheta}$. In the following we will
consider likelihood models (e.g. generalised linear models or parametric
survival models) where the model parameters $\boldsymbol{\vartheta}$ can be
estimated by maximising the log-likelihood
$l((Y,\mathbf{X}),\boldsymbol{\vartheta})$ of those models (e.g. Gaussian log-
likelihood or Weibull log-likelihood) or equivalently by solving the score
equation

$\displaystyle\sum\limits_{i=1}^{N}s((y,\mathbf{x})_{i},{\boldsymbol{\vartheta}})=\mathbf{0}$
(8)

with

$\displaystyle s((y,\mathbf{x})_{i},{\boldsymbol{\vartheta}})=\frac{\partial
l((y,\mathbf{x})_{i},{\boldsymbol{\vartheta}})}{\partial{\boldsymbol{\vartheta}}}.$
(9)

In most applications the model contains only an intercept $\alpha$ and a
treatment effect $\beta$, i.e. $\mathbf{X}=(1,X_{A})$ and
$\boldsymbol{\vartheta}=(\alpha,\beta)^{\top}$ but more parameters are
possible, such as coefficients of additional regressors or scale and shape
parameters for the response distribution. Technically, there can also be more
than two treatment groups. For simplicity, we will focus on the simple case
with intercept and treatment effect and two treatment groups. The method
obtains subgroups $\\{\mathcal{B}_{b=1,\dots,B}\\}$ that differ with regard to
the treatment effect $\beta$ and potentially the intercept $\alpha$. The
subgroups are defined by patient characteristics
$\mathbf{Z}=(Z_{1},\dots,Z_{J})\in\mathcal{Z}$. Hence the intercept and
treatment parameters can be written as a function of the patient
characteristics $\mathbf{z}$.

Conceptually, the partitioned model parameters $\alpha(\mathbf{z})$ and
$\beta(\mathbf{z})$ might depend on $\mathbf{z}$ in a more complex way than a
simple tree structure. Therefore, the model parameters are not step functions,
but rather smooth interaction functions, so that an individual treatment
effect (as in personalised medicine) can be computed for each patient instead
of only for each subgroup of patients (as in stratified medicine). The
function $\beta(\mathbf{z})$ can then be understood as an estimate of the
individual treatment effect of a patient with patient characteristics
$\mathbf{z}$.

The most intuitive step from a tree structure to a more complex structure is
to use a random forest instead of a single tree. Hence we propose a model-
based random forest for personalised medicine, which can be used to predict
the treatment effect of future patients using personalised models.

###  2.1 Random forest

Random forests (Breiman, 2001) compute an ensemble of $T$ trees. The proposed
algorithm draws subsamples $\mathcal{L}_{t},~{}t=1,\dots,T$ of the given $N$
observations and fits a model-based tree to each subsample using a randomly
sampled set of candidate split variables $\mathbf{z}$. The data
$\mathcal{L}_{t}^{c}$ that were not in the learning sample for tree $t$ are
called out-of-bag data.

#### Split procedure

The special feature of our method is the split procedure, which is based on
the empirical estimating function

$\displaystyle\boldsymbol{s}$
$\displaystyle=\begin{pmatrix}s_{\hat{\alpha}}((y,\mathbf{x})_{1},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{1},\hat{\boldsymbol{\vartheta}})\\\
s_{\hat{\alpha}}((y,\mathbf{x})_{2},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{2},\hat{\boldsymbol{\vartheta}})\\\
\vdots&\vdots\\\
s_{\hat{\alpha}}((y,\mathbf{x})_{N},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{N},\hat{\boldsymbol{\vartheta}})\\\
\end{pmatrix}$ (10)

which contains the score contributions
$s_{\hat{\alpha}}((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}})$ and
$s_{\hat{\beta}}((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}})$. The score
contributions are the partial derivatives of the log-likelihood with respect
to $\alpha$ or $\beta$ respectively evaluated at the $N$ observed data points
and the estimated parameters
$\hat{\boldsymbol{\vartheta}}=(\hat{\alpha},\hat{\beta})^{\top}$(Zeileis _et
al._ , 2008). The matrix of score contributions contains information on the
deviation from the model fit for all parameters and observations of a given
model $\mathcal{M}((Y,\mathbf{X}),\boldsymbol{\vartheta})$. The contributions
can thus be seen as residuals.

To obtain a split in model-based recursive partitioning for this setup the
following steps have to be performed:

  * •

Compute the prespecified (parametric) model
$\mathcal{M}((Y,\mathbf{X}),\boldsymbol{\vartheta})$.

  * •

Compute the associated score matrix $\boldsymbol{s}$.

  * •

Perform tests of independence between the residuals (partial score vectors)
and the partitioning variables:

$\displaystyle H_{0}^{\alpha,j}$
$\displaystyle:\quad{s}_{\hat{\alpha}}((Y,\mathbf{X}),\hat{\boldsymbol{\vartheta}})\quad\bot\quad
Z_{j}$ $\displaystyle H_{0}^{\beta,j}$
$\displaystyle:\quad{s}_{\hat{\beta}}((Y,\mathbf{X}),\hat{\boldsymbol{\vartheta}})\quad\bot\quad
Z_{j}\qquad j=1,\dots,J$

The smallest $p$-value corresponds to the greatest deviation from the model
assumption, that intercept and treatment parameter are the same for all
patients in the given node/subgroup.

  * •

If any $p$-value is lower than the significance level, select the partitioning
variable that has the highest association (lowest $p$-value) to any of the
relevant residuals for the split.

  * •

Search for the optimal split point in the selected partitioning variable using
a suitable criterion, such that the models in the resulting daughter-nodes
have as little association between the partitioning variable and the residuals
as possible.

This split procedure is repeated until a stop criterion is met. This can be,
for example, when no $p$-values are lower than the significance level or if
subgroups become too small. In the end a tree is obtained with disjoint
subgroups

$\displaystyle{\bigcup\limits_{b}^{\bullet}}\mathcal{B}_{b}=\mathcal{Z}$ (11)

or a random forest with $T$ trees and disjoint subgroups for each tree

$\displaystyle\bigcup\limits_{b}^{\bullet}\mathcal{B}_{tb}=\mathcal{Z}\quad\forall
t=1,\dots,T.$ (12)

The independence tests can be performed using permutation tests (Hothorn _et
al._ , 2006a, b) or, for reasonably large samples, using M-fluctuation tests
(Zeileis and Hornik, 2007; Zeileis _et al._ , 2008). Note that multiplicity
adjustment, e.g. Bonferroni correction, is recommended. More details on the
algorithm and the test procedures used are documented in the appendix.

###  2.2 Personalised models

In personalised medicine, the goal is to learn how much a person will profit
from a given treatment and what would happen if the standard of care or no
treatment is given. For any patient, it is possible to compute a personalised
model based on the similarity of this observation to the observations in the
training data. The similarity $w_{i}(\mathbf{z})$ of a patient $i$ to any
patient from the training set is given as the number of times the patients are
classified in the same subgroup by the single trees in the random forest

$\displaystyle w_{i}(\mathbf{z})$
$\displaystyle=\sum\limits_{t=1}^{T}\sum\limits_{b=1}^{B_{t}}(\mathbf{z}_{i}\in\mathcal{B}_{tb})\wedge(\mathbf{z}\in\mathcal{B}_{tb}),$
(13)

with $T$ being the number of trees used for the computation of the forest and
$B_{t}$ being the number of subgroups from tree $t$ (Hothorn _et al._ , 2004;
Meinshausen, 2006). If patient $i$ is part of the training set, the weights
can be computed out-of-bag, i.e. the only trees ($t=1,\dots,T$) considered are
those where patient $i$ is not in the subset $\mathcal{L}_{t}$ for the
computation.

To obtain the personalised model for patient $i$ the base model is recomputed
with the weighted training data or equivalently by using the weighted score
function in the score equation

$\displaystyle\sum\limits_{k=1}^{N}w_{i}(\mathbf{z}_{k})\cdot
s((y,\mathbf{x})_{i},{\boldsymbol{\vartheta}}(\mathbf{z}_{i}))=0.$ (14)

In other words every patient $k$ from the training set is included
$w_{i}(\mathbf{z}_{k})$ times in the “new data set” to compute the
personalised model for patient $i$. In the following the parameters estimated
from this model will be denoted by
$\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})=(\hat{\alpha}(\mathbf{z}_{i}),\hat{\beta}(\mathbf{z}_{i}),\dots)$.

Using the personalised models it is possible to obtain a log-likelihood. From
the personalised model for patient $i$, the log-likelihood contribution
$l((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i}))$ for this
observation is computed. The log-likelihood then is

$\displaystyle\sum\limits_{i=1}^{N}l\left((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right),$
(15)

which we refer to as forest log-likelihood.

###  2.3 Improvement through personalised models

To check whether the personalised models actually lead to an improvement of
the base model we test the hypothesis

$\displaystyle H_{0}:$
$\displaystyle\underbrace{\alpha(\mathbf{Z})\equiv\alpha}_{H_{0}^{\alpha}}$
(16) $\displaystyle\quad\bigcap$
$\displaystyle\underbrace{\beta(\mathbf{Z})\equiv\beta}_{H_{0}^{\beta}}.$ (17)

This strict null hypothesis is to be rejected, if any of the patient
characteristics contain information on the outcome or the treatment effect. To
conduct the test, one can proceed as follows:

  * •

Compute the forest log-likelihood and the log-likelihood of the base model and
calculate their difference.

  * •

Draw parametric bootstrap samples from the base model.

  * •

Compute the forest log-likelihood and the log-likelihood of the base model in
the bootstrap samples and again compute the differences.

  * •

The $p$-value is then the proportion of bootstrap samples in which the
difference in log-likelihoods exceeds the observed difference in the original
data. Note, that this $p$-value will be very low or even zero when the patient
characteristics contain information on the outcome or the treatment effect.

In practice, one may be interested in just $H_{0}^{\beta}$, but testing the
sub-hypotheses $H_{0}^{\alpha}$ and $H_{0}^{\beta}$ separately is not straight
forward. An approximation would be to compute the personalised models using a
forest that splits only based on the partial score function with respect to
$\alpha$ or $\beta$ respectively. Patient characteristics, however, are often
not exclusively predictive or prognostic but may be both. Also, if a patient
characteristic is purely prognostic, this still may result in a pattern in
both partial score functions. For more details, see Seibold _et al._ (2015).

###  2.4 Dependence plots

A partial dependence plot describes the dependence of a function (in our case
the treatment effect $\hat{\beta}(\mathbf{z})$) and a variable (in our case a
partitioning variable) (Hastie _et al._ , 2009). The partial dependence plot
resulting from a model-based tree would show a step function. The partial
dependence from a random forest can be smoother for continuous partitioning
variables. It can be obtained by plotting $\hat{\beta}({z}_{j})$ against
$z_{j}$ for each partitioning variable $j=1,\dots,J$.

###  2.5 Variable importance

The variable importance for the random forest is computed based on the tree
log-likelihoods. For a given forest computed with $T$ trees the log-likelihood
is computed as follows:

  * •

Select the out-of-bag data $\mathcal{L}_{t}^{c}$ and determine the terminal
node/subgroup to which each observation $i$ belongs to.

  * •

Compute the log-likelihood contribution of each observation
$i\in\mathcal{L}_{t}^{c}$ based on the respective model in the terminal
node/subgroup with parameters $\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})$.

  * •

Compute the out-of-bag log-likelihood as the sum of the contributions

$\displaystyle l_{t}$
$\displaystyle=\sum\limits_{i\in\mathcal{L}_{t}^{c}}l((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})).$
(18)

To obtain the variable importance of a given variable $z_{j},~{}j=1,\dots,J$,
the variable is permuted. The log-likelihood is computed as above, except that
the column with information about $z_{j}$ in the out-of-bag data is replaced
by the permuted $z_{j}$. We denote the log-likelihood of tree $t$ with
variable $z_{j}$ permuted by $l_{t}^{(j)}$. The variable importance is then

$\displaystyle\text{VI}_{j}=\frac{1}{T}\sum\limits_{t=1}^{T}\left[l_{t}-l_{t}^{(j)}\right].$
(19)

If the variable importance is high the variable is an important predictive
and/or prognostic factor. Note that due to the absolute differences the
variable importances may be negative.

##  3 Results

###  3.1 PRO-ACT data

The Pooled Open Access Clinical Trials (PRO-ACT,
https://nctu.partners.org/ProACT) database contains longitudinal data of ALS
patients that participated in one of 16 phase II and III trials and one
observational study. It is a project initiated by the non-profit organisation
Prize4Life (http://www.prize4life.org/) to enhance knowledge about ALS. It
contains information on a broad variety of patient characteristics, such as
vital signs, the patient’s and family’s history, and treatment information.
Identification criteria, such as study centres, are not included in the
database. Also collected are the survival time and the ALS functional rating
scale (ALSFRS), which is a score measuring the patients ability of living a
normal life (Brooks _et al._ , 1996). The ALSFRS is a sum-score of ten items,
each of which ranges between zero and four, where zero represents complete
inability and four represents normal ability. The items are speech,
salivation, swallowing, hand-writing, cutting food and handling utensils,
dressing and hygiene, turning in bed and adjusting bed clothes, walking,
climbing stairs, and breathing. As outcomes in the study we used both the
survival time (denoted by survival) and the ALSFRS six months after treatment
start (denoted by $\text{ALSFRS}_{6}$) and identified patient characteristics
that influence the effect of Riluzole on these outcomes. For the two outcome
variables, we obtained two different data sets. We only included observations
that contain information on the respective outcome variable and only patient
characteristics that have fewer than 50% missing values. The survival time
data set contains 3306 observations and 18 patient characteristics. The ALSFRS
data set contains 2534 observations and 57 patient characteristics.

Table 1: ALSFRS base model (Gaussian generalised linear model with log-link and offset). Given are the parameter estimates, their standard error and the Wald confidence interval. | Estimate | Std.Error | 2.5 % | 97.5 %
---|---|---|---|---
$\alpha$ | $-0.1595$ | $0.0065$ | $-0.1722$ | $-0.1468$
$\beta$ | $0.0091$ | $0.0077$ | $-0.0060$ | $0.0242$
Table 2: Survival time base model (Weibull model). Given are the parameter estimates, their standard error and the Wald confidence interval. | Estimate | Std.Error | 2.5 % | 97.5 %
---|---|---|---|---
$\alpha_{1}$ | $6.7070$ | $0.0323$ | $6.6437$ | $6.7703$
$\beta$ | $0.1073$ | $0.0387$ | $0.0314$ | $0.1832$
$\log(\alpha_{2})$ | $-0.5833$ | $0.0271$ | $-0.6364$ | $-0.5302$

Tables 1 and 2 show the estimates including standard errors obtained from the
base model for each outcome. For the ALSFRS this base model is given by

$\displaystyle\mathbb{E}\left.\left(\frac{\text{ALSFRS}_{6}}{\text{ALSFRS}_{0}}\right|X=x\right)=\frac{\mathbb{E}(\text{ALSFRS}_{6}|X=x)}{\text{ALSFRS}_{0}}=\exp\\{\alpha+\beta
x_{A}\\},$ (20)

which represents a Gaussian generalised linear model with log-link and offset
$\log(\text{ALSFRS}_{0})$ where $\text{ALSFRS}_{0}$ is the ALSFRS that was
measured at the time of treatment start. The base model for the survival time
is given by the Weibull model

$\displaystyle\mathbb{P}(T\leq\text{survival}|X=x)=F\left(\frac{\log(\text{survival})-\alpha_{1}-\beta
x_{A}}{\alpha_{2}}\right),$ (21)

where $F$ is the cumulative distribution function of the Gompertz
distribution. Note that the Weibull model has a scale parameter in addition to
the intercept, so that both $\alpha_{1}$ and $\alpha_{2}$ control the
appearance of the baseline hazard. In the notation of Equation 4, this leads
to $\mathbf{\vartheta}=(\alpha_{1},\alpha_{2},\beta)^{\top}$.

###  3.2 Personalised models

Figure 1: Kernel density estimates of the personalised parameter estimates for
the ALSFRS. Figure 2: Distribution of the personalised parameter estimates for
the survival time. The baseline hazard functions are given in the left panel;
the kernel density estimate of the treatment effect estimate are given in the
right panel.

We computed personalised models for all observations in the respective
training data, which were used to obtain the random forest. The distribution
of parameter estimates in the personalised models is given in Figure 1 for the
ALSFRS and in Figure 2 for the survival time. Figure 1 shows that all patients
are predicted to have a positive Riluzole effect, i.e. for all patients taking
Riluzole, a higher ALSFRS is achieved compared to those not taking Riluzole.
However, there is a variability in the treatment effects, and the distribution
of the treatment effect is bimodal (as is the distribution of the intercept).
The treatment effect estimated from the base model is between the two modes.
The lowest treatment effect a person in this data set is predicted to have is
0.0027.

For the survival time, the lowest predicted treatment effect is 0.0717.
However, the value of the treatment effect in the personalised survival models
cannot be interpreted in isolation; its meaning depends on the shape of the
baseline hazard, i.e. on $\alpha_{1}$ and $\alpha_{2}$. Instead of depicting
the densities of the two baseline hazard parameters, in Figure 2 we show the
baseline hazard curves. The baseline hazard varies for different patients and
there is a gap in the middle. The baseline hazard estimated from the base
model lies close to that gap.

From the personalised models, we obtained the “forest log-likelihoods” for
both outcomes. For the Gaussian GLM with log-link and offset, the log-
likelihood contribution for observation $i$ is defined as

$\displaystyle l$
$\displaystyle\left((\text{ALSFRS}_{6},\text{ALSFRS}_{0},\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right)$
$\displaystyle\quad=\left(\text{ALSFRS}_{6i}-\exp\left(\mathbf{x}_{i}^{\top}\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right)\cdot\text{ALSFRS}_{0i}\right)^{2}$
(22)

with $\mathbf{x}_{i}=(1,x_{Ai})^{\top}$ and
$\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})=(\hat{\alpha}(\mathbf{z}_{i}),\hat{\beta}(\mathbf{z}_{i}))^{\top}$.
For the Weibull model the log-likelihood contribution for observation $i$ is

$\displaystyle l$
$\displaystyle\left((\text{survival},\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right)$
$\displaystyle\quad=\delta_{i}\log(\hat{\alpha}_{2}(\mathbf{z}_{i}))-\delta_{i}\frac{\text{survival}_{i}-\mathbf{x}_{i}^{\top}\hat{\boldsymbol{\vartheta}}^{*}(\mathbf{z}_{i})}{\hat{\alpha}_{2}(\mathbf{z}_{i})}+\exp\left(\frac{\text{survival}_{i}-\mathbf{x}_{i}^{\top}\hat{\boldsymbol{\vartheta}}^{*}(\mathbf{z}_{i})}{\hat{\alpha}_{2}(\mathbf{z}_{i})}\right)$
(23)

with $\mathbf{x}_{i}=(1,x_{iA})^{\top}$,
$\hat{\boldsymbol{\vartheta}}^{*}(\mathbf{z}_{i})=(\hat{\alpha}_{1}(\mathbf{z}_{i}),\hat{\beta}(\mathbf{z}_{i}))^{\top}$
and $\delta_{i}$ as the censoring indicator.

Figure 3: Difference in log-likelihoods between forest and base model on the
original data (dashed lines; $H_{1}$, the usual forest; $H_{1}^{\alpha}$ the
forest that splits based on $\alpha$; $H_{1}^{\beta}$ the forest that splits
based on $\beta$) and on 50 samples simulated from the base model (density
curve) for the ALSFRS outcome. Figure 4: Difference in log-likelihoods between
forest and base model on the original data (dashed lines; $H_{1}$, the usual
forest; $H_{1}^{\alpha}$ the forest that splits based on $\alpha$;
$H_{1}^{\beta}$ the forest that splits based on $\beta$) and on 50 samples
simulated from the base model (density curve) for the survival outcome.

As can be seen in Figures 3 and 4 the forest log-likelihoods are higher than
the log-likelihoods of the base models for both the ALSFRS and the survival
time. The figures show the difference in log-likelihood between the forest and
the corresponding base model. To show that this difference is not due to
overfitting, we drew 50 samples from the base models, i.e. 50 parametric
bootstrap samples for which the assumption holds that the intercept (or
baseline hazard) and treatment effect are the same for all patients. ALSFRS
values are drawn from a normal distribution truncated at zero to assure
positivity. (The effect of truncation is virtually negligible; only two
observations had a truncation probability of more than 1%.) The survival times
are drawn from a Weibull distribution censored at the originally observed
censoring times (if exceeded). The difference in log-likelihoods for both
ALSFRS and survival time are distributed close to zero, with a slight shift to
the right, for the parametric bootstrap samples. The large difference in the
ALS data supports the assumption that the base models are not ideal and
personalised models are meaningful (the respective p-values are both zero). To
approximately check the sub-hypotheses given in equation 16 and 17, we also
computed log-likelihoods of the two forests that split only with respect to
one of the partial score functions – either intercept (or baseline hazard) or
treatment effect. For the ALSFRS, both the forest under $H_{1}^{\alpha}$
(computed with only splitting based on the partial score function with respect
to the intercept $\alpha$), as well as the forest under $H_{1}^{\beta}$
(computed with only splitting based on partial score function with respect to
the treatment effect $\beta$), greatly improved compared to the base model.
The difference in log-likelihood between the forest under $H_{1}^{\alpha}$ and
the base model is even greater than between the original forest ($H_{1}$) and
the base model. For the survival time, the log-likelihoods of the original
forest and the forest under $H_{1}^{\alpha}$ (based on splits in the partial
score function with respect to the baseline hazard) are very close to each
other. Splitting only based on the partial score function with respect to the
treatment effect ($H_{1}^{\beta}$) already improves the log-likelihood but not
as much as splitting based on both intercept and treatment effect ($H_{1}$).
The good performance of the forests under $H_{1}^{\alpha}$ indicate that (1)
there are no predictive patient characteristics, (2) all predictive patient
characteristics are also prognostic, or (3) the predictive nature of the
predictive patient characteristics are so strong that it has enough impact on
the structure of the partial score function with respect to $\alpha$.

###  3.3 Dependence plots

(a) Dependence plot for the time in days between disease onset and treatment
start.

(b) Dependence plot for the creatinine level in mmol/L.

(c) Dependence plot for the phosphorus level in mmol/L.

(d) Dependence plot for the forced vital capacity (volume of air in litres
that can forcibly be blown out after full inspiration).

Figure 5: Dependence plots for the four patient characteristics with the
highest variable importance from the ALSFRS forest.

(a) Dependence plot for the age.

(b) Dependence plot for the time in days between disease onset and treatment
start.

(c) Dependence plot for the weakness indicator.

(d) Dependence plot for the height.

Figure 6: Dependence plots for the four patient characteristics with the
highest variable importance from the survival time forest.

The dependence plots as shown in Figures 5 and 6 can be obtained for any
partitioning variable. Here we show the dependence plots for the four
variables with the highest variable importance (see Section 3.4). For
continuous variables, such as age, we show a scatter plot, as before. For
categorical variables, such as the variable weakness, which indicates whether
a patient suffers from muscle weakness (yes/no), boxplots giving the variation
of $\beta(z)$ and a square representing $\bar{\beta}({z})$, i.e. the mean, are
a meaningful way of representing the dependence between treatment effect and
the given variable.

The most obvious pattern of the four graphs in Figure 5 is shown in Subfigure
5(d), in which the personalised treatment effects are plotted against the
forced vital capacity (FVC). Patients with a low lung function (low FVC) are
predicted to have a higher treatment effect than those with better lung
function. The graph shows a relatively clear cut at approximately three
litres. This indicates that FVC is a predictive factor. For the time between
disease onset and treatment start, the pattern is less clear. Patients with a
short as well as those with a long time between disease onset and treatment
start seem to benefit most. Also for the creatinine value, which indicates
kidney function, only weak patterns are observed. The phosphorus balance is
slightly negatively associated with the treatment effect.

For the survival time, plotting only the treatment effect against a variable
is not meaningful since the interpretation of the treatment effect depends on
the shape of the baseline hazard. Therefore, we took a different approach in
this case and show on the $y$-axis the difference in median survival between
treatment and control intake. For example, a value of 70 means that based on
the personalised model of this patient, the median survival is prolonged by 70
days if the patient takes Riluzole. The difference in median survival is
denoted by $\Delta_{0.5}$. Any other quantile could be used as well since from
the Weibull model, information on the entire estimated distribution in the two
treatment groups is obtained. Taking the difference in medians makes sense
because it is a measure on the scale of the outcome, just as the treatment
effect in a linear model, which is the difference in means. The shape of
$\Delta_{0.5}$ when plotted against age shows a strong pattern that indicates
that age is a predictive factor (see Figure 6). The treatment efficacy
increases with age until about 55 years and then flattens The difference in
median survival slightly increases with the days between disease onset and
start of treatment in the beginning, but decreases again after about 1000
days. Patients who suffer from weakness have a greater variance in their
benefit from Riluzole. Tall patients are predicted to benefit little on
average.

###  3.4 Variable importance

Figure 7: Variable importances of all split variables used for the ALSFRS
forest. Figure 8: Variable importances of all split variables used for the
survival time forest.

Figures 7 and 8 show the variable importance of each split variable. Figure 7
suggests that the time between disease onset and start of treatment plays the
most important role for the personalised models. The time between disease
onset and start of treatment, the forced vital capacity (FVC), and the
phosphorus balance have been shown to be the most important variables for
stratified models (Seibold _et al._ , 2015) which is underlined by this
analysis. The time between disease onset and start of treatment contains
information on the state of disease progression for patients in the trial. If
the disease onset and the start of treatment are far apart, the patient is
likely to have a slow progression (Hothorn and Jung, 2014). Also Riluzole has
been shown to not be effective when the disease is already far progressed
(European Medicines Agency, 2012). Thus it is not surprising that this
variable is selected as an important variable.

For the Riluzole effect on the survival time the patient’s age and again the
time between onset and treatment start play a role. Both variables have been
identified before (Seibold _et al._ , 2015) as important factors for survival
time.

##  4 Discussion

Model-based forests can find important predictive and prognostic patient
characteristics and – more importantly – via the personalised models provide
the possibility to predict the treatment effect of a future patient. Through
analysis of the PRO-ACT data and simulations (see appendix B), we showed that
personalised models can perform better than the standard global model if there
are differences in treatment effect between patients. If there is no
difference, the performance of the methods is about the same. These results
allow a shift from standardised medicine back to personalised medicine, but
this time in a controlled way by using statistical principles.

The presented methods are based on tree-based subgroup analyses but go a step
further. Not only are subgroups identified and the treatment effect within
each group estimated, but many slightly varying trees are used to retrieve a
measure of similarity between patients. On this basis, a model is computed in
which more similar patients are weighted higher. The personalised models
provide point estimates for the treatment effect. When the individual
treatment effects are plotted against patient characteristics researchers can
determine on whether the patient characteristics are predictive factors and in
what way the patient characteristics and the treatment effect are interacting.
For ALS patients, the FVC value was predictive on the ALSFRS, and the
patient’s age and height were predictive on survival. The next step would be
to generate hypotheses from these findings and plan a study to test these. Our
method offers a promising means of providing individual treatment effect
predictions and can be applied to any clinical trial data where baseline
patient characteristics are available.

All results were obtained solely using open-source implementation software
(see Section 5), which provides easy access to the methods.

##  5 Computational details

The code for data preprocessing of the PRO-ACT data is available in the
TH.data package (Hothorn, 2014). The source code for the full analyses is
available on https://github.com/HeidiSeibold/personalised_medicine.
Implementation of all methods discussed in this paper are based on the R
partykit package (Hothorn and Zeileis, 2015, version 1.0-2). Other R packages
used were sandwich (Zeileis, 2004, 2006, 2.3-3), survival (Therneau, 2015,
2.38-1), eha (Broström, 2014, 2.4-2) and ggplot2 (Wickham, 2009, 2.0.0). All
computations were conducted in the R system for statistical computing (R Core
Team, 2016, version 3.2.0).

## Acknowledgements

We thank the Swiss National Fund for funding this project with grant
205321_163456.

## References

  * Atassi _et al._ (2014) Atassi N, Berry J, Shui A, Zach N, Sherman A, Sinani E, Walker J, Katsovskiy I, Schoenfeld D, Cudkowicz M, Leitner M (2014).  “The PRO-ACT Database: Design, Initial Analyses, and Predictive Features.”  _Neurology_ , 83(19), 1719–1725.  10.1212/WNL.0000000000000951.
  * Breiman (2001) Breiman L (2001).  “Random Forests.”  _Machine Learning_ , 45(1), 5–32.  10.1023/A:1010933404324.
  * Brooks _et al._ (1996) Brooks BR, Sanjak M, Ringel S, England J, Brinkmann J, Pestronk A, Florence J, Mitsumoto H, Szirony K, Wittes J, Charatan M, Stambler N, Cedarbaum JM (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.
  * Broström (2014) Broström G (2014).  _eha: Event History Analysis_.  R package version 2.4-2, URL http://CRAN.R-project.org/package=eha.
  * Ciampi _et al._ (1995) Ciampi A, Negassa A, Lou Z (1995).  “Tree-Structured Prediction for Censored Survival-Data and the Cox Model.”  _Journal of Clinical Epidemiology_ , 48(5), 675–689.  10.1016/0895-4356(94)00164-L.
  * Dusseldorp and Van Mechelen (2013) Dusseldorp E, Van Mechelen I (2013).  “Qualitative Interaction Trees: A Tool to Identify Qualitative Treatment-Subgroup Interactions.”  _Statistics in Medicine_ , 33(2), 219–237.  10.1002/sim.5933.
  * European Medicines Agency (2012) European Medicines Agency (2012).  “Riluzole Zentiva: EPAR Summary for the Public.”  URL http://www.ema.europa.eu/docs/en_GB/document_library/EPAR_-_Summary_for_the_public/human/002622/WC500127609.pdf.
  * Foster _et al._ (2015) Foster JC, Taylor JMG, Kaciroti N, Nan B (2015).  “Simple subgroup approximations to optimal treatment regimes from randomized clinical trial data.”  _Biostatistics_ , 16(2), 368–382.  10.1093/biostatistics/kxu049.
  * Hastie _et al._ (2009) Hastie T, Tibshirani R, Friedman J (2009).  _The Elements of Statistical Learning_.  Springer Series in Statistics. Springer, 2nd edition.
  * Hothorn (2014) Hothorn T (2014).  _TH.data: TH’s Data Archive_.  R package version 1.0-4, URL http://CRAN.R-project.org/package=TH.data.
  * Hothorn _et al._ (2006a) Hothorn T, Hornik K, Van de Wiel MA, Zeileis A (2006a).  “A Lego System for Conditional Inference.”  _The American Statistician_ , 60(3), 257–263.  10.1198/000313006X118430.
  * Hothorn _et al._ (2006b) Hothorn T, Hornik K, Zeileis A (2006b).  “Unbiased Recursive Partitioning: A Conditional Inference Framework.”  _Journal of Computational and Graphical Statistics_ , 15(3), 651–674.  10.1198/106186006X133933.
  * Hothorn and Jung (2014) Hothorn T, Jung HH (2014).  “RandomForest4Life: A Random Forest for Predicting ALS Disease Progression.”  _Amyotrophic Lateral Sclerosis & Frontotemporal Degeneration_, 15(5-6), 444–452.  10.3109/21678421.2014.893361.
  * 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 and Zeileis (2015) Hothorn T, Zeileis A (2015).  “partykit: A Modular Toolkit for Recursive Partytioning in R.”  _Journal of Machine Learning Research_ , 16, 3905–3909.  URL http://jmlr.org/papers/v16/hothorn15a.html.
  * Italiano (2011) Italiano A (2011).  “Prognostic or Predictive? It’s Time to Get Back to Definitions!”  _Journal of Clinical Oncology_ , 29(35), 4718–4718.  10.1200/JCO.2011.38.3729.
  * Kehl and Ulm (2006) Kehl V, Ulm K (2006).  “Responder Identification in Clinical Trials with Censored Data.”  _Computational Statistics & Data Analysis_, 50(5), 1338–1355.  10.1016/j.csda.2004.11.015.
  * Loh _et al._ (2015) Loh WY, He X, Man M (2015).  “A Regression Tree Approach to Identifying Subgroups with Differential Treatment Effects.”  _Statistics in Medicine_.  10.1002/sim.6454.
  * Meinshausen (2006) Meinshausen N (2006).  “Quantile Regression Forests.”  _Journal of Machine Learning Research_ , 7, 983–999.  URL http://www.jmlr.org/papers/v7/meinshausen06a.html.
  * R Core Team (2016) R Core Team (2016).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL https://www.R-project.org/.
  * Seibold _et al._ (2015) Seibold H, Zeileis A, Hothorn T (2015).  “Model-based Recursive Partitioning for Subgroup Analyses.”  _International Journal of Biostatistics_.  10.1515/ijb-2015-0032.  Accepted 2015-10-02.
  * Strasser and Weber (1999) Strasser H, Weber C (1999).  “On the Asymptotic Theory of Permutation Statistics.”  _Mathematical Methods of Statistics_ , 8, 220–250.  URL http://epub.wu.ac.at/102/1/document.pdf.
  * Therneau (2015) Therneau TM (2015).  _A Package for Survival Analysis in S_.  Version 2.38, URL http://CRAN.R-project.org/package=survival.
  * Tian _et al._ (2014) Tian L, Alizadeh AA, Gentles AJ, Tibshirani R (2014).  “A Simple Method for Estimating Interactions Between a Treatment and a Large Number of Covariates.”  _Journal of the American Statistical Association_ , 109(508), 1517–1532.  10.1080/01621459.2014.951443.
  * Weisberg (2015) Weisberg HI (2015).  “What next for randomised clinical trials?”  _Significance_ , 12(1), 22–27.  10.1111/j.1740-9713.2015.00798.x.
  * Wickham (2009) Wickham H (2009).  _ggplot2: Elegant Graphics for Data Analysis_.  Springer-Verlag New York.  ISBN 978-0-387-98140-6.  URL http://had.co.nz/ggplot2/book.
  * Zeileis (2004) Zeileis A (2004).  “Econometric Computing with HC and HAC Covariance Matrix Estimators.”  _Journal of Statistical Software_ , 11(10), 1–17.  10.18637/jss.v011.i10.
  * Zeileis (2006) Zeileis A (2006).  “Object-Oriented Computation of Sandwich Estimators.”  _Journal of Statistical Software_ , 16(9), 1–16.  10.18637/jss.v016.i09.
  * Zeileis and Hornik (2007) Zeileis A, Hornik K (2007).  “Generalized M-fluctuation Tests for Parameter Instability.”  _Statistica Neerlandica_ , 61(4), 488–508.  10.1111/j.1467-9574.2007.00371.x.
  * 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.
  * Zhang _et al._ (2012) Zhang B, Tsiatis AA, Laber EB, Davidian M (2012).  “A Robust Method for Estimating Optimal Treatment Regimes.”  _Biometrics_ , 68(4), 1010–1018.  10.1111/j.1541-0420.2012.01763.x.

##  Appendix A Split algorithm in detail

In the following, the split algorithm in model-based recursive partitioning is
explained. The split procedure starts with all $N$ data points. In nodes other
than the root node, the size of the data set depends on the previous splits.
For notational simplicity, we describe the split procedure in the root node,
i.e. for patients $i=1,\dots,N$.

  * •

Compute prespecified (parametric) model
$\mathcal{M}((Y,\mathbf{X}),\boldsymbol{\vartheta})$. Estimate
$\hat{\boldsymbol{\vartheta}}$ by maximising the log-likelihood

$\hat{\boldsymbol{\vartheta}}=\operatornamewithlimits{argmax}\limits_{\vartheta}l((Y,\mathbf{X}),\boldsymbol{\vartheta})$

or equivalently by solving

$\sum\limits_{i=1}^{N}s((y,\mathbf{x})_{i},{\boldsymbol{\vartheta}})=0$

for $\boldsymbol{\vartheta}$.

  * •

Compute associated empirical estimating function (residuals)

$\boldsymbol{s}=\begin{pmatrix}s_{\hat{\alpha}}((y,\mathbf{x})_{1},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{1},\hat{\boldsymbol{\vartheta}})\\\
s_{\hat{\alpha}}((y,\mathbf{x})_{2},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{2},\hat{\boldsymbol{\vartheta}})\\\
\vdots&\vdots\\\
s_{\hat{\alpha}}((y,\mathbf{x})_{N},\hat{\boldsymbol{\vartheta}})&s_{\hat{\beta}}((y,\mathbf{x})_{N},\hat{\boldsymbol{\vartheta}})\\\
\end{pmatrix}$

  * •

Perform tests of independence between residuals (partial score vectors)
$\boldsymbol{s}_{\hat{\alpha}}$ as well as $\boldsymbol{s}_{\hat{\beta}}$ and
partitioning variables $Z_{j}$.

$\displaystyle H_{0}^{\alpha,j}$
$\displaystyle:\quad\boldsymbol{s}_{\hat{\alpha}}((Y,\mathbf{X}),\hat{\boldsymbol{\vartheta}})\quad\bot\quad
Z_{j}$ $\displaystyle H_{0}^{\beta,j}$
$\displaystyle:\quad\boldsymbol{s}_{\hat{\beta}}((Y,\mathbf{X}),\hat{\boldsymbol{\vartheta}})\quad\bot\quad
Z_{j}\qquad j=1,\dots,J$

For the tests, we use permutation testing with the linear statistic

$T_{j}=\sum\limits_{i\in\mathcal{B}_{b}}g_{j}(Z_{ji})\cdot\mathbf{s}_{i}$

The transformation function $g$ depends on the scale of the variable $Z_{j}$.
If $Z_{j}$ is numeric then $g_{j}(z_{ji})=z_{ji}$. If $Z_{j}$ is categorical
with $K$ categories then
$g_{j}(z_{ji})=e_{K}(z_{ji})=(I(z_{ji}=1),\dots,I(z_{ji}=K))$, i.e., $g_{j}$
is the unit vector of length $K$, where the element, that corresponds to the
value of $z_{ji}$, is one. If there are missing values in $Z_{j}$ the
observations are excluded from the sum so that we actually sum over all
observations $i\in\mathcal{B}_{b}$, except for the observations in
$\mathcal{B}_{b}$, where $Z_{j}$ is missing. The standardised test statistic
is the Pearson correlation coefficient

$c({t}_{j},\mu_{j},\Sigma_{j})=\left|\frac{({t}_{j}-\mu_{j})}{\sqrt{(\Sigma_{j})}}\right|$

if $Z_{j}$ is numeric and otherwise

$c(\mathbf{t}_{j},\mu_{j},\Sigma_{j})=\max\limits_{k=1,\dots,K}\left|\frac{(\mathbf{t}_{j}-\mu_{j})_{k}}{\sqrt{(\Sigma_{j})_{kk}}}\right|.$

The conditional expectation $\mu_{j}$ and covariance $\Sigma_{j}$ can be
derived as in Strasser and Weber (1999). The smallest $p$-value corresponds to
the largest discrepancy from the model assumption, that intercept and
treatment parameter are the same for all patients in the given node/subgroup.

  * •

If any Bonferroni adjusted $p$-value is lower than the significance level,
select the partitioning variable $Z_{j*}$ that has the highest association
(lowest $p$-value) to any of the residuals relevant for the split.

  * •

Select as split point the point that results in the largest discrepancy
between score functions in the two resulting subgroups. The discrepancy can be
measured by the linear statistic

$T_{j*}^{k}=\sum\limits_{i\in\mathcal{B}_{1k}}\mathbf{s}_{i},$

where $\mathcal{B}_{1k}$ here is the first of the two new subgroups that are
defined by splitting in split point $k$ of variable $Z_{j*}$. The split point
is then chosen as follows:

$k*=\operatornamewithlimits{argmin}\limits_{k}c(t_{j*}^{k},\mu_{j*}^{k},\Sigma_{j*}^{k}).$

##  Appendix B Empirical evaluation

To check whether the proposed method can recover smooth treatment effect
functions, we evaluated its performance on artificial data. To do so, we
simulated data from a normal linear regression model. We simulated ten
correlated patient characteristics, where only one is in a non-linear
interaction with the treatment. In the following, we compare the log-
likelihood of our method to the log-likelihood of the true underlying model
and the naive model that assumes an overall applicable treatment effect
(Section B.1) and show the predicted treatment effects in dependence plots
(Section B.2) and the variable importances of the true predictive factor and
the noise variables (Section B.3).

We simulated 600 patients, half of which were treated ($x_{A}=1$) and half of
which were untreated ($x_{A}=0$). The ten partitioning variables $\mathbf{Z}$
are normally distributed

$\displaystyle\mathbf{Z}$
$\displaystyle\sim\mathcal{N}_{10}(\mathbf{0},\mathbf{\Sigma}_{Z})$ (24)

and correlated with the covariance matrix

$\displaystyle\mathbf{\Sigma}_{Z}$
$\displaystyle=\begin{pmatrix}1&0.2&\cdots&0.2\\\ 0.2&1&\cdots&0.2\\\
\vdots&\vdots&\ddots&\vdots\\\ 0.2&0.2&\cdots&1\end{pmatrix}.$

The primary outcome depends on treatment and partitioning variables as follows

$\displaystyle Y|\mathbf{X}=\mathbf{x},\mathbf{Z}=\mathbf{z}$
$\displaystyle\sim\mathcal{N}(1.9+0.2\cdot x_{A}+3\cdot\cos(z_{1})\cdot
x_{A},1).$ (25)

In this example the true model parameters are defined as follows:

$\displaystyle\alpha(\mathbf{z})$ $\displaystyle=1.9$ (26)
$\displaystyle\beta(\mathbf{z})$ $\displaystyle=0.2+0.3\cdot\cos(z_{1}).$

Figure 9: True treatment effect in given simulated data.

This means that the treatment effect depends on the value of $z_{1}$ and this
dependency has the form of a cosinus function (see Figure 9).

###  B.1 Comparison of models by comparing the log-likelihood

To compare our method with (i) a correctly specified model taking into account
the main effects of $x_{A}$ and $\cos(z_{1})$ as well as the interaction of
$x_{A}$ and $\cos(z_{1})$ and (ii) a simple linear model including only the
treatment $x_{A}$ as a covariate, we drew 100 learning samples and 100 test
samples using the data simulation procedure explained above and computed the
out-of-sample log-likelihoods (i.e. based on the test data) for the models
after applying them to each of the 100 learning data sets. The log-likelihood
contributions

$\displaystyle
l\left((y,\mathbf{x})_{i},\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right)=\left(y_{i}-\mathbf{x}_{i}^{\top}\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})\right)^{2}$
(27)

with $\mathbf{x}_{i}=(1,x_{iA})^{\top}$ and
$\hat{\boldsymbol{\vartheta}}(\mathbf{z}_{i})=(\hat{\alpha}(\mathbf{z}_{i}),\hat{\beta}(\mathbf{z}_{i}))^{\top}$
are taken from the personalised models of our method (see Section 2.5). Note
that for the simple linear model the log-likelihood contributions are defined
as above, but only with constant parameters, for the fully specified model
$\mathbf{x}_{i}=(1,x_{A},\cos(z_{1}),x_{A}\cdot\cos(z_{1}))^{\top}$ and
$\mathbf{\hat{\vartheta}}=(\hat{\alpha},\hat{\beta}_{A},\hat{\beta}_{\cos(z_{1})},\hat{\beta}_{A,\cos(z_{1})})^{\top}$.

Figure 10: Out-of-sample log-likelihoods obtained from the three models. Each
line represents one simulated data set.

The log-likelihoods of our method are higher than the log-likelihoods of the
simple and incorrect linear model and lower than the log-likelihoods of the
correctly specified model (Figure 10). Therefore, we conclude that our method
performs reasonably well.

###  B.2 Dependence plots

(a) Dependence plot for $z_{1}$.

(b) Dependence plot for $z_{2}$

Figure 11: Joint dependence plots of all 100 simulations.

For the same 100 simulated test data sets as above, we obtained the dependence
plots. Figure 11 shows two dependence plots in which all 100 simulations are
combined by layering them on top of each other. The dependence plot of
partitioning variable $z_{1}$ (Figure 11(a)) shows a curve that is fairly
similar to that of Figure 9, except that the effect is shrunken towards zero.
Note that with a larger sample or differently tuned parameters (e.g. larger
trees), one could get better results for the extreme treatment effects. As
expected, for partitioning variables $z_{2}$ to $z_{10}$, there is only random
fluctuation around zero (see as an example Figure 11(b), which shows the
dependence plot for $z_{2}$).

###  B.3 Variable importance

Figure 12: Variable importances of the predictive factor $z_{1}$ and noise
variables $z_{2}$ to $z_{10}$.

Variable importances for one simulated data set are shown in Figure 12. As
expected, partitioning variable $z_{1}$ is the only variable with a clearly
positive variable importance. Even though all partitioning variables are
correlated, the method was able to distinguish between the correlation and
predictive effect.