# Trees-Based Models for Correlated Data

Assaf Rabinowicz<EMAIL_ADDRESS>
Department of Statistics and Operations Research
Tel Aviv University
Tel Aviv, Israel Saharon Rosset<EMAIL_ADDRESS>
Department of Statistics and Operations Research
Tel Aviv University
Tel Aviv, Israel

###### Abstract

This paper presents a new approach for regression trees-based models, such as
simple regression tree, random forest and gradient boosting, in settings
involving correlated data. We show the problems that arise when implementing
standard regression trees-based models, which ignore the correlation
structure. Our new approach explicitly takes the correlation structure into
account in the splitting criterion, stopping rules and fitted values in the
leaves, which induces some major modifications of standard methodology. The
superiority of our new approach over trees-based models that do not account
for the correlation is supported by simulation experiments and real data
analyses.

Keywords: random forest, linear mixed models, Gaussian process regression,
prediction error for correlated data, model selection

##  1 Introduction

Trees-based models are widely used for tabular data due to their high
prediction accuracy and their inherent model selection functionality (Hastie
et al., 2009). Commonly, trees-based models are fitted without assuming any
distributional setting on the dependent variable. While the distribution of
the dependent variable is mostly unknown and therefore it is tempting to avoid
distributional assumptions, the correlation structure, which relates to the
sampling mechanism (e.g., clustered data, time-series data, longitudinal data,
spatial data), is frequently known and therefore it is not reasonable to
ignore it. Unlike in trees-based models, the correlation structure is an
essential component in many machine learning models, for example kernel
covariance functions are used in Gaussian processes regression (Rasmussen,
2003), which is frequently implemented for modeling data sets with spatial
correlation structure, such as neuroscience data sets (Caywood et al., 2017)
and climatography data sets (Goovaerts, 1999). Another example is linear mixed
model, which is used for data involving longitudinal correlation structure, as
is common in health (Coull et al., 2001) and trading (Westveld et al., 2011)
applications.

In this paper we develop a method which combines the concepts of random
effects and random fields, which are convenient platforms for analyzing
correlated data, and trees-based models such as: regression tree, random
forest and gradient boosting. The desired result is that the trees-based part
results a high prediction accuracy and model selection capabilities and the
random effects part enables to boost the model performance by using the
correlation structure and even to exract statistical inference. The idea of
integrating between random effects/random field and trees-based methods has
previously been explored (see Sela and Simonoff 2012; Stephan et al. 2015,
discussed in more detail in Section 4). However, we propose a novel approach
which takes advantage of recent developments in model evaluation and selection
methodologies for correlated settings, and yields improved results as
demonstrated below.

Section 2 gives relevant background for the proposed method. The background
contains a brief description of trees-based methods, linear mixed model (which
is based on random effects) and prediction error estimation for correlated
data, which has a key role in our approach. Our new algorithm, REgression Tree
for COorrelated data (RETCO), is presented and discussed in Section 3. Section
4 compares RETCO with other algorithms that were introduced in recent years.
Section 5 presents simulation and real data analyses that support our proposed
algorithm.

##  2 Theoretical Background

This section presents briefly regression trees-based models, linear mixed
models and prediction error estimation for correlated data. Additional
information can be found in Appendix A and in resources which are cited below.

###  2.1 Trees-Based Models

Given a vector of covariates $\boldsymbol{x}^{*}\in\mathbb{R}^{p}$ a
regression tree estimates the corresponding response, $y^{*}\in\mathbb{R},$ as
follows:

$f(\boldsymbol{x}^{*})=\sum_{s=1}^{S}I_{(\boldsymbol{x}^{*}\in
g_{s})}\mu_{s},$

where $g_{s}\subseteq\mathbb{R}^{p}$ and $\mu_{s}\in\mathbb{R},$ for all
$s\in\\{1,...,S\\}.$ $\\{g_{s}\\}_{s=1}^{S}$ define a partition of the
covariate space, that is $g_{s}\cap g_{t}=\emptyset$ for $s\neq t$ and
$\cup_{s=1}^{S}g_{s}$ is the entire covariates space. $\mu_{s}$ is the
predictor for covariate vectors that are in $g_{s}.$ The nodes,
$\mathcal{S}=\\{\mu_{s},g_{s}\\}_{s=1}^{S}$ are selected using a recursive
optimization process that can be diagrammed as a tree, where $\mathcal{S}$ are
the terminal nodes (the leaves). The recursive optimization for selecting
$\mathcal{S}$ is based on the training set,
$\\{y_{i},\boldsymbol{x}_{i}\\}_{i=1}^{n}=\\{\boldsymbol{y},X\\},$ where
$\boldsymbol{y}\in\mathbb{R}^{n}$ and $X\in\mathbb{R}^{n\times p}$ is the
design matrix, and it is implicitly assumed that the prediction points,
$y^{*}$ and $\boldsymbol{x}^{*},$ are drawn from the same distribution as
$y_{i}$ and $\boldsymbol{x}_{i}.$ Each step in the tree’s recursive
optimization process is a model selection of linear models, where a threshold
of one of the available covariates is selected in order to minimize a loss
function, $Loss(\cdot,\cdot):\mathbb{R}^{2}\to\mathbb{R}.$ Also, stopping
rules that follow predefined hyper-parameters, e.g., maximal depth of the tree
and minimal number of training set observations in a node, are enforced on the
recursive optimization and limit the tree’s depth and affect other properties
of the tree’s structure in order to avoid overfitting.

Random forest (RF) and gradient boosting (GB) predictor are based on averaging
an ensemble of trees. More information about regression tree, RF and GB, as
well a formalized regression tree algorithm, is available in Appendix A and
can also be found in Freund et al. (1999); Friedman (2001); Breiman (2001);
Hastie et al. (2009).

###  2.2 Linear Mixed Model

In linear mixed models (LMM) there are two covariate vectors: fixed effects
covariates, $\boldsymbol{x}^{*}\in\mathbb{R}^{p},$ and random effects
covariates, $\boldsymbol{z}^{*}\in\mathbb{R}^{q}.$ Commonly, $y^{*}$ is
assumed to be normally distributed and decomposed as follows:

$\displaystyle
y^{*}=\boldsymbol{\beta}^{t}\boldsymbol{x}^{*}+\boldsymbol{b}^{t}\boldsymbol{z}^{*}+\epsilon^{*},$

where $\boldsymbol{\beta}$ is the fixed effects vector of coefficients,
$\boldsymbol{b}\sim N_{q}(0,G)$ is the random effects vector and
$\epsilon^{*}\sim N(0,\sigma^{2})$ is the residual.
$\boldsymbol{\beta}^{t}\boldsymbol{x}^{*}$ is the marginal mean of $y^{*},$
$\mathbb{E}(y^{*}|\boldsymbol{x}^{*}),$ and
$\boldsymbol{\beta}^{t}\boldsymbol{x}^{*}+\boldsymbol{b}^{t}\boldsymbol{z}^{*}$
is the conditional mean of $y^{*}$ given $\boldsymbol{b},$ which is commonly
denoted as
$\mathbb{E}(y^{*}|\boldsymbol{x}^{*},\boldsymbol{z}^{*},\boldsymbol{b}).$

$\boldsymbol{\beta}$ and $\boldsymbol{b}$ are estimated using the training
sample,
$\\{y_{i},\boldsymbol{x}_{i},\boldsymbol{z}_{i}\\}_{i=1}^{n}=\\{\boldsymbol{y},X,Z\\},$
which follows the same model:

$\boldsymbol{y}=X\boldsymbol{\beta}+Z\boldsymbol{b}+\boldsymbol{\epsilon},$

where $\boldsymbol{\epsilon}\sim N_{n}(0,\sigma^{2}I_{n}),$ and $X,Z$ are the
fixed effects and random effects covariate matrices, respectively. Since
$y^{*}$ and $\boldsymbol{y}$ share the same random effects, $\boldsymbol{b},$
they are correlated, and estimating $\boldsymbol{b}$ by the training sample
can later be used for improving the prediction accuracy of $y^{*}.$ Given $G$
and $V\coloneqq\mathrm{Var}(\boldsymbol{y})=ZGZ^{t}+\sigma^{2}I_{n},$
$\boldsymbol{\beta}$ and $\boldsymbol{b}$ can be estimated as follows:

$\displaystyle\widehat{\boldsymbol{\beta}}$
$\displaystyle=(X^{t}V^{-1}X)^{-1}X^{t}V^{-1}\boldsymbol{y}$
$\displaystyle\widehat{\boldsymbol{b}}$
$\displaystyle=GZV^{-1}(\boldsymbol{y}-X\widehat{\boldsymbol{\beta}}).$

Harville et al. (1976) showed that given the true covariance matrices, $G$ and
$V,$ the estimated conditional mean,
$\widehat{\mathbb{E}}(y^{*}|\boldsymbol{x}^{*},\boldsymbol{z}^{*},\widehat{\boldsymbol{b}})=\widehat{\boldsymbol{\beta}}^{t}\boldsymbol{x}^{*}+\widehat{\boldsymbol{b}}^{t}\boldsymbol{z}^{*},$
is the best linear unbiased predictor (BLUP) of $y^{*}.$ In practice, the
covariance matrices are mostly unknown and therefore are estimated using
maximum likelihood or restricted maximum likelihood (Verbeke, 1997). Note that
given the covariance matrices, LMM is linear in $\boldsymbol{y},$ i.e., the
LMM predictor satisfies:

$\widehat{y}^{*}=h^{*}\boldsymbol{y},$

where $h^{*},$ the hat vector, does not contain the training response vector
$\boldsymbol{y},$ and the element $h^{*}[i]$ is the weight of $y_{i}$ in
predicting $y^{*},\forall i\in[1,...,n].$

####  2.2.1 $\boldsymbol{b}^{*}\neq\boldsymbol{b}$ Scenario

In many cases the random effects of $y^{*}$ are not the same as the random
effects of $\boldsymbol{y},$ i.e.,

$y^{*}=\boldsymbol{\beta}^{t}\boldsymbol{x}^{*}+\boldsymbol{b}^{*{t}}\boldsymbol{z}^{*}+\epsilon^{*},\;\;\boldsymbol{b}^{*}\neq\boldsymbol{b}.$

That means that the correlation between the observations in $\boldsymbol{y}$
is not the same as the correlation between $y^{*}$ and the observations in
$\boldsymbol{y}.$ In case $\boldsymbol{b}^{*}\perp\boldsymbol{b},$ which
implies $\mathrm{Cov}(y^{*},\boldsymbol{y})=0,$ estimating $\boldsymbol{b}$
does not improve the prediction accuracy of $y^{*}.$ Therefore, in this case
$y^{*}$ is predicted by the marginal mean
$,\widehat{\mathbb{E}}(y^{*}|\boldsymbol{x}^{*})=\widehat{\boldsymbol{\beta}}^{t}\boldsymbol{x}^{*}.$
This model, which is a special case of LMM, is also called generalized least
squares model (GLS). A simple example for this scenario, is when the
prediction set contains different clusters than in the training set. Several
prediction tasks that follow this scenario setting are analyzed in Section
5.2. For example, using FIFA data set from Kaggle website, a predictive model
for footballs players’ market-value was trained using GLS, where the
prediction goal is to predict the market-value of players that belong to clubs
that do not appear in the training set. This data set has a clustered
correlation structure, where cluster is the players’ club, i.e., market-values
of players that belong to the same club are correlated, while market-values of
players that belong to different clubs are uncorrelated. Given the prediction
goal of predicting the market-value of players from new clubs, this setting
follows exactly the $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ setting.

Another scenario is when $\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$ (although
$\boldsymbol{b}^{*}\neq\boldsymbol{b}$) and therefore
$\mathrm{Cov}(y^{*},\boldsymbol{y})\neq 0.$ This can happen for example when
some of the random effects of $\boldsymbol{y}$ and $y^{*}$ are the same and
some are not. In this case, the elements in $\boldsymbol{b}$ that are in
$\boldsymbol{b}^{*}$ should be estimated and used for predicting $y^{*}.$
Examples where $\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$ although
$\boldsymbol{b}^{*}\neq\boldsymbol{b},$ and more information about this
scenario can be found in Rabinowicz and Rosset (2020). This scenario, of
$\boldsymbol{b}^{*}\neq\boldsymbol{b},$ is common and should be taken into
account when developing trees-based methods for correlated data, as will be
developed in Section 3.

###  2.3 Prediction Error Estimation and Model Selection for Correlated Data

Once a predictive model is fitted, it is often evaluated by its prediction
error estimator. Moreover, when there is a set of alternative models (e.g.,
for LMM: models with different covariates, for trees-based models: models with
different hyper-parameters) the ’best’ model can be selected based on
minimizing the prediction error estimator. It is important to note that common
prediction error estimators, e.g., AIC (Akaike, 1974), Cp (Mallows, 1973), and
even cross-validation (Stone, 1974, CV) are biased in some settings involving
correlated data. Naturally, their corresponding model selection criteria are
also biased in those scenarios. This bias was studied in the recent years,
mostly for linear models. Here we present Cp, AIC and CV versions for
correlated data. For description of the original Cp and AIC versions, that do
not address correlation structure, see Appendix A.

####  2.3.1 Cp

In Cp, the goal is to estimate the squared prediction error:

$\small\mathbb{E}_{\boldsymbol{y},\boldsymbol{y}^{*}}\frac{1}{n}\|\boldsymbol{y}^{*}-H\boldsymbol{y}\|^{2}_{2},$

where $H$ is the hat matrix, and $\boldsymbol{y}^{*}\in\mathbb{R}^{n}$ is a
vector of new observations measured at the same covariate values as
$\boldsymbol{y},$ $\\{X,Z\\},$ but with new independent noise and potentially
different random effects realizations. This type of prediction error, when
both $\boldsymbol{y}$ and $\boldsymbol{y}^{*}$ relate to the same covariate
points, $\\{X,Z\\},$ is called in-sample prediction error. In this setting, it
is natural to consider $\\{X,Z\\}$ as fixed matrices rather than random
variables. Hodges and Sargent (2001) extended Cp to LMM with
$\boldsymbol{b}^{*}=\boldsymbol{b},$ here we employ a more general formulation
which reduces to Hodges and Sargent (2001) when
$\boldsymbol{b}^{*}=\boldsymbol{b}$ but also covers the case that they are
different:

$\displaystyle\small
Cp=\frac{1}{n}\|\boldsymbol{y}-\widehat{\boldsymbol{y}}\|_{2}^{2}+\frac{2}{n}\mathrm{tr}\Big{(}H\big{(}\mathrm{Var}(\boldsymbol{y})-\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})\big{)}\Big{)}.$

####  2.3.2 AIC

Similarly to Cp, AIC is also an in-sample error, however its loss function is
based on likelihood. Vaida and Blanchard (2005) presented the conditional AIC
(cAIC) and marginal AIC (mAIC) which are suitable for the scenarios where
$\boldsymbol{b}^{*}=\boldsymbol{b}$ and
$\boldsymbol{b}^{*}\perp\boldsymbol{b},$ respectively. Here we will use the
name AIC for our formulation which subsumes cAIC and mAIC, but also covers the
$\boldsymbol{b}^{*}\neq\boldsymbol{b}\cap\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$
scenario:

$\displaystyle\small AIC=$
$\displaystyle-\frac{2\ell\big{(}\boldsymbol{y};\widehat{\mathbb{E}}(\boldsymbol{y}|X,Z,\widehat{\boldsymbol{b}}),V_{c}\big{)}}{n}+\frac{2\mathrm{tr}\Big{(}H\big{(}\mathrm{Var}(\boldsymbol{y})-\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}^{*})\big{)}V^{-1}_{c}\Big{)}}{n},$

where
$\ell(\boldsymbol{y};\widehat{\mathbb{E}}(\boldsymbol{y}|X,Z,\widehat{\boldsymbol{b}}),V_{c})$
is the conditional likelihood of $\boldsymbol{y}$ given
$\widehat{\boldsymbol{b}},$ and
$V_{c}=\mathrm{Var}(\boldsymbol{y}^{*}|\boldsymbol{b}).$

####  2.3.3 Cross-Validation (CV)

Unlike Cp and AIC, CV estimates the generalization error:

$\mathbb{E}_{X,X^{*},Z,Z^{*}}\mathbb{E}_{\boldsymbol{y}^{*},\boldsymbol{y}}\sum_{i=1}^{n}\frac{1}{n}Loss\big{(}y^{*}_{i},\widehat{y}_{i}^{*}(\boldsymbol{x}_{i}^{*},\boldsymbol{z}_{i}^{*};\boldsymbol{y},X,Z)\big{)},$

where $\boldsymbol{y}^{*}$ is drawn from the same marginal distribution as
$\boldsymbol{y},$ but relates to new covariate values,
$\\{\boldsymbol{x}^{*}_{i},\boldsymbol{z}^{*}_{i}\\}_{i=1}^{n}=\\{X^{*},Z^{*}\\},$
which were drawn from the same distribution as $\\{X,Z\\}.$ For simplicity, a
special case of CV algorithm, leave-one-out (LOO), is presented:

  1. 1.

$\forall i\in[1,...,n],$ fit a model using the whole sample besides the
$i^{th}$ observation. Denote the sample without the $i^{th}$ observation as
$\\{\boldsymbol{y}_{-i},X_{-i},Z_{-i}\\}.$

  2. 2.

Predict $y_{i}$ by the fitted model and denote the predictor as
$\widehat{y}_{i}^{-i}=\widehat{y}_{i}(\boldsymbol{x}_{i},\boldsymbol{z}_{i};\boldsymbol{y}_{-i},X_{-i},Z_{-i}).$

For a squared errors loss function and linear predictor of $\boldsymbol{y},$
the CV error is

$CV=\frac{1}{n}\|\boldsymbol{y}-H_{cv}\boldsymbol{y}\|_{2}^{2},$

where $H_{cv},$ the CV hat matrix, is

$H_{cv}=\left[\begin{array}[]{cccc}0&h_{1,2}&...&h_{1,n}\\\
h_{2,1}&0&&h_{2,n}\\\ ...\\\ h_{n,1}&h_{n,2}&...&0\end{array}\right],$

$h_{k,k^{\prime}}\in\mathbb{R}\;\forall k,k^{\prime}\in\\{1,...,n\\}.$ In this
presentation, the vector $[h_{1,2},...,h_{1,n}]$ is the hat matrix of
$\widehat{y}_{1}^{-1}.$ K-fold CV generalizes LOO, by partitioning
$\\{\boldsymbol{y},X,Z\\}$ into K equal size subsets,
$\\{\boldsymbol{y}_{k},X_{k},Z_{k}\\}_{k=1}^{K},$ where $K\leq n$ (rather than
$K=n$ as in LOO)

Rabinowicz and Rosset (2020) presented a generalization of CV, $CV_{c},$ which
is suitable for scenarios involving correlated data:

$\displaystyle\small
CV_{c}=CV+\frac{2}{n}\mathrm{tr}\Big{(}H_{cv}\big{(}\mathrm{Var}(\boldsymbol{y})-\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})\big{)}\Big{)}.$

In the LMM settings, when $\boldsymbol{b}^{*}=\boldsymbol{b},$ $CV_{c}=CV.$
For more information see Rabinowicz and Rosset (2020).

Cp and $CV_{c}$ do not assume a specific distributional setting and can be
applied for any linear model, while AIC is suitable for LMM, however it can
also be adjusted for other linear models that assume normality, such as
Gaussian process regression (GPR). Table 1 summarizes the prediction error
estimators that are described in this section.

Method | Prediction Error Type | Distributional Assumptions
---|---|---
Cp | in-sample error | –
AIC | in-sample error | normal likelihood
$CV_{c}$ | generalization error | –
Table 1: Summary of prediction errors for correlated data

##  3 Trees-Based Models for Correlated Data

This section presents the main algorithm of this paper, REgression Tree for
COrrelated data (RETCO), and discusses the main differences between RETCO and
the standard regression tree algorithm.

###  3.1 RETCO Algorithm

A simple approach for integrating between trees-based methods and random
effects is replacing the marginal mean in LMM,
$\boldsymbol{\beta}^{t}\boldsymbol{x}^{*},$ by a trees-based model,
$f(\boldsymbol{x}^{*}):$

$\displaystyle f(\boldsymbol{x}^{*})+\boldsymbol{b}^{t}\boldsymbol{z}^{*},$
(1)

where $f(\boldsymbol{x}^{*})$ is created in a way that does account for the
correlation structure (unlike in the standard regression tree algorithm). The
power of the model in expression (1) can be perceived from different points of
views. From the LMM point of view, the additive representation of marginal and
conditional means is preserved, however the marginal mean is non-linear and
therefore more expressive than in standard LMM. From the regression tree point
of view, this approach differentiates between the two types of
covariates—fixed effects, which are used for splitting the tree’s nodes, and
random effects that are added linearly to the fitted tree—enabling expressing
and using the correlation structure. This also enables using inference tools
that do not exist for standard regression trees but exists in LMM, for
example, comparing between the variance components, $\sigma^{2}$ and $G.$ Note
that expression (1) assumes that the effect of the random effects on
$\boldsymbol{y}$ is linear, however it can be generalized.

The template in expression (1) was already suggested (for literature review
see Section 4), however here we present a new algorithm, RETCO, for fitting
$f(\boldsymbol{x}^{*})$ that follows the theoretical aspects that were
presented in the previous section. The algorithm formulation is general, and
is not based on a specific prediction error type or distributional setting.
Also, the algorithm refers to a case when $f(\cdot)$ is a single regression
tree, the extension to RF and GB will be discussed in Section 3.1.1.

Algorithm 1, which presents RETCO use the following notations:

  * •

$f(\boldsymbol{x}_{i}|\mathcal{S}),$ the intermediate predictor of $y_{i}$
during the tree fitting.

$f(\boldsymbol{x}_{i}|\mathcal{S})=\sum_{s\in\mathcal{S}}I_{(\boldsymbol{x}_{i}\in
g_{l})}\eta(\mu_{s}),$

where

    * –

$\mu_{s}$ is the GLS predictor of $\\{y_{i}|\boldsymbol{x}_{i}\in
g_{s}\\},\;\forall s\in\mathcal{S}$

    * –

$\eta(\cdot)$ is the identify function for the setting of
$\boldsymbol{b}^{*}\perp\boldsymbol{b},$ and the BLUP for
$\boldsymbol{b}^{*}\not\perp\boldsymbol{b}.$

Of course, after fitting the tree
$f(\boldsymbol{x}_{i})=f(\boldsymbol{x}_{i}|\mathcal{S}).$

  * •

$f\left(\boldsymbol{x}_{i}|j,c,\mathcal{S}/s\right),$ the predictor of $y_{i}$
when splitting node $s$ using covariate $j$ at the threshold $c:$

$f\left(\boldsymbol{x}_{i}|j,c,\mathcal{S}/s\right)=I_{(\boldsymbol{x}_{i}\in
g_{s}\;\cap\;x_{i,j}\leq
c)}\eta\big{(}\mu_{s}^{l}(c)\big{)}+I_{(\boldsymbol{x}_{i}\in
g_{s}\;\cap\;x_{i,j}>c)}\eta\big{(}\mu_{s}^{r}(c)\big{)}+\sum_{l\in\mathcal{S}/s}I_{(\boldsymbol{x}_{i}\in
g_{l})}\eta(\mu_{l}),$

where $\mu_{s}^{l}(c)$ and $\mu_{s}^{r}(c)$ are the GLS predictors of
$\\{y_{i}|\boldsymbol{x}_{i}\in g_{s}\;\cap\;x_{i,j}\leq c\\}$ and
$\\{y_{i}|\boldsymbol{x}_{i}\in g_{s}\;\cap\;x_{i,j}>c\\},$ respectively.

Algorithm 1 REgression Tree for COrrelated Data (RETCO)

Input: $\boldsymbol{y},\;X,\;Z.$

Output: $f(\cdot).$

High level setting:

  * •

select a prediction error estimator loss function—Cp, AIC or $CV_{c}$—from
Table 1

  * •

define the relation between $\boldsymbol{b}^{*}$ and $\boldsymbol{b}$
($\boldsymbol{b}^{*}=\boldsymbol{b}$ versus
$\boldsymbol{b}^{*}\neq\boldsymbol{b}$)

  * •

define the stopping rules

Initialization: $\mathcal{S}=\\{g_{1},\mu_{1}\\},$ where
$g_{1}=\mathbb{R}^{p}$ and $\mu_{1}$ is the GLS estimator.

repeat

Given the predefined stopping rules, find the best node for splitting
($\tilde{s}$), the best covariate ($j_{\tilde{s}}$) and the best threshold
($c_{{\tilde{s}}}$) as follows:
$\displaystyle\tilde{s},j_{\tilde{s}},c_{{\tilde{s}}}$
$\displaystyle=\underset{s\in\mathcal{S},j\in
J_{s},c\in\mathbb{R}}{\mathrm{argmin}}\frac{1}{n}\sum_{i=1}^{n}Loss\big{(}y_{i},f\left(\boldsymbol{x}_{i}|j,c,\mathcal{S}/s\right)\big{)}$
$\displaystyle s.t:$
$\displaystyle\;\;\frac{1}{n}\sum_{i=1}^{n}Loss\big{(}y_{i},f\left(\boldsymbol{x}_{i}|j,c,\mathcal{S}/s\right)\big{)}<\frac{1}{n}\sum_{i=1}^{n}Loss\big{(}y_{i},f\left(\boldsymbol{x}_{i}|\mathcal{S}\right)\big{)},$
(2) where $J_{s}$ is the set of available covariates for splitting node $s.$

if $\\{\tilde{s},j_{\tilde{s}},c_{\tilde{s}}\\}$ exist then

Update $\mathcal{S}:$ replace $(g_{\tilde{s}},\mu_{\tilde{s}})$ by the new two
nodes, $\big{(}g_{\tilde{s}}\cap x_{j_{\tilde{s}}}\leq
c_{\tilde{s}},\mu_{\tilde{s}}^{r}(c_{\tilde{s}})\big{)},$
$\big{(}g_{\tilde{s}}\cap
x_{j_{\tilde{s}}}>c_{\tilde{s}},\mu_{\tilde{s}}^{l}(c_{\tilde{s}})\big{)},$
where $x_{j_{\tilde{s}}}$ is the covariate $j_{\tilde{s}}.$

end if

until $\\{\tilde{s},j_{\tilde{s}},c_{\tilde{s}}\\}$ do not exist or stopping
rules are satisfied $\forall s\in\mathcal{S}.$

As can be seen in Algorithm 1, RETCO covers various settings including
different prediction error measures and correlation structures. Before
analyzing its properties, here are some technical details for RETCO:

  * •

Predefined stopping rules: stopping rules, such as maximal tree’s depth and
minimal number of observations in a node, are commonly applied when fitting
regression trees (for more details, see Section 2.1).

  * •

Variance estimation: the variance components, needed for calculating
$f\left(\boldsymbol{x}_{i}|j,c,\mathcal{S}/s\right),$
$f\left(\boldsymbol{x}_{i}|\mathcal{S}\right)$ and the loss function, can be
estimated in different ways, for example using maximum likelihood or
restricted maximum likelihood. For clustered data, simple closed-form
equations that estimate the variance components are available and presented in
Appendix B.

  * •

Main optimization part:

    * –

the GLS predictors, $\mu_{s}^{l}(c),\;\mu_{s}^{r}(c)$ and
$\\{\mu_{l}\\}_{l\in\mathcal{S}/s}$ are estimated using dummy variables for
the current leaves: $\\{i|\boldsymbol{x}_{i}\in g_{s}\;\cap\;x_{i,j}\leq
c\\},$ $\\{i|\boldsymbol{x}_{i}\in g_{s}\;\cap\;x_{i,j}>c\\}$ and
$\\{i|\boldsymbol{x}_{i}\in g_{l}\\}_{l\in\mathcal{S}/s}$

    * –

for every potential split, the variance components are estimated using the
whole sample (rather than only $\\{y_{i}|\boldsymbol{x}_{i}\in g_{s}\\}$)

    * –

the condition of eq. (2) is required since the loss function is not a training
error (as in the standard regression tree algorithm) and therefore splitting a
node may increase the loss

  * •

Algorithm’s Output: since the tree estimates the marginal mean, its predictors
are $\\{\mu_{l}\\}_{l\in\mathcal{S}}$ rather than
$\\{\eta(\mu_{l})\\}_{l\in\mathcal{S}}.$

There are two main conceptual differences between RETCO and the standard
regression tree algorithm. The first is the use of a prediction error
estimator: Cp, AIC or $CV_{c}$ as loss functions. Unlike in the standard
regression tree, where the loss function is the training error, here the loss
function is prediction error estimator for correlated data. As was mentioned
in Section 2.1, although the regression tree model is a non-linear function of
$\boldsymbol{y},$ each split is a model selection problem of linear models.
Due to the linearity, Cp, AIC and $CV_{c}$ can be implement. More details
about the effect of using these prediction error estimators on the selected
thresholds and variables are given in Section 3.2. The second conceptual
difference is the iterative approach that is used here instead of the
recursive approach that is used in the standard regression tree model. The
iterative approach is expressed by selecting the optimal node, $\tilde{s},$
for splitting, rather than splitting each node independently. The reason for
using an iterative approach is that observations in different paths are
dependent. The dependency is accounted in the GLS predictors as well as in
other expressions in the loss function (e.g., bias correction and likelihood,
depending on the correlation setting and the selected loss function),
therefore splitting one node may affect the splitting of the other. This is in
contrast to i.i.d setting with training error loss function, where a recursive
approach can be used since observations in different paths are independent and
therefore splitting one node does not affect the splitting of the other.

For LMM-like setting, when $\boldsymbol{b}^{*}=\boldsymbol{b},$ once
$f(\cdot)$ was fitted, the random effects can be estimated using the BLUP
formula:

$\widehat{\boldsymbol{b}}=\widehat{\mathrm{Cov}}(\boldsymbol{y}^{*},\boldsymbol{y})\widehat{V}^{-1}\big{(}\boldsymbol{y}-f(X)\big{)},$

where $f(X)=[f(\boldsymbol{x}_{1}),...,f(\boldsymbol{x}_{n})].$

RETCO is formalized in context of LMM: $\\{\mu_{s}\\}_{s=1}^{S},$ are GLS
predictors and the variance is decomposed as in LMM,
$V=ZGZ^{t}+\sigma^{2}\times I_{n}$). However these properties are not
fundamental in the algorithm and can easily be generalized. For example, the
covariance matrices, $\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})$ and
$\mathrm{Var}(\boldsymbol{y}),$ can be expressed using a kernel covariance
function as is common in Gaussian process regression (and will be analyzed
numerically in Section 5). Also, other linear models instead of GLS can be
used for estimating the marginal means, $\\{\mu_{s}\\}_{s=1}^{S}.$ Moreover,
the random effects, $\boldsymbol{b},$ can be estimated in different ways than
using the BLUP formula.

Due to the iterative approach the complexity is higher than the complexity of
a standard regression tree. Assuming $x_{i,j}\neq x_{i^{\prime},j},\forall
i\neq i^{\prime}\in[1,...,n],\;j\in[1,...,p],$ the loss function is evaluated
less than $n\times p\times 2^{d}$ times, where $d$ is the tree’s depth.
Practically, due to hyper-parameters that restrict the number of potential
splitting values (e.g., such as minimum observations at each node) and due to
the constraint in inequality (2), which frequently shortens the depth of
paths, the number of evaluations is smaller. The complexity of each loss
function evaluation depends on the loss function, the predictor type and the
correlation structure. For example, in Cp loss function, GLS predictor and
clustered data, each evaluation complexity is $O(n_{c}^{3}),$ where $n_{c}$ is
the cluster size, and therefore the overall computational complexity is
$O\big{(}n\times p\times 2^{d}\times n_{c}^{3}\big{)}.$ The amount of memory
required is quadratic in $n,$ due to storing of
$\mathrm{Var}(\boldsymbol{y}).$

####  3.1.1 RETCO for RF and GB

Extending RETCO to RF and GB is done by averaging ensemble of implementations
of RETCO, while taking into account the special adjustments that RF and GB
require. The random effects are estimated in the same way. Several aspects
should be noted when implementing RF and GB:

  * •

Training set sampling: sampling with replacement cannot be implemented naively
when the loss function involves calculation of
$\mathrm{Var}(\boldsymbol{y})^{-1}$ (as when GLS predictor or marginal
likelihood loss are used) since in that case $\mathrm{Var}(\boldsymbol{y})$
might be a singular matrix due to duplication. Therefore, a half-sample method
should be used.

  * •

Number of trees: since the trees in RF and GB are correlated, then in order to
reduce the variance the number of trees should be large, especially when the
trees are deep. When correlated data are involved, the trees are even more
correlated due to the correlation between the observations. Therefore, the
number of trees should be even larger than in the i.i.d sample setting. RF
based on RETCO is demonstrated in Section 5.

  * •

Response in GB: in common regression GB the trees are fitted consecutively to
the residual of the previous tree. Therefore, the input of the algorithm in
all the trees except the first one is not $\boldsymbol{y}$. Correspondingly,
the estimated variance matrices relate to the residual of the previous tree,
rather than to $\boldsymbol{y}.$

####  3.1.2 Using CV Loss in Regression Tree

In typical predictive modeling settings, generalization error is the primary
objective of learning, hence CV loss is the natural choice. Surprisingly,
there is not much previous work on using CV loss in trees, even without
correlation. Notable exceptions are the ALOOF algorithm (Painsky and Rosset,
2016) and approximations used in CatBoost (Prokhorenkova et al., 2017). The
main drawback in using CV-based loss function is increasing the computational
cost compares to Cp loss function. Moreover, since RETCO is iterative rather
than recursive, the number of evaluations of the loss function can remain
large for all the splits along the tree.

###  3.2 The Bias Correction Effect

As explained in the previous sections, we suggest to add a bias correction
term to the training error such that the loss function estimates the
prediction error. This section illustrates the effect of the bias correction
on split selection. Extensive numerical analysis is presented in Section 5.

####  3.2.1 $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ Scenario

Observations with positive correlation are similar in higher probability than
uncorrelated observations. Therefore, loss functions that do not take into
account the correlation, tend to split a node based on the correlation
structure of the training set observations rather than their mean. Splitting
based on the correlation of the training data is not useful for predicting
uncorrelated observations. Therefore, when
$\boldsymbol{b}^{*}\perp\boldsymbol{b}$ it is important to fit the regression
tree based on the marginal mean only. The corrections in Cp and $CV_{c},$
which take into account the correlation structure, balance this tendency.
Examples 3.1 and 3.2 demonstrate this mechanism. The code for the examples, as
well as for the numerical part in Section 5, is written in Python and is
available in https://github.com/AssafRab/RETCO.

###### Example 3.1

Consider the setting of $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ and a training
data containing four observations from two clusters with the covariance matrix
$\small\mathrm{Var}(\boldsymbol{y})=\scriptsize\begin{pmatrix}2&1&0&0\\\
1&2&0&0\\\ 0&0&2&1\\\ 0&0&1&2\\\ \end{pmatrix},$ i.e., observations 1 and 2
belong to the first cluster and observations 3 and 4 belong to the second
cluster. The $CV_{c}$ correction is reduced in this setting to
$2\mathrm{tr}\big{(}H_{cv}\mathrm{Var}(\boldsymbol{y})\big{)}/n.$ Two models
with the GLS predictor are tested:

  * •

Model A, which splits the training set into the two clusters. Given the
covariance matrix:

$\small H_{cv}=\begin{pmatrix}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{pmatrix}\to\frac{2}{n}\mathrm{tr}\big{(}H_{cv}\mathrm{Var}(\boldsymbol{y})\big{)}=2.$

  * •

Model B, that mixes between the clusters and selects observations 1 and 3 for
one subset and 2 and 4 for the other subset. Given the covariance matrix:

$\small H_{cv}=\begin{pmatrix}0&.25&1&-.25\\\ .25&0&-.25&1\\\ 1&-.25&0&.25\\\
-.25&1&0.25&0\end{pmatrix}\to\frac{2}{n}\mathrm{tr}\big{(}H_{cv}\mathrm{Var}(\boldsymbol{y})\big{)}=0.5.$

As we can see in Example 3.1, decomposing the penalty,
$2\mathrm{tr}\big{(}H_{cv}\mathrm{Var}(\boldsymbol{y})\big{)}/n,$ shows that
the weights that relate to observations from the same cluster are multiplied
by their positive covariance values and therefore contribute to the penalty,
while weights that relate to observations from different clusters are
multiplied by zero and therefore do not contribute to the penalty. As a result
the penalty of model A, which gives the whole weight for observations from the
same cluster, is larger than for model B. Therefore, while $CV$ selects model
A when $CV(A)<CV(B),$ $CV_{c}$ selects model A when $CV(A)<(CV(B)-1.5).$ In
that way $CV_{c},$ as well as $Cp,$ balance the tendency to split based on
based on the correlation structure of the training set. Obviously, as much
much as the observations in the training set are more correlated, the penalty
effect is stronger, and the superiority of RETCO over the standard algorithm
is more prominent (see also Section 5).

###### Example 3.2

Consider the setting of $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ and
$\boldsymbol{y}\sim N_{100}(0.1\times\boldsymbol{x}_{1},V),$ where

$\small\mathrm{Var}[i,j]=\begin{cases}2,&\text{when }i=j\\\ 1,&\text{when
}i\neq j\text{ and }(i,j\leq 50\;\text{ or }i,j>50)\\\ 0,&\text{ o.w
}\end{cases},$

and

$\displaystyle\small\boldsymbol{x}_{1}[i]$
$\displaystyle=\begin{cases}0.5+\epsilon_{x},&\text{when i is
odd},\quad\epsilon_{x}\sim N(0,0.1)\\\ -0.5+\epsilon_{x},&\text{when i is
even},\end{cases}$

i.e., $\boldsymbol{y}$ contains two clusters of $50$ observations each, and
its mean is not correlated with the clusters.

Two models are tested, model A which uses the threshold $\boldsymbol{x}_{1}=0$
and model B which uses $\boldsymbol{x}_{2}=0,$ where $\boldsymbol{x}_{2}$ is
highly correlated with the clusters:

$\displaystyle\boldsymbol{x}_{2}[i]$
$\displaystyle=\begin{cases}0.5+\epsilon_{x},&\text{when }i\leq
50,\quad\epsilon_{x}\sim N(0,0.1).\\\ -0.5+\epsilon_{x},&\text{when
}50<i.\end{cases}$

A simulation of this setting is visualized in Figure 1. In this simulation
$CV(A)=1.47$ and $CV(B)=1.04,$ while $CV_{c}(A)=2.47$ and $CV_{c}(B)=3.04.$
Therefore, in case $\boldsymbol{b}^{*}\perp\boldsymbol{b},$ while CV selects
model B, $CV_{c}$ selects model A. 111For the setting
$\boldsymbol{b}^{*}=\boldsymbol{b},$ which will be discussed next, $CV_{c}=CV$
and both select model B.

Figure 1: The dashed blue and green lines are the thresholds for model A and
model B respectively.

Unlike Cp and $CV_{c},$ whose penalties depend on
$\mathrm{Var}(\boldsymbol{y}),$ the penalty of AIC in this setting
($\boldsymbol{b}^{*}\perp\boldsymbol{b}$) is fixed regardless of the training
covariance structure at $2p/n.$ Therefore, we can conclude that in AIC the
likelihood, $\ell(\boldsymbol{y};\widehat{\mathbb{E}}(\boldsymbol{y}|X),V),$
is responsible for mixing uncorrelated training set observations in the
different paths, while the penalty only affects the stopping rule.

####  3.2.2 $\boldsymbol{b}^{*}=\boldsymbol{b}$

When $\boldsymbol{b}^{*}=\boldsymbol{b}$ the correlation between
$\boldsymbol{y}^{*}$ and $\boldsymbol{y}$ is the same as the correlation
between observations in $\boldsymbol{y}.$ Therefore, unlike in the
$\boldsymbol{b}^{*}\perp\boldsymbol{b}$ setting, here there is no clear
motivation to restrict the tendency to split the nodes based on the
correlation structure of the training set (as appears in standard regression
trees). Correspondingly, the bias corrections in this setting are also
different than in the $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ setting. For
example, $CV$ is not biased in this setting, i.e., $CV_{c}=CV$ (for more
details see Section 2.3). The penalty in Cp, $2\sigma^{2}\mathrm{tr}(H)/n,$
depends on $\mathrm{Var}(\boldsymbol{y})$ through $H$ for some models (e.g.,
for LMM), however for other models it does not depend implicitly on
$\mathrm{Var}(\boldsymbol{y}).$ In any case, the effect of
$\mathrm{Var}(\boldsymbol{y})$ is much less prominent than in
$\boldsymbol{b}^{*}\perp\boldsymbol{b}$ setting, where the bias is
$2\mathrm{tr}\big{(}H\mathrm{Var}(\boldsymbol{y})\big{)}/n.$ Similarly with
the penalty in AIC, $2\mathrm{tr}(H)/n.$ Therefore, in this setting, both Cp
and AIC penalties mainly affect to the stopping rule rather than mix between
uncorrelated observations. Given additional stopping rules (e.g., tree depth,
minimal number of training set observations in each node), we can conclude
that the effect of using prediction error estimator instead of training error
is limited. A numerical analysis of this scenario is presented in Section
5.1.2. Still, it is important to emphasize that the proposed
model—$f(\boldsymbol{x}^{*})+\boldsymbol{b}^{*^{t}}\boldsymbol{z}^{*}$—is
recommended also because of the inference and the use of a solid statistical
perspective that the random effects framework enables.

####  3.2.3 $\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$ but
$\boldsymbol{b}^{*}\neq\boldsymbol{b}$ Scenario

From a qualitative perspective, this scenario is the same as the
$\boldsymbol{b}^{*}\perp\boldsymbol{b}$ scenario. In both scenarios, the
correlation structure of $\boldsymbol{y}$ is not preserved in the prediction
problem. As a result, the bias correction has a key role in balancing the
tendency of standard regression trees to split based on the correlation
structure of $\boldsymbol{y}.$ The main difference between the scenarios is
quantitative and is explicitly expressed in the bias corrections formulas, for
example the bias correction in Cp is
$2\mathrm{tr}\Big{(}H\big{(}\mathrm{Var}(\boldsymbol{y})-\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}^{*})\big{)}\Big{)}/n.$
This setting of $\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$ but
$\boldsymbol{b}^{*}\neq\boldsymbol{b}$ is demonstrated in Section 5.2.

##  4 Comparison With Other Algorithms

To our knowledge Sela and Simonoff (2012) were the first to propose
integration between LMM and regression tree by introducing the RE-EM
algorithm. The main idea in RE-EM is generating
$\\{g_{s}\\}_{s=1}^{\mathcal{S}}$ using a standard regression tree that is
fitted to the residual, $\boldsymbol{y}-Z\widehat{\boldsymbol{b}}.$ Given
$\\{g_{s}\\}_{s=1}^{\mathcal{S}},$ $\\{\mu_{s}\\}_{s\in\mathcal{S}}$ are
estimated by GLS. Therefore, the correlation is taken into account in
estimating $\\{\mu_{s}\\}_{s\in\mathcal{S}},$ but it is ignored in selecting
$\\{g_{s}\\}_{s=1}^{\mathcal{S}}.$ RE-EM algorithm is presented in Appendix C.
Hajjem et al. (2014) proposed a RF algorithm which is based on the same logic
as RE-EM algorithm. For numerical comparison between RE-EM and RETCO, see
Section 5.1.4.

Stephan et al. (2015) proposed the Mixed Random Forest (MRF), that does not
ignore the correlation structure when selecting $\\{g_{s}\\}_{s=1}^{S},$
however still does not address the correlation correctly. The goal in MRF
(which is also presented in Appendix C) is fitting a model that estimates
accurately the variance components, rather than optimizing prediction
accuracy, as in RETCO. Also, MRF assumes a specific data type and is based on
strong distributional assumptions ($\boldsymbol{y}$ is normally distributed
and $G=\sigma^{2}_{b}I_{q}$). Besides the difference in goals and the assumed
settings, the main difference in the tree fitting approach is that MRF uses
likelihood loss function for finding $\\{g_{s},\mu_{s}\\}_{s=1}^{S},$ rather
than estimated prediction error for correlated data (i.e., training error plus
a bias correction) as in RETCO. Also, MRF is based on a standard recursive
approach, while RETCO is based on an iterative approach (for the motivation of
using an iterative approach see Section 3). For numerical comparison between
MRF and RETCO, see Section 5.1.4.

Extensions of these papers, where the response is binary or count data, as in
generalized linear mixed model (Wolfinger and O’connell, 1993), where proposed
by Fokkema et al. (2018); Hengl et al. (2018); Ngufor et al. (2019); Speiser
et al. (2019).

##  5 Numerical Results

This section compares the performance of RETCO with the standard regression
tree algorithm and relevant modifications of it that will be described. The
analysis is performed using both simulated data and real data sets for
different correlation settings. The simulation part is based on random effects
framework and presents results for $\boldsymbol{b}^{*}=\boldsymbol{b}$ as well
as for $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ correlation settings. The
$\boldsymbol{b}^{*}\neq\boldsymbol{b}$ but
$\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$ correlation setting is analyzed in
the random field context using a real data set with spatial correlation. Also,
different prediction error estimator types (Cp, $CV_{c}$ and AIC) and
different trees-based models (regression tree and RF) are analyzed. The code
is available in https://github.com/AssafRab/RETCO.

###  5.1 Simulation

The training set was generated from the following model:

$\displaystyle\boldsymbol{y}=$ $\displaystyle
I_{(\boldsymbol{x}_{1}>0)}+I_{(\boldsymbol{x}_{2}>0)}+I_{(\boldsymbol{x}_{3}>0)}+I_{(\boldsymbol{x}_{1}>0)}I_{(\boldsymbol{x}_{2}>0)}I_{(\boldsymbol{x}_{3}>0)}+Z\boldsymbol{b}+\boldsymbol{\epsilon},$

where

  * •

$I_{(\boldsymbol{x}_{j}>0)},\;\forall j\in[1,2,3]$ is the indicator vector for
$\;x_{j,i}>0,\;\forall i\in[1,...,n].$

  * •

The sample contains $C$ clusters, each one of size $n_{c}=n/C.$
$Z\in\mathbb{R}^{n\times C}$ indicates the clusters, i.e., for the first
column the first $n_{c}$ elements are $1,$ and the rest are zero, for the last
column the last $n_{c}$ elements are $1$ and the rest are zero.

  * •

$\boldsymbol{b}\in\mathbb{R}^{C}$ is the random effects vector, distributed
$N_{C}(0,\sigma^{2}_{b}I_{C}),$ and $\boldsymbol{\epsilon}\sim
N_{n}(0,I_{n}).$

  * •

$\boldsymbol{x}_{1},\;\boldsymbol{x}_{2}$ and $\boldsymbol{x}_{3}$ are
$Z\boldsymbol{\gamma}+\boldsymbol{\eta},$ where $\boldsymbol{\gamma}\sim
N_{C}(0,\sigma^{2}_{b}I_{C})$ and $\eta_{i}$ are uncorrelated and distributed
uniformly, $U(-1,1),\forall i\in[1,...,n].$

####  5.1.1 New Random Effects ($\boldsymbol{b}^{*}\perp\boldsymbol{b}$)

As was mentioned in Section 2.2, when $\boldsymbol{b}^{*}\perp\boldsymbol{b},$
i.e., when $\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})=0,$ GLS estimator
should be used instead of LMM. Here, Cp prediction error estimator is analyzed
and therefore $\boldsymbol{y}^{*}$ should be related to the same covariate
values as in the training set,
$\Phi=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3},Z\\}.$ In
order to reduce the variance of the prediction error estimate, the test sample
contains $300$ replicates of $\Phi.$ RETCO is compared to a standard
regression tree with squared error loss function, which is the same loss as in
Cp but without the bias correction term. In both algorithms, the stopping
rules are depth of tree smaller than $4,$ and number of observations in the
terminal node greater than $2.$ The relative difference between the RETCO test
error and its alternative:

$\text{error difference[\%]}=\frac{\text{error(RETCO)}-\text{error(standard
tree)}}{\text{error(standard tree)}},$

is calculated repeatedly for $100$ simulation runs. The average simulation run
time is $82.1$ seconds, where RETCO takes on the order of 3-8 fold longer to
run due to its iterative approach. For more details about RETCO’s
computational complexity, see Section 3.1.

Figure 2, left panel, presents boxplots of error difference[%] for different
$\sigma_{b}^{2}.$ As expected, when $\sigma_{b}^{2}$ is larger (i.e., the
correlation is stronger), the improvement in using RETCO over the standard
regression tree algorithm is bigger. However, also for relatively small
$\sigma^{2}_{b}$ values RETCO outperforms the standard regression tree
algorithm. Additional comparisons for different sample sizes ($n$) and cluster
sizes ($n_{c}$), and similar analyses for generalization error setting using
$CV_{c},$ are given in Appendix D.

Figure 2: Boxplots of the error difference[%] for $n=300,\;n_{c}=50$ and
different $\sigma^{2}_{b}$ $(0.5,\;1,\;5)$ for different settings. Left:
$\boldsymbol{b}^{*}\perp\boldsymbol{b}$ setting. The means are $-2.5,\;-9.9$
and $-25.2$. Middle: $\boldsymbol{b}^{*}=\boldsymbol{b}$ setting. The means
are $-2.1,\;-1.1$ and $-10.4.$ Right: RF setting. Two versions of RETCO are
analyzed, in light the version that does not enforce the stopping rule
constraint, and in bold the version that enforces the constraint. The means
are $-2.1,\;-2.9$ $-11.7$ and $-2.2,\;-9.6$ and $-34.7,$ respectively.

As was mentioned in Section 3.2, RETCO balances the tendency of standard tree
to split based on the correlation structure of the training sample. Therefore,
we expect that RETCO mixes training set observations from different clusters
in the leaves more than the standard regression tree. The following measure
quantifies this mixing property:

$\small\text{homogeneity}=\sum_{s\in\mathcal{S}}\sum_{c=1}^{C}|n(c,s)-n(s)/C|,$

where $n(c,s)$ is the number of training set observations in leaf $s$ that
belong to cluster $c$ and $n(s)=\sum_{c=1}^{C}n(c,s).$ Smaller homogeneity
means bigger mixing. Figure 3 plots the training sample homogeneity
difference[%],

$\small\frac{\text{homogeneity(RETCO)}-\text{homogeneity(standard
tree)}}{\text{homogeneity(standard tree)}},$

versus the error difference[%] for different $\sigma^{2}_{b}$ values.

Figure 3: Homogeneity difference[%] versus error difference[%] for
$n=300,\;n_{c}=50$ and different $\sigma_{b}^{2}$ values.

As can be seen in Figure 3, error difference[%] has a positive correlation
with homogeneity difference[%], i.e., the property of RETCO to balance the
tendency of standard regression tree to split based on the correlation
structure of the training sample is essential.

####  5.1.2 Same Random Effects ($\boldsymbol{b}^{*}=\boldsymbol{b}$)

For the $\boldsymbol{b}^{*}=\boldsymbol{b}$ scenario, the training set is the
same as in Section 5.1.1, but the prediction set is different, such that the
random effects realizations from the training set are also used for
constructing $\boldsymbol{y}^{*}.$

For RETCO, AIC loss function with LMM predictor for $\boldsymbol{y}$ is used
for splitting. As presented in Section 3, the predictor in the tree’s leaves
is GLS and the random effects term is added after fitting the tree. For the
standard regression tree, normal likelihood loss function is used with no
distinction between random and fixed effects, i.e., all the covariates,
including the cluster, can be selected for splitting. Correspondingly, the
minus log likelihood of i.i.d normal distribution (which is effectively the
same as squared error loss) is used as a loss function for the standard
regression tree algorithm. The middle panel of Figure 2 compares between the
algorithms for different $\sigma_{b}^{2}$ values (when $n=300$ and
$n_{c}=50$). As we can seen, RETCO outperforms the standard tree algorithm.
For $\sigma_{b}^{2}=0.5$ and $\sigma_{b}^{2}=1,$ the average error
difference[%] is relatively small. As was mentioned in Section 3.2.2, this
phenomenon is expected.

####  5.1.3 RF - New Random Effects

RF is analyzed for $\boldsymbol{b}^{*}\perp\boldsymbol{b}$ and Cp loss
function setting. The training sample model is:

$\displaystyle\boldsymbol{y}=$ $\displaystyle
I_{(\boldsymbol{x}_{1}>0)}+I_{(\boldsymbol{x}_{2}>0)}+I_{(\boldsymbol{x}_{3}>0)}+I_{(\boldsymbol{x}_{1}>0)}I_{(\boldsymbol{x}_{2}>0)}I_{(\boldsymbol{x}_{3}>0)}+I_{(\boldsymbol{x}_{4}>0)}I_{(\boldsymbol{x}_{5}>0)}I_{(\boldsymbol{x}_{6}>0)}+Z\boldsymbol{b}+\boldsymbol{\epsilon},$

where $\boldsymbol{x}_{i}\;\forall i\in\\{1,...,6\\},\;Z\boldsymbol{b}$ and
$\boldsymbol{\epsilon}$ have the same distribution as in Section 5.1.1.
Additional parameters that are relevant for RF are:

  * •

the maximal tree depth is $10$

  * •

The number of regression trees is $T=100$

  * •

A random half-sample method is used for sampling the training set for each
tree (i.e., the training sample size for each tree is 250 without duplicates)

  * •

Three covariates are randomly selected at each split, following the rule of
thumb of selecting randomly $round\big{(}\log_{2}(p)\big{)}$ potential
covariates at each split.

Also, two versions of RETCO are tested: the first uses the stopping rule
constraint as presented in RETCO inequality (2), the second does not enforce
the constraint, and therefore results in deeper trees. The prediction set
contains new random effects realizations, such that
$\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})=0.$ The covariates of the
prediction set are $200$ replicates of the covariates of the training set. The
analysis was repeated $50$ times. The right panel in Figure 2 presents
boxplots of the error difference[%] for different $\sigma^{2}_{b}.$ As we can
seen, both versions of RETCO outperform the standard algorithm. Also, forcing
the stopping rule gives better results.

####  5.1.4 Comparison With Previous Algorithms

The competitors in the left and the right graphs in Figure 2 preserve the main
characteristics of the MRF algorithm: taking into account the correlations
structure by using GLS estimator and differentiation between random and fixed
effects. The exact MRF algorithm, which was designed for a genetic
application, is not implemented here since some of its technical details are
specific for genetic applications, which are not our main use case. A
comparison between RETCO and RE-EM algorithm is presented in Appendix D. As
expected, RETCO’s performance is uniformly superior to both algorithms due to
its use of prediction error estimates for splitting and the careful
consideration of correlation structures in splitting and prediction.

###  5.2 Real Data Analysis

This section presents real data analyses comparing the performance of the
standard regression tree and RF algorithms to their RETCO versions for six
different data sets with various correlation structures. The data sets and the
prediction problems are briefly described in Section 5.2.1, additional
technical information can be found in Appendix D.

Table 2 summarizes the results. As can be seen, for all the six analyses the
test errors of the standard regression tree and RF algorithms are greater than
the test errors of their RETCO versions (negative error difference [%]),
moreover in several analyses the improvement of RETCO over the standard
algorithm is very large.

Correlation Structure | Data Set Name | Regression Tree | RF
---|---|---|---
Clusters | FIFA | $-4.9\%$ | $-7.1\%$
Crimes | $-8.1\%$ | $-4.4\%$
Spatial | Korea Temperature | $-14.5\%$ | $-13.6\%$
California Housing | $-14.2\%$ | $-3.3\%$
Longitudinal | Parkinson’s Disease | $-32.9\%$ | $-35.4\%$
Wages | $-20.1\%$ | $-7.9\%$
Table 2: Error difference [%] between RETCO and standard regression tree, and
between RETCO and standard RF.

####  5.2.1 Prediction Problems Description

  * •

FIFA – This data set contains football players’ market-values. The data set
has a clustered correlation structure, where the cluster variable is the
player’s club, such that market-values of players from the same club are
correlated but from different clubs are not correlated. The prediction goal is
to predict the market-values of new players from new clubs. In order to
satisfy this prediction goal, the training set contains the observations of
players from $20$ clubs that were randomly sampled ($548$ observations), and
the test set contains the observations of the other clubs ($17,939$
observations). Since the covariate values of the prediction set are not the
same as the covariate values of the training set, CV-type loss function is
used in the algorithms’ splitting criterion – $CV_{c}$ for RETCO and CV for
the standard regression tree/RF algorithms. The prediction error is estimated
by the average squared errors of the test set. The data set is publicly
available on Kaggle.

  * •

Communities and Crime – This data set presents the number of violent crimes
per population size in US communities. The data set has a clustered
correlation structure, where the clusters are the US states (each state
contains many communities). The training set contains $15$ clusters that were
randomly sampled ($790$ observations), where the test set contains the other
clusters ($1,204$ observations). For the same reason as in the FIFA data set,
CV-type loss function is used in the algorithms’ splitting criterion. The data
set was introduced by Redmond and Baveja (2002), and is publicly available on
the UCI repository.

  * •

South Korea Temperature (’bias correction of numerical prediction model
temperature forecast’) – This data set contains daily maximal temperature
measurements (collected in August between the years $2013-2017$) at several
sites in South Korea. The prediction goal is to predict the maximal
temperature of new days. Measurements of the first two years were selected
($575$ observations) in order to predict the maximal temperature of the same
set of days in the next years ($2,325$ observations). Due to the spatial
correlation structure, exponential kernel covariance function was used for
modeling. Since all the records are measured at the same sites, in-sample
error type is used in the algorithms’ splitting criterion ($Cp$ for RETCO and
squared error loss for the standard regression tree/RF). The data set was
introduced by Cho et al. (2020), and is publicly available on UCI repository.

  * •

California Housing Prices – This data set contains the median house value
within a block for different blocks in California. Some of the blocks belong
to the same clusters (same coordinate values), therefore the data set has a
clustered-spatial correlation structure, which can be represented by the
following kernel covariance function:

$\displaystyle\mathrm{Cov}(y_{i},y_{j})=$
$\displaystyle\mathcal{K}(|\boldsymbol{z}_{i}-\boldsymbol{z}_{j}|)+\sigma_{b}^{2}I_{cluster(i)=cluster(j)}+\sigma_{\epsilon}^{2}I_{i=j},$

where $\mathcal{K}(\cdot)$ is the exponential kernel covariance function,
$\boldsymbol{z}_{i},\boldsymbol{z}_{j}$ are the coordinates and $cluster(i),$
$cluster(j)$ are the clusters of $y_{i},y_{j}.$ The prediction goal is to
predict the median house value of new blocks from new clusters. Therefore,
$100$ clusters are randomly sampled ($279$ observations) for the training set
and the other clusters are used as the test set ($12,124$ observations). Since
the prediction goal is to predict median house values from new clusters, then
$\boldsymbol{b}^{*}\neq\boldsymbol{b}.$ However, due to the spatial
correlation, the prediction set median house values are correlated with the
training set median house values, and therefore this setting satisfies the
$\boldsymbol{b}^{*}\neq\boldsymbol{b}\cap\boldsymbol{b}^{*}\not\perp\boldsymbol{b}$
scenario (see Section 2.2.1). CV-type loss function is used in the algorithms’
splitting criterion since the covariate points of the training set and the
prediction set are different. The data set was introduced by Pace and Barry
(1997), and is publicly available on Kaggle.

  * •

Parkinson’s Disease Telemonitoring – This longitudinal data set contains
Parkinson’s disease symptom scores of $42$ individuals along six-months trial.
The clustered-temporal correlation structure, where the clusters refer to the
individuals, can be modeled by LMM with random intercept for the cluster and
random slope for the time variable. The prediction goal is predicting the
score of new individuals, therefore five individuals were randomly sampled
($742$ observations) for the training set and the others were designated as
the test set ($5,179$ observations). The covariates in this data set are
biomedical voice measurements and their values are approximately the same for
all the individuals. Therefore, $Cp$ and squared error loss loss functions are
used in the splitting criterion for RETCO and standard regression tree/RF,
respectively. The data set was introduced by Tsanas et al. (2009), and is
publicly available on the UCI repository.

  * •

Wages – This longitudinal data set presents the average hourly wages by year
of $888$ employees. As in the Parkinson’s Disease Telemonitoring data set,
this data set can be modeled by LMM with random intercept and random slope,
where the employee is the cluster variable. The prediction goal is to predict
the average hourly wage of new employees. $50$ individuals were randomly
sampled ($331$ observations) for the training set and the other are used as
the test set ($6,071$ observations). Since the covariate values of the
training set and the prediction set are different, CV-type loss function is
used in the algorithms’ splitting criterion. The data set was introduced by
Singer et al. (2003), and is publicly available in brolgar package in R
software.

##  6 Conclusions

This paper presents a new algorithm, RETCO, for fitting regression trees-based
models for correlated data. Analyzing various settings with different
correlation structures lead to the conclusion that RETCO substantially
improves prediction performance in settings involving correlated data.

Unlike standard regression trees-based models, which ignore the correlation
structure of the data, RETCO accounts for the correlation structure in various
ways, such as using prediction error estimates for correlated data as the loss
function in the splitting criterion. As discussed and demonstrated, using
prediction error estimators for correlated data instead of training error
neutralizes the tendency to fit a tree that divides the training set based on
its correlation structure, as is likely to happen in standard regression
trees-based models.

Extensive data analysis, including analysis of six different real data sets,
shows the superiority of RETCO over standard regression trees-based model, as
well as its generality that enables to implement it under various settings.

Acknowledgments

This work was supported by the Israel Science Foundation, grant 1804/16 and by
the European Union Seventh Framework Programme grant agreement no. 785907
(Human Brain Project).

##  A Theoretical Background

This appendix extends the theoretical background that is given in Section 2.

###  A.1 Regression Tree, Random Forest and Gradient Boosting

Algorithm 2 presents a typical regression tree algorithm.

Algorithm 2 A typical tree-based algorithm

Input: $\boldsymbol{y},\;X.$

Output: $f(\cdot).$

High-level setting: select a training error loss function and define stopping
rules

Initialization: $\mathcal{S}=\\{g_{1},\mu_{1}\\},$ where
$g_{1}=\mathbb{R}^{p},\;\mu_{1}=\sum_{i=1}^{n}y_{i}/n.$

repeat

  1. 1.

Given the predefined stopping rules, for each node $s\in\mathcal{S}$ solve the
following optimization problem:

$\displaystyle c_{s},j_{s}=\underset{c\in\mathbb{R},j\in
J_{s}}{\mathrm{argmin}}\frac{1}{|I_{s}|}\sum_{i\in
I_{s}}Loss\big{(}y_{i},I_{(x_{i,j}\leq
c)}\mu_{s}^{l}(c)+I_{(x_{i,j}>c)}\mu_{s}^{r}(c)\big{)},$

where $J_{s}$ is the set of available covariates for splitting node $s,$
$I_{s}=\\{i|\boldsymbol{x}_{i}\in g_{s}\\},$ $\mu_{s}^{l}(c)$ and
$\mu_{s}^{r}(c)$ are the mean estimators of $\\{y_{i}|\boldsymbol{x}_{i}\in
g_{s}\cap x_{i,j}\leq c\\}$ and $\\{y_{i}|\boldsymbol{x}_{i}\in g_{s}\cap
x_{i,j}>c\\}$ respectively.

  2. 2.

Update $\mathcal{S}$ by replacing $(g_{s},\mu_{s})$ by the new two nodes:
$\big{(}g_{s}\cap x_{j_{s}}\leq c_{s},\mu_{s}^{r}(c_{s})\big{)},$
$\big{(}g_{s}\cap x_{j_{s}}>c_{s},\mu_{s}^{l}(c_{s})\big{)},$ where
$x_{j_{s}}$ is the covariate $j_{s}$ and
$\mu_{s}^{r}(c_{s}),\mu_{s}^{l}(c_{s})$ are the related mean predictors.

until Stopping rules are satisfied $\forall s\in\mathcal{S}$

Random forest (RF) and gradient boosting (GB) predictor are based on averaging
an ensemble of trees:

$f(\boldsymbol{x}^{*})=\sum_{t=1}^{T}\lambda_{t}\sum_{s=1}^{S_{t}}I_{(\boldsymbol{x}^{*}\in
g_{t,s})}\mu_{t,s},$

where $T$ is the number of trees and $\lambda_{t}\in(0,1]$ is the learning
rate (for RF $\lambda_{t}=1/T,\;\forall t$).

The regression trees in RF and GB are fitted in different ways than in a
standard regression tree. In RF, the training set of each tree is sampled from
the original sample (e.g., sampling with replacement of size n, half-sample),
and the set of the potential covariates of each split is a random sample of
$J_{s}.$ In GB the trees are dependent and created consecutively, where the
dependent variable of each tree is the residual of the previous tree. Also,
there are many techniques for reducing over-fitting and model variance which
are relevant for RF and GB, but not relevant for standard regression tree
model. For more information about RF and GB see Freund et al. (1999); Friedman
(2001); Breiman (2001); Hastie et al. (2009).

Note, unlike in regression tree model, which tends to over-fit and therefore
suffers from high variance, RF has relatively low variance due to the
averaging over the $T$ trees. This property affects the optimal structure of
trees in RF. While the tree depth in regression tree model should be
restricted in order to avoid over-fitting, the trees in RF can be large
whenever $T$ is respectively large (Criminisi et al., 2011). Since the trees
in RF are correlated, the RF variance decreases in a smaller rate than $T.$
Commonly the trees’ depth in RF is also restricted for various reasons, such
as computational cost that RF with deep trees (and consequently large $T$)
requires.

###  A.2 Prediction Error Estimation and Model Selection for Correlated Data

####  A.2.1 Cp

The original Cp, when
$\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})=\sigma^{2}\times I_{n}$ and
$\mathrm{Cov}(\boldsymbol{y}^{*},\boldsymbol{y})=0,$ was introduced by Mallows
(1973) is:

$\displaystyle
Cp=\frac{1}{n}\|\boldsymbol{y}-\widehat{\boldsymbol{y}}\|_{2}^{2}+\frac{2\sigma^{2}}{n}p.$

####  A.2.2 AIC

The standard AIC under normality and i.i.d assumptions, that was introduced by
Akaike (1974) is:

$AIC=-\frac{2\ell\big{(}\boldsymbol{y};\widehat{\mathbb{E}}(\boldsymbol{y}|X),\sigma^{2}\big{)}}{n}+\frac{2p}{n},$

where
$\ell\big{(}\boldsymbol{y};\widehat{\mathbb{E}}(\boldsymbol{y}|X),\sigma^{2}\big{)}$
is the log-likelihood of $\boldsymbol{y}.$

##  B Estimating Variance Components for Clustered Data

When $\boldsymbol{y}$ has a clustered correlation structure, i.e., its
covariance matrix follows:

$\mathrm{Cov}[i,j]=\begin{cases}\sigma^{2}_{\epsilon}+\sigma^{2}_{b},&\text{when
}i=j\\\ \sigma^{2}_{b},&\text{when }i\neq j\text{ but }c(i)=c(j)\\\ 0,&\text{
o.w }\end{cases},$

where $c(i)$ is the cluster that observation $i$ belongs to and
$\sigma_{\epsilon},\;\sigma_{b}$ are in $R^{+},$ then $\sigma^{2}_{\epsilon}$
and $\sigma^{2}_{b}$ can be estimated in a closed-form way. In order to
simplify the equations let us assume $\mathbb{E}\boldsymbol{y}=0.$ In this
case:

$\widehat{\sigma}_{\epsilon}^{2}=\frac{\sum_{i=1}^{n}\big{(}y_{i}-\bar{y}(i)\big{)}^{2}}{(n-C)},$

where $\bar{y}(i)$ is the average of the cluster that $y_{i}$ belongs to, and
$C$ is the number of clusters,

$\widehat{\sigma}_{b}^{2}=\Big{(}\frac{\sum_{i=1}^{n}(\bar{y}\big{(}i)-\bar{y}\big{)}^{2}}{C-1}-\widehat{\sigma}_{\epsilon}^{2}\Big{)}\times\frac{C-1}{n-\sum_{j=1}^{C}n_{j}^{2}/n},$

where $\bar{y}$ is the average of $\boldsymbol{y}$ and $n_{j}$ is the number
of observations in cluster $j.$ When $\mathbb{E}(\boldsymbol{y})\neq 0,$ the
variance parameters should be calculated for the residual,
$\boldsymbol{y}-\widehat{E}{\boldsymbol{y}}.$

##  C Comparison with Other Methods

Algorithm 3 presents the RE-EM algorithm, which was proposed by Sela and
Simonoff (2012). Algorithm 4 presents MRF algorithm for a single tree, that
was proposed by Stephan et al. (2015).222Stephan et al. (2015) do not supply
an organized algorithm, Algorithm 4 tries to formalize their approach as given
in their supplementary material. For MRF, in order to simplify notations,
denote $\boldsymbol{y},\;X,\;Z$ as the bootstrap sample and ignore the
features sampling at each split.

Algorithm 3 RE-EM Algorithm

Input: $\boldsymbol{y},\;X,\;Z.$

Output: $f(\cdot),\;\widehat{\boldsymbol{b}}.$

Initialization: $\widehat{\boldsymbol{b}}^{(0)}=0.$

repeat

  1. 1.

Using the fixed effects covariates, fit CART algorithm (Breiman et al., 1984)
for $(y-Z\widehat{\boldsymbol{b}}^{(k-1)})$ and extract
$\\{g^{(k)}_{s}\\}_{s=1}^{S^{(k)}}$ from the fitted tree.

  2. 2.

Estimate $\\{\mu_{s}\\}_{s=1}^{S},$ $G$ and $\sigma^{2}$ using the following
LMM model:

$y_{i}=\sum_{s=1}^{S}I_{(\boldsymbol{x}_{i}\in
g_{s}^{(k)})}\mu_{s}+\boldsymbol{z}_{i}\boldsymbol{b}+\epsilon_{i},$

where $\boldsymbol{b}\sim N_{q}(0,G),\;\epsilon_{i}\sim N(0,\sigma^{2}).$

  3. 3.

Given $\\{\mu_{s}^{(k)}\\}_{s=1}^{S^{(k)}},\;\widehat{G}^{(k)}$ and
$\widehat{V}^{(k)}=Z\widehat{G}^{(k)}Z^{t}+\widehat{\sigma}^{2,(k)}I_{n}:$
estimate $\widehat{\boldsymbol{b}}^{(k)}$ using the BLUP formula.

until Convergence of $\widehat{\boldsymbol{b}}^{(k)}.$

Algorithm 4 Mixed Random Forest (algorithm for fitting a single tree)

Input: $\boldsymbol{y},\;X,\;Z.$

Output: $f(\cdot),\;\widehat{\sigma}^{2}_{b},\widehat{\sigma}^{2}.$

Initialization: $g_{1}=\\{0,1\\}^{J}.$ Also, estimate
$\delta=\sigma^{2}/\sigma^{2}_{b}.$

repeat

  1. 1.

Given the predefined stopping rules, for each node $s\in\mathcal{S},$ find the
following parameters:

$\displaystyle\small\tilde{j_{s}},\tilde{\sigma}^{2}=$
$\displaystyle\underset{j_{s}\in
J_{s},\sigma^{2}\in\mathbb{R}^{+}}{\text{argmax}}\;\ell(\boldsymbol{y};\sigma^{2},\mu_{s}^{l},\mu_{s}^{r},\\{\mu_{l}\\}_{l\in\mathcal{S}/s},g_{s}\cap(x_{j_{s}}=0),g_{s}\cap(x_{j_{s}}=1),\\{g_{l}\\}_{l\in\mathcal{S}/s}),$

where $\mu_{s}^{l},\;\mu_{s}^{r}$ and $\\{\mu_{l}\\}_{l\in\mathcal{S}/s}$ are
the GLS estimators for $\\{y_{i}|\boldsymbol{x}_{i}\in
g_{s}\cap(x_{i,j_{s}}=0)\\},$ $\\{y_{i}|\boldsymbol{x}_{i}\in
g_{s}\cap(x_{i,j_{s}}=1)\\},$ and $\\{y_{i}|\boldsymbol{x}_{i}\in
g_{l}\\}_{l\in\mathcal{S}/s}$ subsets respectively.

  2. 2.

Update $f(\cdot):$ replace each node $s\in\mathcal{S}$ by
$\\{g_{s}\cap(x_{\tilde{j}_{s}}=0),\mu_{s}^{l}\\}$ and
$\\{g_{s}\cap(x_{\tilde{j}_{s}}=1),\mu_{s}^{r}\\}.$

until Stopping rules are satisfied $\forall s\in\mathcal{S}.$

##  D Numerical Results

This appendix presents additional results that are related to Section 5.1, as
well as detailed information relating the settings in Section 5.2.

###  D.1 Regression Tree - New Random Effects
($\boldsymbol{b}^{*}\perp\boldsymbol{b}$)

####  D.1.1 In-Sample Error Setting:

Given the simulation setting that is described in Section 5.1.1, Figure 4,
left panel, presents the effect of the cluster size ($n_{c}$), on the
performance of RETCO. As can be seen in the figure, for larger block size
($n_{c}=150$) the average error difference[%] is smaller.

Figure 4: Boxplots of the error difference[%]. Left: In-sample error setting,
$\boldsymbol{b}^{*}\perp\boldsymbol{b},$ $n=300,\;\sigma_{b}^{2}=1$ and
different $n_{c}.$ The means are $-9.9,\;-8.0$ and $-4.3.$ Middle: In-sample
error setting, $\boldsymbol{b}^{*}\perp\boldsymbol{b},$
$\sigma_{b}^{2}=1\;n_{c}=50$ and different $n.$ The means are $-2.9,\;-9.9$
and $-5.0.$ Right: Generalization error setting,
$\boldsymbol{b}^{*}\perp\boldsymbol{b},$ $n=300,\;n_{c}=50$ and different
$\sigma_{b}^{2}.$ The means are $-2.8,\;-11.0$ and $-7.3.$

Figure 4, middle panel, presents the effect of the sample size on the
performance of RETCO. As can be seen, RETCO performs better for all the
settings. Also, as expected, the error difference[%] variance is smaller for
larger sample sizes. The variance depends on the maximal depth of the tree,
which was set to three. Tree with three levels has potentially eight
predictors, which is a large amount of predictors when $n=100,$ but small when
$n=1000.$ Therefore, when $n=100$ the trees are noisy for both algorithms, and
their relative difference is noisy as well.

####  D.1.2 Generalization Prediction Error Setting:

In order to analyze RETCO performance in generalization prediction error
setting, the test set setting was changed such that the prediction set
covariates are nonidentical to the training sample covariates (but are sampled
from the same distribution). As was described in the paper, $CV_{c}$ loss
function estimates the generalization error unbiasedly by correcting the
standard CV error. Therefore, $CV_{c}$ loss function is used in RETCO and CV
loss function is used for the standard regression tree algorithm. Figure 4,
right panel, presents the error difference[%] for different $\sigma^{2}_{b}.$

###  D.2 Comparison With RE-EM

Figure 5, left figure, presents the error difference between RETCO and RE-EM
for the scenario when $\boldsymbol{b}^{*}\perp\boldsymbol{b}.$ All the other
setting details are the same as in Section 5.1.1. Similarly, Figure 5, right
figure, presents the error difference between RETCO and RE-EM for the scenario
when $\boldsymbol{b}^{*}=\boldsymbol{b}.$ All the other setting details are
the same as in Section 5.1.2. As can be seen in Figure 5, RETCO performs
better than RE-EM.

Figure 5: Boxplot of Error difference[%] between RETCO and RE-EM. Left figure:
$\boldsymbol{b}^{*}\perp\boldsymbol{b}^{*}$ setting. Right Figure:
$\boldsymbol{b}^{*}=\boldsymbol{b}^{*}$ setting.

###  D.3 Real Data Analysis

Detailed information about the settings of the real data analyses is given
below.

General information:

  * •

The number of covariates that were sampled at each split in the RF
implementations is $round\big{(}log_{2}(p)\big{)}.$

  * •

The number of trees that are used for the RF is not fixed. The fitting process
was stooped once the overall RF error was converged in both algorithms, RETCO
and the standard regression tree.

Specific information for each data set analysis:

  * •

FIFA

    * –

Dependent variable: Player’s market values

    * –

Cluster variable: Player’s club

    * –

Loss function type: CV

    * –

Number of covariates: $11$

    * –

Regression tree depth: $5$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $60$

    * –

Comments: The covariates that are used in this analysis are: ’Age’, ’Overall’,
’Potential’, ’Wage’, ’Special’, ’Preferred Foot’, ’International Reputation’,
’Weak Foot’, ’Skill Moves’, ’Height’, ’Weight’. Other variables have many
missing values or are irrelevant.

  * •

Crimes (Communities and Crime in US)

    * –

Dependent variable: Violent crimes in US communities per population size

    * –

Cluster variable: State

    * –

Loss function type: CV

    * –

Number of covariates: $100$ (all the available covariates were used)

    * –

Regression tree depth: $5$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $55$

    * –

Comments: -

  * •

Korea Temperature (’bias correction of numerical prediction model temperature
forecast’)

    * –

Dependent variable: Daily maximum temperature at several sites in South Korea

    * –

Cluster variable: Day

    * –

Loss function type: Cp (since training and test set measurements are sampled
from the same sites)

    * –

Number of covariates: $4$

    * –

Regression tree depth: $5$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $70$

    * –

Comments:

      * *

The covariates that are used in this analysis are: ’Present_Tmin’,’DEM’,
’Slope’, ’Solar radiation’. Other variables in this data set are models’s
scores of the data set supplier, which are based on previous dependent
variable measurements (and therefore cannot be used in LMM framework)

      * *

Exponential kernel covariance function was used. Maximal temperature of
different days are assumed to be uncorrelated.

      * *

The original data set contains records from July and August. Due to many
missing values in July along the years, only records from August are analyzed.
Also, sites with missing values along the years were omitted.

  * •

California Housing

    * –

Dependent variable: Values of houses in California

    * –

Cluster variable: Block’s cluster (blocks with the same coordinate values)

    * –

Loss function type: CV

    * –

Number of covariates: $6$ (all the available covariates were used)

    * –

Regression tree depth: $5$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $50$

    * –

Comments: -

  * •

Parkinson’s Disease Telemonitoring

    * –

Dependent variable: Total UPDRS score, which is a score of Parkinson’s Disease
progression

    * –

Cluster variable: Patient

    * –

Loss function type: Cp (all the individuals receive approximately the same
covariate values)

    * –

Number of covariates: $18$

    * –

Regression tree depth: $5$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $80$

    * –

Comments: All the supplied covariates were used except the motor_UPDRS (which
its relation with the dependent variable is not fully clear to us)

  * •

Wages

    * –

Dependent variable: Average hourly wages

    * –

Cluster variable: Employee

    * –

Loss function type: CV

    * –

Number of covariates: $6$ (all the available covariates were used)

    * –

Regression tree depth: $4$

    * –

Minimum number of observations in a node: $3$

    * –

Number of trees for RF: $80$

    * –

Comments: -

## References

  * Akaike (1974) Hirotugu Akaike.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6):716–723, 1974.
  * Breiman (2001) Leo Breiman.  Random forests.  _Machine learning_ , 45(1):5–32, 2001.
  * Breiman et al. (1984) Leo Breiman, Jerome Friedman, Charles J Stone, and Richard A Olshen.  _Classification and regression trees_.  CRC press, 1984.
  * Caywood et al. (2017) Matthew S Caywood, Daniel M Roberts, Jeffrey B Colombe, Hal S Greenwald, and Monica Z Weiland.  Gaussian process regression for predictive but interpretable machine learning models: An example of predicting mental workload across tasks.  _Frontiers in human neuroscience_ , 10:647, 2017.
  * Cho et al. (2020) Dongjin Cho, Cheolhee Yoo, Jungho Im, and Dong-Hyun Cha.  Comparative assessment of various machine learning-based bias correction methods for numerical weather prediction model forecasts of extreme air temperatures in urban areas.  _Earth and Space Science_ , 7(4), 2020.
  * Coull et al. (2001) Brent A Coull, Joel Schwartz, and MP Wand.  Respiratory health and air pollution: additive mixed model analyses.  _Biostatistics_ , 2(3):337–349, 2001.
  * Criminisi et al. (2011) Antonio Criminisi, Jamie Shotton, and Ender Konukoglu.  Decision forests for classification, regression, density estimation, manifold learning and semi-supervised learning.  _Microsoft Research Cambridge, Tech. Rep. MSRTR-2011-114_ , 5(6):12, 2011.
  * Fokkema et al. (2018) Marjolein Fokkema, Niels Smits, Achim Zeileis, Torsten Hothorn, and Henk Kelderman.  Detecting treatment-subgroup interactions in clustered data with generalized linear mixed-effects model trees.  _Behavior research methods_ , 50(5):2016–2034, 2018.
  * Freund et al. (1999) Yoav Freund, Robert Schapire, and Naoki Abe.  A short introduction to boosting.  _Journal-Japanese Society For Artificial Intelligence_ , 14(771-780):1612, 1999.
  * Friedman (2001) Jerome H Friedman.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , pages 1189–1232, 2001.
  * Goovaerts (1999) Pierre Goovaerts.  Geostatistics in soil science: state-of-the-art and perspectives.  _Geoderma_ , 89(1-2):1–45, 1999.
  * Hajjem et al. (2014) Ahlem Hajjem, François Bellavance, and Denis Larocque.  Mixed-effects random forest for clustered data.  _Journal of Statistical Computation and Simulation_ , 84(6):1313–1328, 2014.
  * Harville et al. (1976) David Harville et al.  Extension of the Gauss-Markov theorem to include the estimation of random effects.  _The Annals of Statistics_ , 4(2):384–395, 1976\.
  * Hastie et al. (2009) Trevor Hastie, Robert Tibshirani, and JH Friedman.  _The elements of statistical learning: data mining, inference, and prediction_.  New York, NY: Springer, 2009.
  * Hengl et al. (2018) Tomislav Hengl, Madlene Nussbaum, Marvin N Wright, Gerard BM Heuvelink, and Benedikt Gräler.  Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables.  _PeerJ_ , 6, 2018.
  * Hodges and Sargent (2001) James S Hodges and Daniel J Sargent.  Counting degrees of freedom in hierarchical and other richly-parameterised models.  _Biometrika_ , 88(2):367–379, 2001.
  * Mallows (1973) Colin L Mallows.  Some comments on c p.  _Technometrics_ , 15(4):661–675, 1973.
  * Ngufor et al. (2019) Che Ngufor, Holly Van Houten, Brian S Caffo, Nilay D Shah, and Rozalina G McCoy.  Mixed effect machine learning: A framework for predicting longitudinal change in hemoglobin a1c.  _Journal of biomedical informatics_ , 89:56–67, 2019.
  * Pace and Barry (1997) R Kelley Pace and Ronald Barry.  Sparse spatial autoregressions.  _Statistics & Probability Letters_, 33(3):291–297, 1997.
  * Painsky and Rosset (2016) Amichai Painsky and Saharon Rosset.  Cross-validated variable selection in tree-based methods improves predictive performance.  _IEEE transactions on pattern analysis and machine intelligence_ , 39(11):2142–2153, 2016.
  * Prokhorenkova et al. (2017) Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, and Andrey Gulin.  Catboost: unbiased boosting with categorical features.  _arXiv preprint arXiv:1706.09516_ , 2017.
  * Rabinowicz and Rosset (2020) Assaf Rabinowicz and Saharon Rosset.  Cross-validation for correlated data.  _Journal of the American Statistical Association_ , pages 1–14, 2020\.
  * Rasmussen (2003) Carl Edward Rasmussen.  Gaussian processes in machine learning.  In _Summer school on machine learning_ , pages 63–71. Springer, 2003\.
  * Redmond and Baveja (2002) Michael Redmond and Alok Baveja.  A data-driven software tool for enabling cooperative information sharing among police departments.  _European Journal of Operational Research_ , 141(3):660–678, 2002.
  * Sela and Simonoff (2012) Rebecca J Sela and Jeffrey S Simonoff.  Re-em trees: a data mining approach for longitudinal and clustered data.  _Machine learning_ , 86(2):169–207, 2012.
  * Singer et al. (2003) Judith D Singer, John B Willett, John B Willett, et al.  _Applied longitudinal data analysis: Modeling change and event occurrence_.  Oxford university press, 2003.
  * Speiser et al. (2019) Jaime Lynn Speiser, Bethany J Wolf, Dongjun Chung, Constantine J Karvellas, David G Koch, and Valerie L Durkalski.  Bimm forest: A random forest method for modeling clustered and longitudinal binary outcomes.  _Chemometrics and Intelligent Laboratory Systems_ , 185:122–134, 2019.
  * Stephan et al. (2015) Johannes Stephan, Oliver Stegle, and Andreas Beyer.  A random forest approach to capture genetic effects in the presence of population structure.  _Nature communications_ , 6:7432, 2015.
  * Stone (1974) Mervyn Stone.  Cross-validatory choice and assessment of statistical predictions.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 36(2):111–133, 1974.
  * Tsanas et al. (2009) Athanasios Tsanas, Max Little, Patrick McSharry, and Lorraine Ramig.  Accurate telemonitoring of parkinson’s disease progression by non-invasive speech tests.  _Nature Precedings_ , pages 1–1, 2009.
  * Vaida and Blanchard (2005) Florin Vaida and Suzette Blanchard.  Conditional akaike information for mixed-effects models.  _Biometrika_ , 92(2):351–370, 2005.
  * Verbeke (1997) Geert Verbeke.  Linear mixed models for longitudinal data.  In _Linear mixed models in practice_ , pages 63–153. Springer, 1997\.
  * Westveld et al. (2011) Anton H Westveld, Peter D Hoff, et al.  A mixed effects model for longitudinal relational and network data, with applications to international trade and conflict.  _The Annals of Applied Statistics_ , 5(2A):843–872, 2011.
  * Wolfinger and O’connell (1993) Russ Wolfinger and Michael O’connell.  Generalized linear mixed models a pseudo-likelihood approach.  _Journal of statistical Computation and Simulation_ , 48(3-4):233–243, 1993.