SMAppendix References

# Analyze Additive and Interaction Effects via Collaborative Trees

Chien-Ming Chi
Academia Sinica  Chien-Ming Chi is Assistant Research Fellow, Institute of
Statistical Science, Academia Sinica, Taipei 11529, Taiwan (Email:
xbbchi@stat.sinica.edu.tw). This work was supported by grant
111-2118-M-001-012-MY2 from the National Science and Technology Council,
Taiwan.

###### Abstract

We present Collaborative Trees, a novel tree model designed for regression
prediction, along with its bagging version, which aims to analyze complex
statistical associations between features and uncover potential patterns
inherent in the data. We decompose the mean decrease in impurity from the
proposed tree model to analyze the additive and interaction effects of
features on the response variable. Additionally, we introduce network diagrams
to visually depict how each feature contributes additively to the response and
how pairs of features contribute interaction effects. Through a detailed
demonstration using an embryo growth dataset, we illustrate how the new
statistical tools aid data analysis, both visually and numerically. Moreover,
we delve into critical aspects of tree modeling, such as prediction
performance, inference stability, and bias in feature importance measures,
leveraging real datasets and simulation experiments for comprehensive
discussions. On the theory side, we show that Collaborative Trees, built upon
a “sum of trees” approach with our own innovative tree model regularization,
exhibit characteristics akin to matching pursuit, under the assumption of
high-dimensional independent binary input features (or one-hot feature
groups). This newfound link sheds light on the superior capability of our tree
model in estimating additive effects of features, a crucial factor for
accurate interaction effect estimation.

Running title: CTE

Key words: Sensitivity analysis; Feature interaction; Visual network diagram;
Matching pursuit; Sum of trees.

##  1 Introduction

Insightful observations derived from statistical analyses of variable
associations often lead to the formulation of good data-driven questions. For
instance, when investigating how variables such as gender, college majors, and
job categories influence the wage function of individual workers, analyzing
the additive and interaction effects of these variables prompts the
development of hypotheses for further exploration. Suppose we conduct an
analysis and find that the influence of gender on wage rates is pronounced and
primarily additive, suggesting $\textnormal{Wage}=f_{1}(\textnormal{Other
Variables})+f_{2}(\textnormal{Gender})$, where $f_{1}$ and $f_{2}$ are some
functions. When other variables fully characterize workers’ capacity, this
analysis raises the possibility of unequal pay for equal work, highlighting a
social issue deserving attention. The example demonstrates how important
subsequent analysis and inferences can stem from a thorough examination of
variable effects, thereby motivating our efforts to establish a comprehensive
analysis of the impact of variables on responses.

This paper introduces a new tree model called Collaborative Trees and its
bagging version to estimate effects of variables. Our estimation framework
shares similarities with the Mean Decrease in Impurity (MDI) feature
importance measure of Random Forests [7]. Notably, we make a significant
advancement by offering a more detailed analysis of both additive and
interaction effects of input features, providing insights beyond the overall
importance of each feature. Among various works [30, 35, 3, 2] concerning
interactions, to name just a few, our work is closely aligned with Sobol’
indices [35, 33] for sensitivity analysis, which decompose the variance of
response into the contribution shares of each feature and the interaction
between all features. This decomposition is often referred to as the ANOVA
decomposition [17, 2]. Unlike Sobol’ indices, our analysis emphasizes the
consistency of effect estimation and is attentive to applications involving
high-dimensional and correlated input features. Additionally, to enhance the
presentation of results from our analysis, which may yield a substantial
number of estimated effects, we have introduced network diagrams for visually
summarizing and screening the essential signals of interest. These diagrams
aim to provide a comprehensive assessment of the statistical associations
among all features, distinguishing them from plots of partial dependence [17]
and accumulated local effects [3], which typically focus on analysis of
specific features.

An example of a network diagram is provided in Figure 1 for the analysis of
temperature-dependent sex determination in animals, which is crucial for
conserving endangered species impacted by climate change [19, 1]. In Figure 1,
the ratio of females to total births (z-axis) for newborns from five distinct
species (y-axis), including various crocodile and turtle types, are shown. The
x-axis indicates the incubation temperature. From the right panel of Figure 1,
most species exhibit higher ratios of females to total births at higher
incubation temperatures. Particularly, for species Lepidochelys olivacea,
Chelonia mydas, Caretta caretta, the sex ratios are well approximated by a
logistic function with a single variable temperature. In fact, a majority of
species in the data has a similar temperature-dependent sex determination,
suggesting that additive effect may account for a significant proportion of
the overall importance measure of temperature in determining sex ratios. On
the other hand, we observe from the samples of Chelydra serpentina and Emys
orbicularis that sex ratios varies across different species, which is clearly
an interaction effect between temperature and species on the sex ratio.

Our network diagram summarizes the above observation. In the diagram of Figure
1, large circles represent features with high importance measures, thick edges
between circles indicate strong interactions, and blue (red) circles highlight
features with relatively strong additive (interaction) effects. Further
diagram details are provided in Section 3.2. Our diagram suggests the observed
pattern—the strong interaction between species and temperature for the 5
selected species—is consistent across all species in the dataset.
Additionally, it illustrates the strong additive effect of temperature on the
sex ratio. Moreover, the diagram prompts further analysis on the importance of
area and incubation periods (days), as discussed in Section 6.1. Notably, this
new diagram can be used independently from Collaborative Trees, provided the
required feature importance measures are available. However, a key advantage
of using Collaborative Trees with the diagram is their ability to effectively
distinguish truly important effects from unimportant ones. This capability is
crucial for visually screening and highlighting relevant information in the
diagram. This merit of Collaborative Trees is further discussed in Section 3.1
and Section 5.1.

The novel visualization tool arises from our in-depth analysis of feature
association due to Collaborative Trees. Trees, known for their ability to
reveal interactions between features, are utilized for analyzing both additive
and interaction effects. However, trees are acknowledged to have limitations
in capturing additive signals effectively [22, 39]. To augment Collaborative
Trees’ capability in handling additive components, several tree-based
techniques have been proposed, including incorporating linear regression into
tree nodes [15] and training a sum of additive trees. In this work, we
leverage the ”sum of trees” approach. Details of this enhancement are provided
in the following section.

(a)

(b)

Figure 1: Embryo growth dataset (right panel) and the network diagram (left
panel) for analyzing temperature-dependent sex determination. The z-axis,
y-axis, and x-axis are respectively the sex ratio, species name, and
incubation temperature.

###  1.1 Learning numerous additive components

A single decision tree is known to have limited capacity for fitting linear
combinations of features [22], which cannot be improved by merely applying
methods like bagging or column subsampling [7]. To overcome this problem, the
proposed Collaborative Trees, adopts a specific tree-growing approach where
new trees in the model are grown from the root nodes, utilizing residuals from
the current round of updates. Intuitively, this tree-growing approach allows
trees to extract the remaining additive signals in the residual by splitting
on shallow nodes of new trees, which is often more efficient than learning
additive signals by splitting deep nodes of the existing trees. This idea of
learning from residuals is widely used for building tree models, and the
resulting trees are often referred to a sum of trees [18, 40, 12]. See, for
example, BART [12], FIGS [40], and Gradient Boosting Decision Trees (GBDT)
such as XGBoost [9], lightGBM [23], and MART [18]. Among them, Collaborative
Trees insist on growing a new root node after splitting each root node, GBDT
insists on finishing the current tree before growing a new root node, while
BART and FIGS rely on a sample-driven decision for each update, allowing them
to potentially split on a new root node at any round; see Figure 2 for a
graphical illustration.

Figure 2: Collaborative Trees grow a new root node after splitting each root
node, GBDT finishes the current tree before growing a new root node, while
BART and FIGS rely on a sample-driven decision for each update, allowing them
to potentially split on a new root node at any round. The number of potential
root nodes is a hyperparameter.

Our update priority and other model regularization techniques (see Section 2.2
for the detail) enhance Collaborative Trees as a specific variant of matching
pursuit [32, 13]. These enhancements ensure the model’s capability in
capturing numerous additive components. Notably, our matching pursuit approach
involves tree heredity and projections onto multiple feature groups to assess
feature interactions, thereby differing from the standard matching pursuit
methods [32] in the literature. We formally show that the sequential decisions
of split coordinates of Collaborative Trees can be viewed as a feature
selection path of our matching pursuit for applications with independent
binary features or feature groups with one-hot indicators (e.g., resulting
from binning continuous features). Additionally, our analysis shows that: 1)
Extracting all additive signals (i.e., sure screening [14]) is crucial for
subsequently learning interaction effects. 2) To sure screen important
components, the numbers of significant additive and interaction components are
permitted to increase at some polynomial order with respect to the sample
size. 3) The dimension of the input high-dimensional feature vector is
permitted to increase at a rate on the order of the sample size raised to an
arbitrarily large constant. To the best of our knowledge, we are among the
first to formally establish this result, providing an interesting insight into
the tree model building method “learning from residuals,” as depicted in
Figure 2. Moreover, our theory allows us to develop new feature importance
measures for both additive effects and interaction effects of features. These
measures serve as crucial components in constructing our network diagrams,
offering a comprehensive visualization of the relationships and contributions
of individual features in the model.

The rest of the paper is organized as follows. The Collaborative Trees
algorithm is introduced in Section 2. The additive and interaction effects and
their estimation are presented in Section 3. Theoretical results for
ollaborative Trees and matching pursuit are in Section 4. Simulation
experiments and real data studies are in Sections 5–6.

###  1.2 Notation

Let $(\Omega,\mathcal{F},\mathbb{P})$ be the probability triple, and
$\boldsymbol{X}=(X_{1},\cdots,X_{p})^{{}^{T}}$ represent the $p$-dimensional
random feature vector in the probability space. For
$H\subset\\{1,\cdots,p\\}$, $\boldsymbol{X}_{H}=(X_{j},j\in H)^{{}^{T}}$ and
$\boldsymbol{X}_{-H}=(X_{j},j\not\in H)^{{}^{T}}$. Similar notation is used
for their sample counterparts. Define $a\vee b=\max\\{a,b\\}$ for
$(a,b)\in\mathbb{R}^{2}$. Let $\boldsymbol{1}_{A}=\boldsymbol{1}\\{A\\}$ be
the indicator function, where $\boldsymbol{1}_{A}=1$ if $A$ is true, and
$\boldsymbol{1}_{A}=0$ otherwise. Summation over an empty set is defined as
zero. Sets contain unique elements, e.g., $\\{j,j,k\\}=\\{j,k\\}$. Integer
parameters, like tree depth, are rounded to the nearest integer when
necessary. For two real sequences $a_{n}$’s and $b_{n}$’s, define
$a_{n}=O(b_{n})$ to mean that
$\limsup_{n\rightarrow\infty}|a_{n}/b_{n}|<\infty$; $a_{n}=o(b_{n})$ implies
that $\lim_{n\rightarrow\infty}|a_{n}/b_{n}|=0$. Lastly, $\log a=\log_{2}a$
for real $a>0$.

##  2 Collaborative Trees

The proposed tree model is termed as Collaborative Trees, with its bootstrap
aggregation ensemble version called Collaborative Trees Ensemble (CTE). The
name reflects the development of our tree model, which is briefly introduced
here. 1.) Our tree model is a sum of predictive trees, which can be seen as
collaboration between trees in making improved prediction. 2.) Child nodes
with the same parent node collaborate, in the sense that these child nodes are
updated together, as detailed in Section 2.2. 3.) Features within the same
group collaborate, meaning that Collaborative Trees split on feature groups
rather than individual features when updating nodes. The concept of feature
groups is introduced in Section 2.1. Our tree model regularization is
introduced in detail in (2).

###  2.1 Feature groups of one-hot indicators

Collaborative Trees apply to general features, but our primary inference
applications consider binary features and indicator feature groups, in which
the feature groups may be formed by transforming categorical features and
binning continuous features into one-hot vectors. The use of feature groups
significantly broadens the applications of our inference method. For example,
in the embryo growth data set, the categorical variable species is represented
by an one-hot vector, and the continuous variables incubation time and
temperature are respectively binned into two one-hot vectors. See Section 6.1
for details.

Specifically, let set $\mathcal{X}(m)\subset\\{1,\cdots,p\\}$ with
$\texttt{\\#}\mathcal{X}(m)>1$ denote a group of features of interest for each
$m\in\\{1,\cdots,M_{0}\\}$ for some $0\leq M_{0}\leq M$, where $M$ is the
number of groups. Here, $\mathcal{X}(l)\cap\mathcal{X}(k)=\emptyset$ for every
$1\leq l<k\leq M_{0}$ and
$\sum_{j\in\mathcal{X}(m)}\boldsymbol{1}_{X_{j}=1}=1$ for $1\leq m\leq M_{0}$.
For notation completeness, each $\mathcal{X}(m)$ has a single binary feature
(e.g., the gender variable) for $M_{0}<m\leq M$ such that
$\cup_{m=1}^{M}\mathcal{X}(m)=\\{1,\cdots,p\\}$. To clarify, two special cases
are $(M_{0},M)=(0,p)$ and $M_{0}=M$, corresponding respectively to the
scenario with only single features and the scenario without any single
features. Throughout this paper, “feature groups” is used with awareness that
they may include single features.

###  2.2 Collaborative Trees

Collaborative Trees employ $K$ collaborative decision trees, denoted as
$T_{k}:\mathbb{R}^{p}\mapsto\mathbb{R}$ for $k\in\\{1,\cdots,K\\}$. These
trees are trained on the sample
$\\{\boldsymbol{X}_{i}=(X_{i1},\cdots,X_{ip})^{{}^{T}},Y_{i}\\}_{i=1}^{n}$,
where each $(\boldsymbol{X}_{i},Y_{i})$ and $(\boldsymbol{X},Y)$ are
independently and identically distributed (i.i.d.). For any test point
$\vv{x}\in\mathbb{R}^{p}$, Collaborative Trees generate predictions by summing
up individual tree predictions: $\sum_{k=1}^{K}T_{k}(\vv{x})$, which is often
referred to a sum of trees [18, 40, 12]. In what follows, we will introduce
the recursive updating algorithm for training the $K$ collaborative decision
trees $T_{k}(\cdot)$.

Let us consider a hyper rectangle $C_{0}\subset\mathbb{R}^{p}$ , also known as
a node. When splitting $C_{0}$ on the $m$th feature group, the child nodes are
denoted by $C_{1}=\\{\vv{x}\in C_{0}:x_{j}>\widehat{c}\\}$ and
$C_{2}=\\{\vv{x}\in C_{0}:x_{j}\leq\widehat{c}\\}$ if
$\mathcal{X}(m)=\\{j\\}$, where the split point $\widehat{c}$ is decided by
CART [7] as in (4) below. If $\texttt{\\#}\mathcal{X}(m)>1$, the child nodes
are given by $C_{j}=\\{\vv{x}\in C_{0}:x_{j}>0\\}$ for $j\in\mathcal{X}(m)$.
In the latter case, the indices of child nodes $C_{j}$’s may not be
consecutive integers. Child nodes with the same parent node are referred to as
associated nodes. For example, in Figure 3, $\\{a_{1},a_{2}\\}$,
$\\{b_{1},b_{2},b_{3}\\}$, and $\\{c_{1},c_{2},c_{3}\\}$ are three sets of
associated nodes. For completeness, each root node (the feature space) is also
seen as a set of associated nodes. Collaborative Trees are trained using
recursive splitting. During each iteration of updating Collaborative Trees,
the algorithm decides a set of associated nodes to be updated, where each node
in this set of associated nodes is split on some feature group. In addition,
the feature group to be split should be the same across these associated
nodes. For a graphical example, the set of associated nodes $a_{1}$ and
$a_{2}$ are split on the same $m_{1}$th feature group at STEP 1 in Figure 3.

Figure 3: Here, $K=1$ and $W=\\{A\\}$ at STEP 1, where
$A=\\{(a_{1},1),(a_{2},1)\\}$; elements in $A$ are pairs of (node, index of
tree having the node). An update on $A$ is made in three steps: At STEP 1,
$\widehat{Q}=A$, and $\widehat{m}=m_{1}$ is decided by (2), where
$\texttt{\\#}\mathcal{X}(m_{1})=3$. At STEP 2, the tree is updated as in Step
12 of Algorithm 1, and $A$ is removed from $W$, leading to $W=\emptyset$. At
STEP 3, two new sets of associated nodes,
$B=\\{(b_{1},1),(b_{2},1),(b_{3},1)\\}$ and
$C=\\{(c_{1},1),(c_{2},1),(c_{3},1)\\}$, are added to $W$, leading to
$W=\\{B,C\\}$.

In this work, $W$ denotes the set of all eligible node sets for splitting at
each round, along with the indices of their corresponding trees, which are
needed for Algorithm 1. For example, at the first round,

$W=\Big{\\{}\big{\\{}(\textnormal{Root
Node},1)\big{\\}},\cdots,\big{\\{}(\textnormal{Root
Node},K)\big{\\}}\Big{\\}},$ (1)

where the $k$th element in $W$ represents a set containing only the root node
of the $k$th tree, with the tree index indicated in parentheses. As the
algorithm progresses, an element in $W$ may encompass multiple associated
nodes; associated nodes in each element are updated at the same round. A
graphical example illustrating how each element in $W$ is removed and new
elements are included at each round is shown in Figure 3. During each round,
given $W$ and the residuals
$\\{\widehat{r}_{i}=Y_{i}-\sum_{k=1}^{K}T_{k}(\boldsymbol{X}_{i})\\}_{i=1}^{n}$,
the optimization in (2) decides the best update. At this round, we obtain a
set of associated nodes $\widehat{Q}\in W$ to be updated and the
$\widehat{m}$th feature group to be split as follows, with ties broken
randomly.

$\begin{split}(\widehat{Q},\widehat{m})&=\operatorname*{arg\,max}_{Q\in W,\
m\in\\{1,\cdots,M\\}}\Big{[}\Big{(}\sum_{(C,k)\in Q}\textnormal{
SplitScore}(C,m)\Big{)}-\lambda(W,Q)\Big{]},\end{split}$ (2)

where $\lambda(W,Q)$ is defined as follows: $\lambda(W,Q)=\infty$ if some
$Q^{\dagger}\in W$ includes a root node (at depth zero) but $Q$ does not, or
some $Q^{\dagger}\in W$ includes associated nodes at depth $1$ while $Q$ has
deeper associated nodes. Otherwise, $\lambda(W,Q)=0$. The term $\lambda(W,Q)$
ensures the updating priority: 1. Root nodes. 2. Associated nodes at depth
one. 3. Other associated nodes with depth deeper than one. The rationale
behind our tree regularization aligns with the idea of tree boosting [17],
which also emphasizes learning from shallow nodes.

In addition, the split score, denoted as $\textnormal{SplitScore}(C_{0},m)$,
is the result of splitting a node $C_{0}$ using the $m$th feature group, and
it is given by

$\begin{cases}\max_{a\in N_{0}}\left[\left(\sum_{i\in
N_{0}}\widehat{r}_{i}^{2}\right)-\sum_{l=1}^{2}\sum_{i\in
N_{l}(a)}(\widehat{r}_{i}-\frac{1}{1\vee\texttt{\\#}N_{l}(a)}\sum_{i\in
N_{l}(a)}\widehat{r}_{i})^{2}\right]\ \ \textnormal{ if
}\texttt{\\#}\mathcal{X}(m)=1,\\\ \left(\sum_{i\in
N_{0}}\widehat{r}_{i}^{2}\right)-\sum_{j\in\mathcal{X}(m)}\sum_{i\in
N_{j}}(\widehat{r}_{i}-\frac{1}{1\vee\texttt{\\#}N_{j}}\sum_{i\in
N_{j}}\widehat{r}_{i})^{2}\hskip 79.6678pt\textnormal{ if
}\texttt{\\#}\mathcal{X}(m)>1,\end{cases}$ (3)

where $N_{j}=\\{i:\boldsymbol{X}_{i}\in C_{0},X_{ij}>0\\}$ if
$\texttt{\\#}\mathcal{X}(m)>1$, $N_{1}(a)=\\{i:\boldsymbol{X}_{i}\in
C_{0},X_{iq}>X_{aq}\\}$ and $N_{2}(a)=\\{i:\boldsymbol{X}_{i}\in
C_{0},X_{iq}\leq X_{aq}\\}$ if $\mathcal{X}(m)=\\{q\\}$ for some
$q\in\\{1,\cdots,p\\}$, and $N_{0}=\\{i:\boldsymbol{X}_{i}\in C_{0}\\}$. See
Algorithm 1 below for further details, where the validity of nodes for
splitting and tree update is introduced in Section 2.3.

For clarity in Step 9 of Algorithm 1, when splitting $C_{0}$ with respect to
$\mathcal{X}(\widehat{m})=\\{q\\}$ for some $q\in\\{1,\cdots,p\\}$, we denote
$\widehat{a}$ as the sample index in $N_{0}$ that maximizes (3) (note that
$\texttt{\\#}\mathcal{X}(\widehat{m})=1$). In this case, the two corresponding
child nodes are given by

$C_{1}=\\{\vv{x}\in C_{0}:x_{q}>\widehat{c}\\}\textnormal{ and
}C_{2}=\\{\vv{x}\in C_{0}:x_{q}\leq\widehat{c}\\},\textnormal{ where
}\widehat{c}=X_{\widehat{a}q}.$ (4)

Input: Initial decision trees with
$\sup_{\vv{x}\in\mathbb{R}^{p}}|T_{k}(\vv{x})|=0$, and $W$ as in (1).

Output: $K$ collaborative decision trees $T_{1}(\cdot),\cdots,T_{K}(\cdot)$.

1

2while  _$W\not=\emptyset$_ do

3 Get $\widehat{Q}$ and $\widehat{m}$ from (2) with
$\widehat{r}_{i}=Y_{i}-\sum_{k=1}^{K}T_{k}(\boldsymbol{X}_{i})$ and $W$.

4 Remove $\widehat{Q}$ from $W$.

5 for  _$(C_{0},\widehat{k})\in\widehat{Q}$_ do

6

7 if  _$\texttt{\\#}\mathcal{X}(\widehat{m}) >1$_ then

8 Denote the $\texttt{\\#}\mathcal{X}(\widehat{m})$ child nodes by
$C_{j}=\\{\vv{x}\in C_{0}:x_{j}>0\\}$ for $j\in\mathcal{X}(\widehat{m})$

9 Append $\\{(C_{j},\widehat{k}):j\in\mathcal{X}(\widehat{m})\textnormal{ and
}C_{j}\textnormal{ is a valid node for splitting}\\}$ to $W$.

10 else

11

12 Let $C_{1}$ and $C_{2}$ be child nodes obtained as in (4).

13 Append $\\{(C_{l},\widehat{k}):l\in\\{1,2\\}\textnormal{ and
}C_{l}\textnormal{ is a valid node for splitting}\\}$ to $W$.

14 end if

15 $T_{\widehat{k}}(\vv{x})\leftarrow
T_{\widehat{k}}(\vv{x})+\sum_{j\in\\{l:\textnormal{ valid child nodes
}C_{l}\textnormal{ for tree update}\\}}\frac{\boldsymbol{1}\\{\vv{x}\in
C_{j}\\}}{1\vee\texttt{\\#}\\{i:\boldsymbol{X}_{i}\in
C_{j}\\}}\sum_{i\in\\{\boldsymbol{X}_{i}\in C_{j}\\}}\widehat{r}_{i}$.

16 end for

17

18 end while

Algorithm 1 Recursive Splitting for Growing $K$ Collaborative Trees

###  2.3 Collaborative Trees Ensemble

In this section, we introduce additional algorithm implementations in the
following four points for Collaborative Trees Ensemble. 1.) To enhance
prediction, we employ Collaborative Trees with bootstrap aggregating, or
bagging. 2.) The validity of a node for splitting and/or tree updating depends
on the number of subsamples in the node and its current depth. Our tree model
employs hyperparameters min_samples_split, min_samples_leaf, and max_depth for
node validation. 3.) To optimize computation time, after the initial $2K$
updates, we only update a subset of nodes in the update node list $W$. 4.) In
order to reduce estimation variance, Collaborative Trees make update decisions
based on probability weights derived from split scores, which are related to
the split sampling decision for BART [12].

Furthermore, when required for inference applications, we group each
continuous variable, such as the temperature variable in the embryo growth
dataset, into n_bins bins of equal sizes, creating a feature group with one-
hot indicators representing data in each bin. The detail of hyperparameter
tuning and an algorithm runtime discussion are given in the Supplementary
Material. All source codes for this paper are available upon request.

##  3 Extended mean decrease in impurity

In this section, we introduce how to calculate the extended mean decrease in
impurity (XMDI) for each Collaborative Trees model grown in Algorithm 1,
assuming a centered response. To calculate the XMDI for Collaborative Trees
Ensemble (see Section 2.3), the values of XMDIs are the respective averages
over all bagging tree models.

At the end of the $s$th round of update (i.e., Step 13 of Algorithm 1) with
$s>0$, let $\widehat{U}_{s}(\vv{x})=\sum_{k=1}^{K}T_{k}(\vv{x})$ for each
$\vv{x}\in\mathbb{R}^{p}$, with $\widehat{U}_{0}(\vv{x})=0$, and let
$\widehat{Q}_{s}$ denote the set of associated nodes to be updated. These
associated nodes will be split on the $\widehat{m}_{s}$th feature group for
some $1\leq\widehat{m}_{s}\leq M$. Let $1\leq\widehat{e}_{s}\leq s-1$ denote
the index of round when the parent node of $\widehat{Q}_{s}$ is updated, in
which $\widehat{e}_{s}=\emptyset$ if $\widehat{Q}_{s}$ consists of a root
node. In addition, for each $s>0$, define $\widehat{u}_{s}$ such that
$\widehat{u}_{s}=\widehat{m}_{s}$ if $\widehat{e}_{s}=\emptyset$,
$\widehat{u}_{s}=\widehat{u}_{\widehat{e}_{s}}$ if
$\widehat{m}_{s}=\widehat{m}_{\widehat{e}_{s}}$, and
$\widehat{u}_{s}=\widehat{m}_{\widehat{e}_{s}}$ if
$\widehat{m}_{s}\not=\widehat{m}_{\widehat{e}_{s}}$. Then, the XMDI is defined
to be such that for every $1\leq l\leq k\leq M$,

$\begin{split}\textnormal{XMDI}_{lk}&=\sum_{s\in\\{q:(\widehat{m}_{q},\widehat{u}_{q})=(l,k)\textnormal{
or
}(\widehat{m}_{q},\widehat{u}_{q})=(k,l)\\}}\left[\left(\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s}(\boldsymbol{X}_{i}))^{2}\right]\frac{1}{n}.\end{split}$

The sample responses are centered before calculating the XMDI in practice.
Notably, the definition of XMDI also accounts for applications with continuous
and categorical variables, where consecutive splits on the same feature are
possible. Meanwhile, the overall feature importance is given by

$\textnormal{XMDI}_{m}=\sum_{1\leq l\leq M}\textnormal{XMDI}_{lm}$ (5)

for $m\in\\{1,\cdots,M\\}$. Let us demonstrate the calculation of the XMDI
with Example 1.

######  Example 1.

In an application involving binary features and feature groups, child nodes
are not split on the same feature group as their parent nodes to avoid trivial
splits yielding zero split scores. Therefore, the calculation of XMDI involves
two key steps: 1.) $\textnormal{XMDI}_{mm}$ aggregates all split scores
resulting from splitting the $m$th feature group at the root nodes. 2.)
$\textnormal{XMDI}_{lk}$ is the summation of split scores from splitting the
$k$th (or $l$th) feature group on certain nodes, provided their parent nodes
have been split on the $l$th (or $k$th) group. For a visual demonstration of
computing XMDIs in this context, refer to Figure 4.

We now introduce the additive effects (or first-order effects or main effects)
and interaction effects at the population level. The overall prediction power
of $\boldsymbol{X}_{\mathcal{X}(m)}$ can be measured by the variance of the
residual
$f(\boldsymbol{X})-\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(m)})$,
where we recall $\boldsymbol{X}_{-\mathcal{X}(m)}$ is a random vector
consisting of all features but those in $\mathcal{X}(m)$. The additive effect
of $\boldsymbol{X}_{\mathcal{X}(m)}$ is measured by
$\textnormal{Var}(g_{m}(\boldsymbol{X}))$, while the two-way interaction
effect of $\boldsymbol{X}_{\mathcal{X}(l)}$ and
$\boldsymbol{X}_{\mathcal{X}(k)}$ is measured by
$\textnormal{Var}(g_{lk}(\boldsymbol{X})-g_{l}(\boldsymbol{X})-g_{k}(\boldsymbol{X}))$.
Here,

$\begin{split}g_{J}(\boldsymbol{X})&=\mathbb{E}(f(\boldsymbol{X})-\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(J)})|\boldsymbol{X}_{\mathcal{X}(J)}),\end{split}$
(6)

for each $J\subset\\{1,\cdots,M\\}$, with $\mathcal{X}(J)=\cup_{m\in
J}\mathcal{X}(m)$. We subtract $g_{l}(\boldsymbol{X})$ and
$g_{k}(\boldsymbol{X})$ from $g_{lk}(\boldsymbol{X})$ to better focus on
estimating the two-way interaction effect between
$\boldsymbol{X}_{\mathcal{X}(l)}$ and $\boldsymbol{X}_{\mathcal{X}(k)}$. For
simplicity, we write $g_{J}(\boldsymbol{X})$ as $g_{m}(\boldsymbol{X})$ for
univariate $J=\\{m\\}$, and $g_{lk}(\boldsymbol{X})$ as
$g_{J}(\boldsymbol{X})$ for $J=\\{l,k\\}$.

An example of $f(\boldsymbol{X})$ consisting of linear and XOR interaction
components is given in Example 2 of Section 4.1. In Theorem 2, we demonstrate
that XMDImm and XMDIlk = XMDIkl are respectively consistent estimators of
$\textnormal{Var}(g_{m}(\boldsymbol{X}))$ and
$\textnormal{Var}(g_{lk}(\boldsymbol{X})-g_{l}(\boldsymbol{X})-g_{k}(\boldsymbol{X}))$
when feature groups are independent and regularity conditions hold. Therefore,
XMDIm defined as in (5) evaluates the sum of the additive effect of the $m$th
feature group and its interactions with the other $M-1$ feature groups.

Figure 4: In this graphical example with $2$ collaborative decision trees,
each set of associated nodes is highlighted in blue. At the $1$st round, the
root node of the left tree is split on group $\mathcal{X}(m_{1})$ with
$\texttt{\\#}\mathcal{X}(m_{1})=3$, resulting in three child nodes. At the
$4$th round, a set of associated nodes on the left tree is about to be updated
on some group $\mathcal{X}(m_{2})$ with $\texttt{\\#}\mathcal{X}(m_{2})=2$. In
addition, after each update, the split score is attributed to the
corresponding group feature importance, as illustrated in this graphical
example.

###  3.1 Comparison with Sobol’ indices and conventional two-way interaction
models

To illustrate with an example, let us consider $f(\boldsymbol{X})=\beta
X_{1}X_{2}+X_{1}+X_{2}$ with $\mathbb{E}(f(\boldsymbol{X}))=0$. In this
example, the conventional two-way interaction inference [6] aims at testing
whether $\beta=0$. Meanwhile, Sobol’ indices measure the additive effect of
$X_{j}$ by $\textnormal{Var}[\mathbb{E}(f(\boldsymbol{X})|X_{j})]$ for
$j\in\\{1,2\\}$, while the interaction effect between $X_{1}$ and $X_{2}$ is
measured by
$\textnormal{Var}[\mathbb{E}(f(\boldsymbol{X})|X_{1},X_{2})]-\textnormal{Var}[[\mathbb{E}(f(\boldsymbol{X})|X_{1})]-\textnormal{Var}[\mathbb{E}(f(\boldsymbol{X})|X_{2})]$.
As a result, when features are independent, there is no interaction effect if
and only if $\beta=0$, regardless of which of the three approaches mentioned
here is used. Among these approaches, a major restriction of the conventional
one is the assumption of a parametric model, compared to Sobol’ indices and
XMDIs.

On the other hand, the additive and interaction Sobol’ indices often exhibit a
bias towards irrelevant features that are correlated with relevant ones, and
hence may not be apt for statistical applications. For example, given
$f(\boldsymbol{X})=\gamma X_{1}+X_{3}$ with a sufficiently large $\gamma>0$,
it can be shown that
$\textnormal{Var}[\mathbb{E}(f(\boldsymbol{X})|X_{2})]>\textnormal{Var}[\mathbb{E}(f(\boldsymbol{X})|X_{3})]$
if $X_{1}$ and $X_{2}$ are highly correlated, and $X_{3}$ is independent of
$(X_{1},X_{2})$. In contrast, the population XMDI correctly identifies $X_{2}$
as an irrelevant feature, yielding $\textnormal{Var}(g_{2}(\boldsymbol{X}))=0$
in this example, regardless of the arbitrary correlation between $X_{1}$ and
$X_{2}$.

In fact, XMDIs are closely related to the total Sobol’ index [5], which
quantifies the importance of $X_{\mathcal{X}(m)}$ using
$\mathbb{E}[\textnormal{Var}(f(\boldsymbol{X}|\boldsymbol{X}_{-\mathcal{X}(m)})]=\mathbb{E}\\{[f(\boldsymbol{X})-\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(m)})]^{2}\\}$,
addressing the bias concern by considering the variance of the residual
function
$f(\boldsymbol{X})-\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(m)})$.
By way of comparison, $\textnormal{Var}(g_{m}(\boldsymbol{X}))$ given in (6)
quantifies the variance of the residual function
$f(\boldsymbol{X})-\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(m)})$
projected onto the $m$th feature group. Similarly, our measure of two-way
interaction is also based on decompositions of the residual functions.
Therefore, XMDIs may offer enhanced reliability in the presence of correlated
features compared to Sobol’ indices [35].

Our theory in Section 4.2 provides a solid guarantee that positive XMDI values
are exclusive to significant components asymptotically when feature groups are
independent. Meanwhile, our empirical analysis in Section 5.1 demonstrates the
advantages of XMDIs over Sobol’ indices in distinguishing significant
components from insignificant ones, particularly when features are highly
correlated. It is noteworthy that XMDIs may benefit from incorporating
debiasing techniques, as in [26, 41, 5, 37], although the details are beyond
the scope of this paper.

###  3.2 Network diagrams

Let us introduce the four key elements in network diagrams (see Figure 1) as
follows: 1.) The sizes of the circles represent the overall feature importance
$\textnormal{XMDI}_{i}$’s. 2.) The thickness of each edge depends on
$\textnormal{XMDI}_{ij}$. 3.) The colors (blue or red) of the features depend
on the ratio $\frac{\textnormal{XMDI}_{ii}}{\textnormal{XMDI}_{i}}$. A feature
with high value of this ratio is represented in a deep blue shade, indicating
that it is more likely to be an additive component in the model. Otherwise,
this feature tends to be an interaction components, since
$\textnormal{XMDI}_{i}=\sum_{1\leq l\leq M}\textnormal{XMDI}_{il}$. 4.) The
grayness of edges and the direction of arrows are determined by
$\frac{\textnormal{XMDI}_{ij}}{\textnormal{XMDI}_{i}}$, reflecting the
strength of interaction between $X_{i}$ and $X_{j}$ based on the overall
importance of $X_{i}$. The arrow points towards the feature in the
denominator. In each diagram, we are interested in features represented by
larger circles in a deep blue shade, and interactions depicted with greater
thickness and darker edges. We also present the maximum values of
$\frac{\textnormal{XMDI}_{ij}}{\textnormal{XMDI}_{i}}$ and
$\textnormal{XMDI}_{i}\times\left[\textnormal{Sample variance of the
response}\right]$. These values, expressed as percentages, are displayed atop
the corresponding measures in the diagrams.

##  4 Asymptotic analysis of Collaborative Trees

In this section, we show the sure screening property [14] of Collaborative
Trees in Algorithm 1 and consistent XMDI estimation under proper model
assumptions, allowing the number of important additive and interaction
components to increase polynomially with the sample size. In contrast, it is
well-known that decision trees cannot effectively capture smooth signals in
additive models [22, 39] when the number of significant components is of
polynomial order in the sample size.

A key observation here is that our tree model is a type of matching pursuit
[32] (see (7) in Section 4.1) within the first $2K$ updates. Therefore, from a
theoretical viewpoint, we are interested in some important topics of
asymptotic analysis for trees and/or matching pursuit, including sure
screening property [14], model consistency and convergence rates for matching
pursuit [4, 27, 25] and for trees [11, 24], and estimation consistency of
feature importance measures [34, 31]. In this paper, we formally analyze the
sure screening property of Collaborative Trees in Theorem 1 and the estimation
consistency of XMDIs in Theorem 2, while leaving the tree model consistency
for future work.

Let us summarize the analysis challenges for our tree model. We will
demonstrate that the first $K$ updates of Collaborative Trees are a standard
matching pursuit [32], while the subsequent $K$ updates involve tree heredity
(see Condition 3) and projections onto multiple feature groups. Hence, the
existing model consistency results [32, 4, 24, 27] developed for the standard
matching pursuit are not directly applicable to our tree model. In addition,
the updates after the first $2K$ steps in Collaborative Trees are more akin to
sums of trees [12] than matching pursuit, making the analysis more involved
and may further require analysis techniques for high-dimensional trees [11,
24]. Therefore, we shall focus on analyzing our tree inference (the sure
screening property and XMDIs consistency) and postpone the consistency
analysis of Collaborative Trees to future work.

On the other hand, sure screening at the early updates within the first $2K$
updates ensures reliable assessment of important features. The sure screening
property has never been formally established for trees, to the best of our
knowledge. Our matching pursuit with the non-standard projection scheme under
tree heredity makes a rigorous analysis of sure screening a challenging task.
Furthermore, the XMDI estimation of
$\textnormal{Var}(g_{J}(\boldsymbol{X}))$’s (see (6)) involves the analysis of
the prediction function
$\mathbb{E}(f(\boldsymbol{X})|\boldsymbol{X}_{-\mathcal{X}(J)})$, making our
technical analysis of the XMDI estimation consistency especially difficult
when considering correlation between feature groups. We overcome these
technical challenges by controlling the distributional dependence between
feature groups (see Condition 2 in Section 4.1) and assuming that the feature
selection path has no repeated elements (see Theorem 1). Despite having
theoretical limitations, our theory is among the first formal results that
establish equivalence between trees and matching pursuit, which has been
discussed and inspected by many authors [25, 40, 12].

###  4.1 Sure screening additive and interaction components

In this section, we show that Collaborative Trees introduced in Section 2.2
can identify all significant components in the data generating models. Since
individual binary feature can be seen as a group of a single feature, our
analysis is established without loss of generality for feature groups
$\mathcal{X}(1),\cdots,\mathcal{X}(M)$ defined in Section 2.1. The set of
significant additive components is given by
$S_{1}^{*}=\\{m:\textnormal{Var}(g_{m}(\boldsymbol{X}))>0,1\leq m\leq M\\}$,
while the set of significant interaction components is defined to be
$S_{2}^{*}=\\{\\{l,k\\}:\textnormal{Var}(g_{lk}(\boldsymbol{X})-g_{l}(\boldsymbol{X})-g_{k}(\boldsymbol{X}))>0,1\leq
l<k\leq M\\}$. The function $g_{J}(\boldsymbol{X})$’s are defined in (6).

Our asymptotic analysis is based on a crucial observation: When features are
either binary or groups of one-hot indicators, Collaborative Trees exhibit
asymptotic similarity to the matching pursuit method. To explain this
connection, we begin with a recap of Algorithm 1 from Section 2.2, and provide
some definitions of notation. In each round of tree update, a set of
associated nodes is updated together w.r.t. some group $\mathcal{X}(m)$.
Consider the $s$th round of update on some set of associated nodes in the
$k$th tree. Let $\widehat{J}_{s}\subset\\{1,\cdots,M\\}$ denote the set of
group indices representing the update history from the root node of the $k$th
tree to the current set of associated nodes. Particularly, the index of the
feature group to be split at the $s$th round of update is also in
$\widehat{J}_{s}$. For example,
$\widehat{J}_{1}=\\{m_{1}\\},\widehat{J}_{2}=\\{m_{3}\\},\widehat{J}_{3}=\\{m_{3},m_{1}\\}$,
and $\widehat{J}_{4}=\\{m_{1},m_{2}\\}$ in Figure 4.

To gain insights into the asymptotic equivalence between Collaborative Trees
and matching pursuit, we consider functions
$R_{l}^{\star}:\mathbb{R}^{p}\mapsto\mathbb{R}$ defined recursively for
$0<l\leq 2K$:

$R_{l}^{\star}(\boldsymbol{X})=R_{l-1}^{\star}(\boldsymbol{X})+\mathbb{E}(f(\boldsymbol{X})-R_{l-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(\widehat{J}_{l})}),$
(7)

where $R_{0}^{\star}(\vv{x})=0$ and $\mathcal{X}(J)=\cup_{m\in
J}\mathcal{X}(m)$ for $J\subset\\{1,\cdots,M\\}$. In our technical proofs, we
show that setting $J=\widehat{J}_{s}$ maximizes the sample counterpart of the
function
$L(R_{s-1}^{\star},J)=\textnormal{Var}\left[\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)})\right]$,
subject to $J\in\Theta_{s}$ for some constraints $\Theta_{s}$’s. These
constraints are such that
$\Theta_{1}=\cdots=\Theta_{K}=\\{\\{1\\},\cdots,\\{M\\}\\}$, and
$\Theta_{K+1},\cdots,\Theta_{2K}$ are related to tree heredity. Furthermore,
$L(R_{s-1}^{\star},\widehat{J}_{s})$ is asymptotically equivalent to the
sample split score at the $s$th split. In essence, given our tree model
regularization (2), the first $K$ updates of Collaborative Trees follow a
standard matching pursuit [32], while the subsequent $K$ updates involve tree
heredity and projections onto multiple feature groups. In Theorem 1, our
analysis relies on this connection to establish results ensuring the sure
screening property of $\widehat{J}_{1},\cdots,\widehat{J}_{2K}$, as stated in
the same theorem. The following technical conditions are required for Theorem
1.

######  Condition 1.

Assume $\boldsymbol{X}$ is a binary random vector. If
$\texttt{\\#}\mathcal{X}(m)>1$ for some $1\leq m\leq M$, then
$\sum_{j\in\mathcal{X}(m)}\boldsymbol{1}_{X_{j}=1}=1$ and
$\texttt{\\#}\mathcal{X}(m)\leq M_{\mathcal{X}}$ for some constant $2\leq
M_{\mathcal{X}}<\infty$. Let $Y=f(\boldsymbol{X})+\varepsilon$ with
$\mathbb{E}(f(\boldsymbol{X}))=0$,
$\max_{\vv{c}\in\\{0,1\\}^{p}}|f(\vv{c})|\leq M_{f}$ for some constant
$1<M_{f}<\infty$, and zero mean model error $\mathbb{E}(\varepsilon)=0$ with
$\mathbb{E}(\varepsilon\boldsymbol{1}\\{|\varepsilon|\leq
c\\}|\boldsymbol{X})=0$ for each $c>0$ and
$\mathbb{E}|\varepsilon|^{q_{0}}<M_{\varepsilon}$ for some constants $q_{0}>0$
and $M_{\varepsilon}<\infty$.

######  Condition 2.

Assume
$\textnormal{Var}(g_{m}(\boldsymbol{X}))>2d_{n}+42\delta_{0}p_{\min}^{-2}\mathbb{E}[(f(\boldsymbol{X}))^{2}]$
and $\textnormal{Var}(g_{J}(\boldsymbol{X})-(\sum_{j\in
J}g_{j}(\boldsymbol{X})))>2d_{n}+7616\delta_{0}(\texttt{\\#}S_{2}^{*})^{4}p_{\min}^{-2}\mathbb{E}[(f(\boldsymbol{X}))^{2}]$
for some $d_{n}>0$, each $m\in S_{1}^{*}$, and each $J\in S_{2}^{*}$, where
$p_{\min}=\min_{1\leq l<k\leq
M,i\in\mathcal{X}(l),j\in\mathcal{X}(k),(a,b)\in\\{0,1\\}^{2}}\mathbb{P}(X_{i}=a,X_{j}=b)$,
and

$\begin{split}\delta_{0}&=\max_{1\leq m\leq M}\inf\left\\{\delta\geq
0:\mathbb{P}\left(\max_{i\in\mathcal{X}(m)}|\mathbb{P}(X_{i}=1|\boldsymbol{X}_{-\mathcal{X}(m)})-\mathbb{P}(X_{i}=1)|\leq\delta\right)=1\right\\}.\end{split}$

######  Condition 3 (Tree heredity).

For every $S^{\dagger}\subset S_{1}^{*}$ and $J\in S_{2}^{*}$ such that each
element in $S^{\dagger}$ belongs to a distinct set in $(S_{2}^{*}\backslash
J)$, it holds that $J\cap(S_{1}^{*}\backslash S^{\dagger})\not=\emptyset$.

The regularity conditions on $f(\boldsymbol{X})$ and model error in Condition
1 are assumed to simplify the application of standard concentration
inequalities. No specific heteroskedastic assumption on
$\textnormal{Var}(\varepsilon|\boldsymbol{X})$ is required. The condition
$p_{\min}\geq p_{01}>0$ assumed in Theorem 1 holds if feature groups are
independent and that $\mathbb{P}(X_{j}=1)$ is bounded away from one and zero
for all binary features. Additionally, when all components in $S_{1}^{*}$ are
selected, Condition 3 ensures tree heredity, allowing for the selection of all
components in $S_{2}^{*}$. For example, the case where $S_{1}^{*}=\\{2,4\\}$
and $S_{2}^{*}=\\{\\{2,3\\},\\{2,4\\}\\}$ is not permitted, whereas the case
where $S_{1}^{*}=\\{3,4\\}$ and $S_{2}^{*}=\\{\\{2,3\\},\\{2,4\\}\\}$
satisfies Condition 3. Hence, Condition 3 is more restrictive than weak
heredity [6], but it is natural for sum of trees.

Condition 2 requires the signal strength of significant components to surpass
the sum of $2d_{n}$, representing estimation variance, and an additional term
capturing bias due to feature group distributional dependence, measured by
$\delta_{0}$. This condition ensures that all important components are
sufficiently strong and can be effectively screened by Collaborative Trees, as
stated in Theorem 1. The technical bias terms may benefit from improvement
through a more careful analysis.

In Theorem 1, we assume that all nodes in the initial $2K$ updates are valid
in Algorithm 1, and define
$\mathbb{E}(h(\boldsymbol{X})|\boldsymbol{X}=\vv{c})=0$ when
$\mathbb{P}(\boldsymbol{X}=\vv{c})=0$ for every $\vv{c}\in\\{0,1\\}^{p}$ and
every measurable function $h:\mathbb{R}^{p}\mapsto\mathbb{R}$. Recall
$R_{l}^{\star}(\boldsymbol{X})$ and $\widehat{J}_{l}$ from equation (7). The
technical proofs are provided in the Supplementary Material.

######  Theorem 1.

Assume $\widehat{J}_{l}\not=\widehat{J}_{k}$ for each
$\\{l,k\\}\subset\\{1,\cdots,2K\\}$, Condition 1, and that
$(\boldsymbol{X}_{1},Y_{1}),\cdots,(\boldsymbol{X}_{n},Y_{n}),(\boldsymbol{X},Y)$
are i.i.d. Assume Conditions 2–3 hold with $p_{\min}\geq
p_{01}\vee\frac{2\log{(n)}}{\sqrt{n}}$, and
$d_{n}=L_{1}n^{-\frac{1}{2}}(\log{n})(n^{\frac{1}{q_{1}}}\sqrt{2K}+\sqrt{\sum_{l=1}^{2K}\iota_{l-1}^{2}})(\iota_{2K-1}^{2}+\iota_{2K-1}+3)$,
in which $0<p_{01}<1$ is a constant, $L_{1}>0$ is a sufficiently large
constant, $q_{1}$ is an arbitrary constant with $q_{0}>q_{1}>0$, and that
$\iota_{-1}=\iota_{0}=0$ and

$\begin{split}\iota_{s}&=\sum_{l=1}^{s}\max_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(\widehat{J}_{l})}}\left|\mathbb{E}[(f(\boldsymbol{X})-R_{l-1}^{\star}(\boldsymbol{X}))|\boldsymbol{X}_{\mathcal{X}(\widehat{J}_{l})}=\vv{c}]\right|\qquad\textnormal{
for }s>0.\end{split}$ (8)

Assume $M\geq K\geq\texttt{\\#}S_{1}^{*}\vee\texttt{\\#}S_{2}^{*}$, $d_{n}\leq
1$, and that $M=O(n^{K_{01}})$ for some arbitrary constant $K_{01}>0$. Then,
$\mathbb{P}(\\{S_{1}^{*}\not\subset\cup_{l=1}^{K}\widehat{J}_{l}\\}\cup(S_{2}^{*}\not\subset\\{\widehat{J}_{K+1},\cdots,\widehat{J}_{K+\texttt{\\#}S_{2}^{*}}\\}))=o(1)$.

Theorem 1 investigates Collaborative Trees’ capability to select numerous
significant components from high-dimensional input features; note that the
tree model is not necessarily consistent after sure screening important
components. According to the current theory, a necessary condition for the
estimation error $d_{n}$ to diminish asymptotically is $K\leq
n^{\frac{1}{7}}$, which is derived from Supplementary Material, where we show
that $\iota_{s}\leq sc_{0}$ for some constant $c_{0}>0$. Corollary 1 below
relaxes the requirement on $K$ by assuming independent features, although the
results still do not achieve optimality in the number of significant
components allowed in the model. For example, the Lasso allows
$o(\sqrt{n/\log{p}})$ important features given proper regularity conditions
[8]; this rate may be further improved under stronger conditions. Despite
technical limitations that prevent optimal results, Theorem 1 demonstrates
Collaborative Trees’ ability to handle well beyond $\log{(n)}$ significant
components, which is advantageous than decision trees, as have introduced at
the beginning of Section 4. Such a result is further illustrated in Corollary
1 below, where notation follows that in Theorem 1, while
$S_{1}\subset\\{1,\cdots,M\\}$ and $S_{2}$ takes the form
$\\{\\{m_{1},m_{2}\\},\\{m_{3},m_{4}\\},\cdots\\}$.

######  Corollary 1.

Let $M_{f}>1$ be an arbitrary constant, and $S_{1}$ and $S_{2}$ be two sets
satisfying Condition 3 with $\texttt{\\#}S_{1}\vee\texttt{\\#}S_{2}\leq
K\leq(\log{n})^{-1}n^{\frac{1}{5}-\frac{2}{5q_{1}}}$. Let
$f(\boldsymbol{X})=\sum_{m\in
S_{1}}h_{m}(\boldsymbol{X}_{\mathcal{X}(m)})+\sum_{J\in
S_{2}}h_{J}(\boldsymbol{X}_{\mathcal{X}(J)})$ with $\delta_{0}=0$ (i.e.,
independent feature groups), in which $\boldsymbol{X}$ satisfies Condition 1
with $p_{\min}\geq p_{01}$,
$\mathbb{E}(h_{m}(\boldsymbol{X}_{\mathcal{X}(m)}))=0$ and
$\mathbb{E}(h_{J}(\boldsymbol{X}_{\mathcal{X}(J)})|\boldsymbol{X}_{\mathcal{X}(l)})=0$
for each $m\in S_{1}$ and each $l\in J\in S_{2}$, and that
$\textnormal{Var}(h_{J}(\boldsymbol{X}_{\mathcal{X}(J)}))=16L_{1}n^{-\frac{1}{2}}(\log{n})n^{\frac{1}{q_{1}}}\sqrt{2K}\eqqcolon
w_{n}$ for each $J\in\\{\\{m\\}:m\in S_{1}\\}\cup S_{2}$. Then,
$o(1)=w_{n}\geq 2d_{n}$ and $\max_{\vv{c}\in\\{0,1\\}^{p}}|f(\vv{c})|<M_{f}$
for all large $n$. In addition, if other required conditions for Theorem 1 are
assumed, then
$\mathbb{P}(\\{S_{1}\not\subset\cup_{l=1}^{K}\widehat{J}_{l}\\}\cup(S_{2}\not\subset\\{\widehat{J}_{K+1},\cdots,\widehat{J}_{K+\texttt{\\#}S_{2}}\\}))=o(1)$.

Given the assumption of a bounded regression function, the variance of each
component $w_{n}$ has to decrease asymptotically to accommodate an increasing
(polynomially with respect to $n$) $\texttt{\\#}S_{1}\vee\texttt{\\#}S_{2}$,
as stated in Corollary 1. Example 2 below is an instance of
$f(\boldsymbol{X})$ in Corollary 1, where $S_{1}$ and $S_{2}$ respectively
consist of linear components and XOR interaction components.

######  Example 2.

Assume each feature group has a single feature with $M=p$. Let
$f(\boldsymbol{X})=\big{[}\sum_{m\in
S_{1}}\beta_{m}(X_{m}-\mathbb{E}(X_{m}))\big{]}+\sum_{\\{l,k\\}\in
S_{2}}[\beta_{lk}\boldsymbol{1}_{(X_{l}-\mathbb{E}(X_{l}))(X_{k}-\mathbb{E}(X_{k}))>0}-\beta_{lk}\boldsymbol{1}_{(X_{l}-\mathbb{E}(X_{l}))(X_{k}-\mathbb{E}(X_{k}))<0}]$,
in which $\beta_{m}$’s and $\beta_{lk}$’s are real coefficients, $X_{j}$’s are
i.i.d. Bernoulli random variables.

When dealing with non-identically distributed features, groups with more than
one feature, and a small $\delta_{0}>0$, one cannot anticipate a clean model
setting as in Example 2. Nevertheless, Theorem 1 guarantees Collaborative
Trees identifying all important components, whose count potentially grows
polynomially with sample size given regularity conditions and small
$\delta_{0}$, in high-dimensional nonlinear settings.

###  4.2 Analysis of additive and interaction effects

In this section, we analyze the XMDIs defined in Section 3. Recall
$\sum_{k=1}^{K}T_{k}(\boldsymbol{X})$ from Algorithm 1,
$R_{l}^{\star}(\boldsymbol{X})$ from (7), and $g_{J}(\boldsymbol{X})$ from
Section 3. The technical proof of Theorem 2 is provided in Section B of the
Supplementary Material.

######  Theorem 2.

Assume all conditions in Theorem 1. Then, the following results hold with
probability approaching one. For each feature group $m\in\\{1,\cdots,M\\}$,

$\begin{split}|\textnormal{XMDI}_{mm}-\textnormal{Var}(g_{m}(\boldsymbol{X}))|&\leq
d_{n}+L_{3}\delta_{0},\\\
\sum_{k\not=m}|\textnormal{XMDI}_{mk}-\textnormal{Var}(g_{mk}(\boldsymbol{X})-g_{m}(\boldsymbol{X})-g_{k}(\boldsymbol{X}))|&\leq\Delta+\widehat{\Gamma}+(\texttt{\\#}S_{2}^{*}+1)d_{n}+L_{3}\delta_{0}(\texttt{\\#}S_{2}^{*})^{5},\end{split}$

where $d_{n}$ and $\delta_{0}$ are given in Theorem 1, $L_{3}>0$ is a
sufficiently large constant,
$\Delta=\textnormal{Var}(f(\boldsymbol{X})-R_{K+\texttt{\\#}S_{2}^{*}}^{\star}(\boldsymbol{X}))$,
and
$\widehat{\Gamma}=(n^{-1}\sum_{i=1}^{n}\varepsilon_{i}^{2})-n^{-1}\sum_{i=1}^{n}(Y_{i}-\sum_{k=1}^{K}T_{k}(\boldsymbol{X}_{i}))^{2}$.

In the context of modern MDI feature importance analysis [34, 31, 26], Theorem
2 is the first to formally separate feature importance of trees into additive
and interaction effects. Our results ensure the reliable inference when
feature dependence is mild and the number of significant interaction
components is limited. In Theorem 2, the second upper bound involves
$\texttt{\\#}S_{2}^{*}\times d_{n}$ since we sum up all estimation errors
resulting from interaction estimation. The term $(\texttt{\\#}S_{2}^{*})^{5}$
may be improved through a more careful analysis.

######  Remark 1.

As mentioned at the beginning of Section 4, the analysis of $\Delta$ and
$\widehat{\Gamma}$, which is related to tree model consistency, could be too
involved for the current paper. Instead, we empirically assess the consistency
quality of our tree model in Section 6.2. On the other hand, it is well-
documented that impurity-based measures exhibit bias towards irrelevant
features correlated with relevant ones [31, 36]. Despite ongoing efforts to
mitigate this bias [26, 29, 41, 5, 10], a formal analysis of bias in feature
importance derived from sum-of-trees models is still lacking. While our
Theorem 2 is limited in assessing the bias issue due to the stringent
technical requirements on $\delta_{0}$, we provide empirical experiments in
Section 5.1 demonstrating that the XMDIs are bias-resistant.

##  5 Simulation experiments

This section demonstrates XMDIs’ numerical stability and resistance to bias
from irrelevant correlated features. We also investigate key hyperparameters’
impact on Collaborative Trees Ensemble’s inferences using artificial datasets
generated from Model (9), where each feature group contains a single feature.
Groups are therefore ignored for simplicity.

$\displaystyle Y$
$\displaystyle=5X_{1}+20(X_{3}-0.5)^{2}+15X_{5}+2X_{9}+10\sin{(\pi(X_{9}-0.5)(X_{10}-0.5))}+\varepsilon,$
(9)

where $\varepsilon$ is an independent standard Gaussian model error. The
covariate distribution of $(X_{1},\cdots,X_{p})^{T}$ is chosen as a Gaussian
copula, ensuring that $0\leq X_{j}\leq 1$ almost surely, with an AR(1)
covariance matrix $\Sigma=(\lambda^{|k-l|})$ for each
$\\{k,l\\}\subset\\{1,\cdots,p\\}$, in which $0\leq\lambda<1$ and
$p\in\\{10,100\\}$. Specifically, $(X_{1},\cdots,X_{p})$ has a joint CDF
function $\Phi_{\Sigma}(\Phi^{-1}(x_{1}),\cdots,\Phi^{-1}(x_{p}))$, where
$\Phi_{\Sigma}$ is the CDF of multivariate Gaussian with mean zero and
covariance matrix $\Sigma$, and $\Phi^{-1}$ is the inverse of the univariate
standard Gaussian CDF. Model (9) is a variant of the Friedman regression
function [16], which is used for testing the selection power for the linear
main effect, quadratic effect, and interaction effect. A weak additive
component $2X_{9}$ is included in the model to ensure the tree heredity and
stabilize the numerical results.

The hyperparameter tuning of Collaborative Trees Ensemble (see Section 2.3) is
performed as introduced in the Supplementary Material, with $20$ runs of
optimization. Some important tuned hyperparameters are $K=12$ and
$\textsf{n\\_bins}\in\\{None,5\\}$ for the simulation experiments in Section
5. Here, the original features $X_{j}$’s are used when
$\textsf{n\\_bins}=None$. Otherwise, every continuous feature is binned into
n_bins bins, and becomes a feature group. We set the sample size $n=500$ for
all simulation experiments.

###  5.1 XMDIs are bias-resistant and numerically stable

In this section, we empirically demonstrate that the XMDIs of Collaborative
Trees Ensemble can effectively separate significant and insignificant
components in learning environments with high correlation between features. To
illustrate the bias resistance of XMDIs, we conduct simulations based on Model
(9) with $\lambda\in\\{0.1,0.8\\}$ and $p=100$. The results, presented in
Table 1, indicate the rates at which a threshold exists to distinguish the
importance measures of significant and insignificant components in Model (9).
The significant features are defined within the index set
$S=\\{1,3,5,9,10\\}$, where four additive components are respectively based on
$(X_{1},X_{3},X_{5},X_{9})$ and an interaction component is based on
$(X_{9},X_{10})$. Meanwhile, the noisy features are in
$S^{c}\coloneqq\\{1,\cdots,100\\}\backslash S$. It is important to note that
our focus here is not on estimating the separation threshold. We compare the
separation rates obtained from 100 independent and repeated simulation
experiments with those obtained using Random Forests MDI, Sobol’ indices [35,
33], and conditional permutation importance (CPI [36]). Among many methods [5,
26, 29, 41, 10], CPI is a well-established and representative measure proposed
for addressing the bias issue of Random Forests MDI. The total (overall),
first-order (additive), and second-order (interaction) Sobol’ indices [35, 33]
are calculated with Random Forests as the prediction model. See the
Supplementary Material for additional details of Sobol’ indices and
hyperparameter tuning.

In column (I) of Table 1, XMDI and CPI consistently distinguish between
important features and noisy variables across 100 experiments, showcasing
their numerical stability and robustness in the presence of highly correlated
features ($\lambda=0.8$) and numerous noisy variables ($p=100$). The column
(III) shows that the estimated additive effects of $X_{1}$, $X_{3}$, and
$X_{5}$ based on the XMDI indeed consistently dominate those from
insignificant features in $S^{c}$. However, in addition to results in Table 1,
we report here that the additive effect of $X_{9}$ is too weak so the
separation rates of $X_{9}$ to $S^{c}$ based on the XMDI are roughly 80% for
$\lambda\in\\{0.1,0.8\\}$; more samples and a stronger signal can increase
these rates. Moving forward to column (II), the interaction effect between
$X_{9}$ and $X_{10}$ is accurately identified by the XMDI, further
underscoring its precision in capturing meaningful interactions amidst complex
data scenarios. On the other hand, even with mild correlation $\lambda=0.1$,
Sobol’ indices struggle to effectively distinguish the interaction effect
between $X_{9}$ and $X_{10}$. The scenario with $\lambda=0.8$ illustrates that
Sobol’ indices exhibit some resistance to bias, except when it comes to
selecting interaction effects. It is worth noting that the inferior
performance of MDI in the high correlation scenario underscores the
considerable challenges posed by the learning environment for the methods
under investigation, thereby emphasizing the bias-resistant nature of XMDI and
CPI, despite the fact that CPI does not apply in cases (II) and (III). These
results reinforce the reliability and informativeness of our network diagrams
for drawing insightful inferences.

Method | $\lambda=0.1$ | $\lambda=0.8$
---|---|---
(I) | (II) | (III) | (I) | (II) | (III)
XMDI of Collaborative Trees Ensemble | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Sobol’ indices [35] | 0.93 | 0.10 | 0.83 | 0.85 | 0.00 | 0.70
MDI [7] of Random Forests | 1.00 | Not Applicable | 0.37 | Not Applicable
CPI [36] of Random Forests | 1.00 | Not Applicable | 0.98 | Not Applicable
Table 1: (I) Success rates of distinguishing significant and insignificant
features in Model (9) based on overall feature importance measures. (II) Rates
of the interaction effect between $X_{9}$ and $X_{10}$ being the maximum among
all pairs. (III) Rates of the additive effects of $(X_{1},X_{3},X_{5})$
dominating those from insignificant features in $S^{c}$. All rates are
obtained from 100 simulation experiments with $p=100$ features and $n=500$.
MDI and CPI do not infer feature interaction or additive effects.

###  5.2 Number of trees $K$, data binning, and feature correlation

Figure 5 presents network diagrams (see Section 3.2) for data simulated from
Model (9) with $p=10$, in which other varying parameters are indicated in the
corresponding panels. From Figure 5, XMDI successfully distinguishes
significant features from noisy ones in Model (9) across all cases, except for
panel (d) with $K=1$. In this particular case, the importance of $X_{9}$ and
$X_{10}$ might be underestimated, and $X_{3}$ and $X_{1}$ are incorrectly
identified as interaction components. The results in panel (d) underscore the
crucial role of $K$ in our tree model for making accurate inferences,
consistent with our theoretical expectations. Meanwhile, the major distinction
of the other three panels lies in the inference about the interaction
components based on $X_{9}$ and $X_{10}$. A comparison between panels (a) and
(b) shows that whether to use binned features has only a mild impact on the
inferences of our tree models in this experiment. This result suggests the
inferences of the XMDI may still be informative without binning.

The adverse effects of feature correlation on our inferences are clearly seen
in panel (c). Notably, the importance of features $X_{9}$ and $X_{10}$ appears
slightly underestimated in this scenario compared to panels (a) and (b),
despite the excellent separation rates reported in Table 1 (the same setup but
with $p=10$ features here). Interestingly, the interaction effect account for
roughly only $47\%$ for their overall feature importance, which is far lower
than panels (a) and (b). The high feature correlation between $X_{9}$ and
$X_{10}$ weakens the contribution of their interaction effects to their
overall importance. This observation highlights the nuanced nature of studying
interaction effects, emphasizing the importance of understanding both the
regression function and feature correlation, namely, the joint distribution of
$(Y,\boldsymbol{X})$. Collaborative Trees represent our attempt to address
these complexities. Our finding echoes the empirical observations made by [12]
regarding the importance of a sum of trees for BART. However, their simulation
study does not investigate complex feature effects or the potential bias of
effect measures.

These simulation experiments enhance our understanding of the proposed tree
model, and raise intriguing research questions. Key areas for further
exploration include developing more efficient methods for binning features,
exploring the feasibility of making valid inferences from general continuous
features without binning, and advancing the XMDI to better resist bias by
incorporating debiasing techniques, such as CPI [36] and related works [5,
26].

(a)

(b)

(c)

(d)

Figure 5: Network diagrams for data simulated from Model (9) with $p=10$,
$n=500$, and other varying parameters indicated in the corresponding panels.

##  6 Real data inferences and prediction evaluation

###  6.1 Embryo growth data study

Temperature-dependent sex determination, or TSD, is a common occurrence among
vertebrates, including reptiles such as crocodilians, turtles, fish, and some
lizards. TSD is crucial for species resilience to climate change impacts. Our
analysis focuses on TSD using the“2022-03-29” version of the embryo growth
dataset comprising $874$ samples (with missing values in sex-ratio removed)
from the R package embryogrowth [19], covering data from 60 species, primarily
reptiles.

Each hatch of eggs comprises an observation in the sample. For every hatch of
eggs incubated in temperature-regulated chambers, the following information
has been recorded: the sex ratio of females to total births in this hatch
(response), species category (species and subspecies), the location where the
eggs come from (area and RMU, short for Regional Management Unit), incubation
period (incubation periods (days)), incubation temperature (temperature), and
incubation temperature amplitude (amplitude). The temperature amplitude is the
maximum variation of the incubation temperature, while temperatures are
recorded in Celsius. Other details of the dataset can be found in [19, 1].
There are $60$ species, $10$ subspecies, $81$ areas, and $14$ RMUs, which are
all grouped into feature groups with one-hot indicators. We apply
Collaborative Trees Ensemble based on these features for analyzing the sex
ratio. The resulting network diagram is shown in the left panel of Figure 1.
Hyperparameter tuning involves 100 optimization runs, yielding $K=11$ and
$\textsf{n\\_bins}=40$, where continuous features are binned into n_bins bins.

Different species are known to exhibit various patterns of TSD [28], implying
an interaction between species and temperature (each TSD pattern naturally
depends on temperature). Our XMDI inferences correctly identify the
interaction between species and temperature, as illustrated in the left
diagram of Figure 1. Additionally, the diagram reveals a distinct additive
effect associated with temperature, suggesting the prevalence of a specific
TSD pattern in the dataset. This observation aligns with our experience, as
the logistic-shaped TSD pattern, depicted in the right panel of Figure 1, is
frequently observed in the dataset. Now, let us investigate the importance of
the variable area. Considering that same species may inhabit same areas, a
plausible explanation for the clear interaction effect (the second strongest)
between area and temperature is the potential synonymy of area with species.
The distance correlation [38] between area and species is as high as $0.67$,
providing support for our explanation. Despite this, the overall importance of
area and its interaction between temperature are weaker compared to species.
This could be attributed to the presence of multiple distinct species in the
same region, diminishing the importance of area relative to species.

On the other hand, the diagram in Figure 1 suggests that the effects of
incubation periods (days) are not related to species, as the statistical
association pattern for incubation periods (days) is markedly different from
that of species. While we acknowledge that we are not experts in embryo growth
analysis, the variable incubation periods (days) has less often been
considered an important factor in TSD studies, making our results potentially
interesting for further investigation. Now, let us consider two potential
explanations: 1.) Uncontrolled environmental factors, such as humidity, might
influence incubation periods (days) and, consequently, the sex ratio. A more
detailed examination of the samples is essential for a thorough understanding.
2.) The detected signal of incubation periods (days) in Figure 1 may be
considered noise. Significance tests are required for further analysis.

Using the embryo growth dataset, we have demonstrated how to leverage
Collaborative Trees Ensemble for advanced data analysis. We remark that the
technical tree heredity requirement Condition 3 does not prevent the XMDI from
correctly accounting for multiple interactions stemming from the same feature
group in practice. A discussion, along with a larger version of Figure 1 and
detailed numerical results, is provided in the Supplementary Material to save
space.

###  6.2 Prediction accuracy evaluation

We evaluate the predictive performance of Collaborative Trees Ensemble and
compare it with Random Forests [7] and XGBoost [9]. These models are known for
their competitiveness when compared to modern deep learning models, as noted
in [20]. The benchmark comparison follows the methodology outlined in [20], a
well-established work that provides a state-of-the-art benchmark procedure and
simulation environments.

Following the approach outlined in [20], we utilize 19 selective datasets from
OpenML, each containing no more than $10,000$ samples, for prediction
evaluation based on the goodness-of-fit measure R2 value. The calculation of
each R2 value involves splitting the full sample into training, validation,
and test sets, constituting 48%, 32%, and 20% of the full dataset,
respectively. This division is employed for hyperparameter tuning based on the
training and validation samples, as well as prediction loss calculation based
on the test sample. We perform 100 optimization runs for the hyperparameter
tuning procedure, without binning any features (i.e., n_bins = $None$). To
ensure result stability, we conducted each experiment on the 19 datasets
independently 10 times.

To summarize the comparison, we calculated the average (adjusted) win rates
[20] for Collaborative Trees Ensemble, XGBoost, and Random Forests, which are
0.8069, 0.5986, 0.2439, respectively. While XGBoost has been recognized for
its strong performance against modern deep learning models in tabular
regression [20], our tree model surpassed its performance, highlighting the
competitive predictive capabilities of our proposed approach. For the detailed
R2 scores for each of the 19 datasets, along with other experiment specifics,
please refer to the Supplementary Material.

##  7 Discussion

We introduce Collaborative Trees, a novel tree model designed for regression
prediction, along with its bagging version, to analyze complex statistical
associations between features and uncover potential patterns in the data. In
conclusion, we highlight three challenges related to Collaborative Trees.
Firstly, the study of bias in XMDI and rigorous significance tests are
relevant for making more stable inferences. Moving forward, applying our
analysis tools to scenarios with lots of features, like in genetic data
analysis, poses a significant challenge. Lastly, investigating the consistency
of Collaborative Trees and its bagging version is crucial for the applications
of Collaborative Trees.

## References

  * Abreu-Grobois et al. [2020] Abreu-Grobois, F. A., B. A. Morales-Mérida, C. E. Hart, J.-M. Guillon, M. H. Godfrey, E. Navarro, and M. Girondot (2020).  Recent advances on the estimation of the thermal reaction norm for sex ratios.  PeerJ 8, e8451.
  * Antoniadis et al. [2021] Antoniadis, A., S. Lambert-Lacroix, and J.-M. Poggi (2021).  Random forests for global sensitivity analysis: A selective review.  Reliability Engineering & System Safety 206, 107312.
  * Apley and Zhu [2020] Apley, D. W. and J. Zhu (2020).  Visualizing the effects of predictor variables in black box supervised learning models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 82(4), 1059–1086.
  * Barron et al. [2008] Barron, A. R., A. Cohen, W. Dahmen, and R. A. DeVore (2008).  Approximation and learning by greedy algorithms.
  * Bénard et al. [2022] Bénard, C., S. Da Veiga, and E. Scornet (2022).  Mean decrease accuracy for random forests: inconsistency, and a practical solution via the sobol-mda.  Biometrika 109(4), 881–900.
  * Bien et al. [2013] Bien, J., J. Taylor, and R. Tibshirani (2013).  A lasso for hierarchical interactions.  Annals of statistics 41(3), 1111.
  * Breiman [2001] Breiman, L. (2001).  Random forests.  Machine learning 45(1), 5–32.
  * Bühlmann and Van De Geer [2011] Bühlmann, P. and S. Van De Geer (2011).  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * Chen and Guestrin [2016] Chen, T. and C. Guestrin (2016).  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pp. 785–794.
  * Chi et al. [2022] Chi, C.-M., Y. Fan, and J. Lv (2022).  Fact: High-dimensional random forests inference.  arXiv preprint arXiv:2207.01678.
  * Chi et al. [2022] Chi, C.-M., P. Vossler, Y. Fan, and J. Lv (2022).  Asymptotic properties of high-dimensional random forests.  The Annals of Statistics 50(6), 3415–3438.
  * Chipman et al. [2010] Chipman, H. A., E. I. George, and R. E. McCulloch (2010).  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics 4(1), 266–298.
  * Davis et al. [1997] Davis, G., S. Mallat, and M. Avellaneda (1997).  Adaptive greedy approximations.  Constructive approximation 13, 57–98.
  * Fan and Lv [2008] Fan, J. and J. Lv (2008).  Sure independence screening for ultrahigh dimensional feature space.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(5), 849–911.
  * Friedberg et al. [2020] Friedberg, R., J. Tibshirani, S. Athey, and S. Wager (2020).  Local linear forests.  Journal of Computational and Graphical Statistics 30(2), 503–517.
  * Friedman [1991] Friedman, J. H. (1991).  Multivariate adaptive regression splines.  The annals of statistics, 1–67.
  * Friedman [2001] Friedman, J. H. (2001).  Greedy function approximation: a gradient boosting machine.  Annals of statistics, 1189–1232.
  * Friedman and Meulman [2003] Friedman, J. H. and J. J. Meulman (2003).  Multiple additive regression trees with application in epidemiology.  Statistics in medicine 22(9), 1365–1381.
  * Girondot [2019] Girondot, M. (2019).  embryogrowth: Tools to analyze the thermal reaction norm of embryo growth.  The Comprehensive R Archive Network.
  * Grinsztajn et al. [2022] Grinsztajn, L., E. Oyallon, and G. Varoquaux (2022).  Why do tree-based models still outperform deep learning on typical tabular data?  Advances in neural information processing systems 35, 507–520.
  * Gwanyama [2004] Gwanyama, P. W. (2004).  The hm-gm-am-qm inequalities.  College Mathematics Journal, 47–50.
  * Hastie et al. [2009] Hastie, T., R. Tibshirani, J. H. Friedman, and J. H. Friedman (2009).  The elements of statistical learning: data mining, inference, and prediction, Volume 2.  Springer.
  * Ke et al. [2017] Ke, G., Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu (2017).  Lightgbm: A highly efficient gradient boosting decision tree.  Advances in neural information processing systems 30.
  * Klusowski [2021] Klusowski, J. M. (2021).  Universal consistency of decision trees in high dimensions.  arXiv preprint arXiv:2104.13881.
  * Klusowski and Siegel [2023] Klusowski, J. M. and J. W. Siegel (2023).  Sharp convergence rates for matching pursuit.  arXiv preprint arXiv:2307.07679.
  * Li et al. [2019] Li, X., Y. Wang, S. Basu, K. Kumbier, and B. Yu (2019).  A debiased mdi feature importance measure for random forests.  Advances in Neural Information Processing Systems 32.
  * Livshitz and Temlyakov [2003] Livshitz, E. and V. Temlyakov (2003).  Two lower estimates in greedy approximation.  Constructive approximation 19, 509–523.
  * Lockley and Eizaguirre [2021] Lockley, E. C. and C. Eizaguirre (2021).  Effects of global warming on species with temperature-dependent sex determination: Bridging the gap between empirical research and management.  Evolutionary Applications 14(10), 2361–2377.
  * Loecher [2020] Loecher, M. (2020).  Unbiased variable importance for random forests.  Communications in Statistics - Theory and Methods 51, 1413–1425.
  * Lou et al. [2007] Lou, X.-Y., G.-B. Chen, L. Yan, J. Z. Ma, J. Zhu, R. C. Elston, and M. D. Li (2007).  A generalized combinatorial approach for detecting gene-by-gene and gene-by-environment interactions with application to nicotine dependence.  The American Journal of Human Genetics 80(6), 1125–1137.
  * Louppe et al. [2013] Louppe, G., L. Wehenkel, A. Sutera, and P. Geurts (2013).  Understanding variable importances in forests of randomized trees.  In Advances in neural information processing systems, pp. 431–439.
  * Mallat and Zhang [1993] Mallat, S. G. and Z. Zhang (1993).  Matching pursuits with time-frequency dictionaries.  IEEE Transactions on signal processing 41(12), 3397–3415.
  * Saltelli et al. [2010] Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola (2010).  Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index.  Computer physics communications 181(2), 259–270.
  * Scornet [2023] Scornet, E. (2023).  Trees, forests, and impurity-based variable importance in regression.  In Annales de l’Institut Henri Poincare (B) Probabilites et statistiques, Volume 59, pp. 21–52. Institut Henri Poincaré.
  * Sobolprime [1993] Sobolprime, I. (1993).  Sensitivity analysis for nonlinear mathematical models.  Mathematical Modeling & Computational Experiment 1, 407–414.
  * Strobl et al. [2008] Strobl, C., A.-L. Boulesteix, T. Kneib, T. Augustin, and A. Zeileis (2008).  Conditional variable importance for random forests.  BMC bioinformatics 9(1), 1–11.
  * Strobl et al. [2007] Strobl, C., A.-L. Boulesteix, A. Zeileis, and T. Hothorn (2007).  Bias in random forest variable importance measures: Illustrations, sources and a solution.  BMC bioinformatics 8(1), 1–21.
  * Székely and Rizzo [2009] Székely, G. J. and M. L. Rizzo (2009).  Brownian distance covariance.  The annals of applied statistics, 1236–1265.
  * Tan et al. [2022] Tan, Y. S., A. Agarwal, and B. Yu (2022).  A cautionary tale on fitting decision trees to data from additive models: generalization lower bounds.  In International Conference on Artificial Intelligence and Statistics, pp. 9663–9685. PMLR.
  * Tan et al. [2022] Tan, Y. S., C. Singh, K. Nasseri, A. Agarwal, and B. Yu (2022).  Fast interpretable greedy-tree sums (figs).  arXiv preprint arXiv:2201.11931.
  * Zhou and Hooker [2021] Zhou, Z. and G. Hooker (2021).  Unbiased measurement of feature importance in tree-based methods.  ACM Transactions on Knowledge Discovery from Data (TKDD) 15, 1–21.

Supplementary Material to “Analyze Additive and Interaction Effects via
Collaborative Trees”

Chien-Ming Chi

The Supplementary Material provides supplementary content related to
hyperparameter tuning (Section A.1), algorithm implementation (Section A.2), a
runtime analysis (Section A.3), details of the Sobol’ indices (Section A.4), a
larger version of Figure 1 and detailed numerical results for Figure 1
(Section A.5.1), and simulation experiments (Section A.5.2 and Section A.6).
Additionally, Section B includes the proofs of Theorems 1–2, Corollary 1, and
several technical lemmas. The notation used throughout remains consistent with
that of the main paper. All Python codes used in the paper, including the
implementation of Collaborative Trees Ensemble, are available upon request.

##  Appendix A Supplementary Material

###  A.1 Hyperparameter tuning

Table 2 displays all hyperparameters of Collaborative Trees Ensemble
introduced in Section 2.3, where the number of aggregated predictors is
denoted by n_estimators. In Table 2, n_bins is required when transforming
certain continuous or categorical variables into feature groups. Additionally,
when $\textsf{random\\_update}=0$, Collaborative Trees randomly update one
node sampled from the update node list $W$ at each round.

We use the Python package hyperopt from https://hyperopt.github.io/hyperopt/
for tuning all predictive models. The tuning procedure is introduced as
follows. The full dataset is split into three subsamples, representing the
training, validation, and test sets, which constitute 48%, 32%, and 20% of the
full dataset, respectively. To determine the optimal hyperparameters for our
predictive models, we assess $B$ distinct sets. Utilizing the Python package
hyperopt, we perform hyperparameter sampling in each round. The training set
is used to train tree models based on the sampled hyperparameters.
Subsequently, the model’s performance is evaluated on the validation set to
score each set of hyperparameters. After conducting this process for $B$
rounds, we select the set with the highest evaluation scores as the best
hyperparameters. A final tree model is then trained using this optimal set,
taking into account both the training and validation samples, which represent
80% of the full dataset. Finally, we compute a R2 score on the test set for
assessing the model’s predictive performance. When prediction evaluation is
unnecessary, we split the sample into 60% for training and 40% for validation.

The hyperparameter spaces of Collaborative Trees, XGBoost, and Random Forests
are respectively given in Tables 2–4. In Table 2, n_estimators is the number
of boostrap aggregated trees, n_trees is the number of collaborative trees
$K$, alpha is the probability weights ciefficients used in the softmax
function, n_bins is the number of bins for binning continuous and categorical
features, and that random_update, min_samples_split, min_samples_leaf,
max_depth have been introduced in Section 2.3. The choice of the n_trees
parameter depends on the number of significant components in the data
generating models. Based on our experience in the data study, we recommend an
upper limit of $12$ for n_trees. Meanwhile, Tables 3–4 respectively present
the hyperparameters of XGBoost and Random Forests. Additional details about
these hyperparameters can be found on the respective package websites.

Parameter name | Parameter type | Parameter space
---|---|---
n_estimators (# baggings) | integer | 100
n_trees ($K$ in Algorithm 1) | integer | $\\{6,\cdots,12\\}$
alpha | non-negative number | $\\{0.001,1,10,100,10000,\infty\\}$
min_samples_split | integer | $\\{5,10,15,20,30\\}$
min_samples_leaf | integer | $\\{0,5,10,15,20,30\\}$
n_bins | integer | $\\{5,7,10,15,20,40\\}$
random_update | non-negative number | $\\{0,0.0001,0.001,0.01,0.1,1\\}$
max_depth | integer | $\\{3,5,10,20,30,\infty\\}$
Table 2: Hyperparameter space for tuning Collaborative Trees Ensemble. Note that even when random_update is less than 1, the first $2K$ updates of Collaborative Trees still adhere to our tree model regularization (2) Parameter name | Parameter type | Parameter space
---|---|---
n_estimators | integer | 1000
gamma | non-negative real number |  $\exp{(Z)}$, $Z\sim\textnormal{Uniform}[-8\log_{e}(10),\log_{e}(7)]$
reg_alpha | non-negative real number |  $\exp{(Z)}$, $Z\sim\textnormal{Uniform}[-8\log_{e}(10),\log_{e}(100)]$
reg_lambda | non-negative real number |  $\exp{(Z)}$, $Z\sim\textnormal{Uniform}[\log_{e}(0.8),\log_{e}(4)]$
learning_rate | non-negative real number |  $\exp{(Z)}$, $Z\sim\textnormal{Uniform}[-5\log_{e}(10),\log_{e}(0.7)]$
subsample | positive real number | $\textnormal{Uniform}[0.5,1]$
colsample_bytree | positive real number | $\textnormal{Uniform}[0.5,1]$
colsample_bylevel | positive real number | $\textnormal{Uniform}[0.5,1]$
min_child_weight | non-negative integer | Uniform$\\{0,\cdots,20\\}$
max_depth | non-negative integer | Uniform$\\{2,\cdots,15\\}$
Table 3: Hyperparameter space for tuning XGBoost. Details of these parameters can be found on the website https://xgboost.readthedocs.io/en/stable/parameter.html.  Parameter name | Parameter type | Parameter space
---|---|---
n_estimators | integer | 1000
gamma | non-negative real number | $\textnormal{Uniform}[0,1]$
min_samples_split | integer | $\\{1,\cdots,20\\}$
min_samples_leaf | integer | $\\{2,\cdots,20\\}$
min_impurity_decrease | non-negative real number | $\textnormal{Uniform}\\{0,0.01,0.02,0.05\\}$
max_depth | non-negative integer | Uniform$\\{5,10,20,\infty\\}$
criterion | string | $[\textnormal{squared\\_error},\textnormal{absolute\\_error}]$
Table 4: Hyperparameter space for tuning Random Forests. Details of these
parameters can be found on the website https://scikit-
learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html

###  A.2 Collaborative Trees Ensemble implementation

Collaborative Trees Ensemble includes the following implementations in
addition to Algorithm 1, with all hyperparameters given in Table 2. 1.) To
enhance prediction, we build Collaborative Trees with bootstrap aggregating,
or bagging. 2.) A node of depth less than max_depth is valid for splitting if
its subsample size is larger than
$\max\\{\textsf{min\\_samples\\_split},\textsf{min\\_samples\\_leaf}\\}$. If a
node has at least two child nodes whose subsample sizes are larger than
min_samples_leaf, then these child nodes are valid for tree update. 3.) To
reduce computation time, after the initial $2K$ updates, we update nodes in
$W^{\dagger}$, randomly sampled from $W$ with
$\texttt{\\#}W^{\dagger}=\max\\{\textsf{random\\_update}\times\texttt{\\#}W,1\\}$.
See Section 2.2 for details of the update node list $W$. 4.) To reduce
estimation variance, Collaborative Trees make update decisions based on
probability weights derived from split scores. At each round of update with
some update node list $W$, we compute a vector of split scores
$\big{(}\sum_{(C,k)\in Q}\textnormal{SplitScore}(C,m)\big{)}-\lambda(W,Q)$ as
in (2) for each $Q\in W$, denoted by $\vv{SS}$. We update the split sampled
from the probability vector
$\textnormal{softmax}(\textsf{alpha}\times\vv{SS})$, where
$\textnormal{softmax}(\vv{z})=\big{[}\sum_{j=1}^{k}\exp{(z_{j})}\big{]}^{-1}\times(\exp{(z_{1})},\cdots,\exp{(z_{k})})$
for $\vv{z}\in\mathbb{R}^{k}$. The case with $\textsf{alpha}=\infty$
corresponds to maximization as in (2). This split sampling approach draws
inspiration from BART [12], known to enhance tree prediction accuracy in
certain applications. Here, to streamline the sampling process, we adopt a
softmax transform on the split scores.

###  A.3 Computational time analysis

We extracted a subsample of size $n=10,000$ from the “superconduct” dataset in
[20], sourced from OpenML (https://www.openml.org/). This dataset comprises 79
features. This sample is also used in Table 9. The runtime evaluation for
Collaborative Trees Ensemble with a single bagging tree is conducted with the
following hyperparameters: $K\in\\{1,3,5,7,9,11\\}$, min_samples_split = 5,
min_samples_leaf = 5, $\textsf{random\\_update}\in\\{0,0.001,0.01,0.1,1\\}$,
$\textsf{alpha}=\infty$, $\textsf{max\\_depth}=20$, $\textsf{n\\_bins}=None$.
The specifics of the hyperparameters for Collaborative Trees Ensemble are
outlined in Section 2.3 and Section A.2.

The evaluation was performed on a Mac laptop with a 2.6 GHz 6-Core Intel Core
i7 processor and macOS Sonoma Version 14.1.2. The environment details are as
follows: Spyder version 5.4.3 (conda), Python version 3.10.12 64-bit, Qt
version 5.15.2, PyQt5 version 5.15.7, Operating System: Darwin 23.1.0.

Although the parameters related to the tree structure (min_samples_split,
min_samples_leaf, max_depth) can significantly impact the model’s training
runtime, the evaluation focuses on the number of collaborative trees $K$ and
random_update, which determine the number of associated nodes considered at
each update round. The results of runtime are reported in Table 5.

From Table 5, it takes approximately 5.5 hours
($200*\frac{1}{60}*\frac{1}{60}*100\approx 5.5$) to train a Collaborative
Trees Ensemble model with full updates (i.e., random_update = 1), $K=11$, and
$100$ bagging trees on a sample of size $10,000$ with $79$ features. It is
important to note that while the current implementation of Collaborative Trees
Ensemble in Python is not yet optimized for computation (codes are written
using the Python numpy package), we have implemented parallel training to
expedite processing.

Additionally, it is worth mentioning that the number of bagging trees could
potentially be much smaller than our default setting of 100. Furthermore,
setting $\textsf{random\\_update}=1$ may not be necessary across the entire
hyperparameter space in Table 2, which could significantly reduce
computational time, as indicated in Table 5.

$K$ $\Big{\backslash}$ | random_update | 0 | 0.01 | 0.1 | 1
---|---|---|---|---|---
1 |  | 2 | 2 | 4 | 14
3 |  | 5 | 5 | 11 | 43
5 |  | 8 | 9 | 19 | 93
7 |  | 13 | 12 | 30 | 115
9 |  | 15 | 16 | 35 | 142
11 |  | 21 | 22 | 51 | 203
Table 5: Table showing the computational time evaluation for Collaborative
Trees Ensemble with a single bagging tree (i.e., the number of bagging trees =
1). The table presents the computational time for various hyperparameter
configurations, including the number of collaborative trees ($K$) and
random_update settings. Time measurements are displayed in seconds, rounded to
the nearest second. Each entry represents the average runtime over 5
repetitions.

###  A.4 Sobol’ indices in Section 5.1

We implemented the Sobol’ indices using the R package sensitivity. This
package requires a trained predictive model for calculating various importance
measures. Similar to our approach for MDI and CPI, we trained a Random Forests
model for the Sobol’ indices. It is worth noting that there are several
implementations available for calculating the Sobol’ indices with a predictive
model. We opted for the ‘sobolSalt‘ function in the package due to its user-
friendly interface for calculating up to second-order importance measures.

For the Sobol’ indices, we estimated the total Sobol’ indices [35, 5], first-
order indices, which are similar to XMDIjj’s, and the second-order indices for
each pair of features.

It is important to mention that the Random Forests algorithm in R differs
slightly from its Python counterpart. We followed a similar hyperparameter
tuning procedure as in Python as in Section A.1, but for simplicity, we do not
provide detailed reporting here.

###  A.5 Simulation experiment details for Section 6.1

####  A.5.1 Numerical analysis results for Section 6.1

Figure 6 provides an enlarged version of the image displayed in Figure 1.
Additionally, Table 6 presents detailed numerical values used in Figure 1,
rounded to three decimal places. In Table 6, standardized feature importance
measures are computed as $\frac{\textnormal{XMDI}_{m}}{0.0167}$, where
$0.0167$ represents the sample variance of the response variable (sex ratio of
female to total births). Standardized additive effects are calculated as
$\frac{\textnormal{XMDI}_{mm}}{\textnormal{XMDI}_{m}}$ for each of the $7$
feature groups. For the interaction matrix, values are obtained as
$\frac{\textnormal{XMDI}_{ij}}{\textnormal{XMDI}_{j}}$, where $j$ is the
column index of the interaction effect matrix. The names of the feature groups
are indicated in the table.

Let us illustrate how to use the table with an example of computing
$\textnormal{XMDI}_{73}$, which represents the interaction effect between area
and temperature. To compute $\textnormal{XMDI}_{73}$ from Table 6, we follow
these steps: 1.) Obtain the value of
$\frac{\textnormal{XMDI}_{73}}{\textnormal{XMDI}_{3}}$ from the (7, 3) entry
of the interaction matrix. 2.) Get $\frac{\textnormal{XMDI}_{3}}{0.0167}$ from
the row of the standardized feature importance. 3.) Multiply the results of
steps 1 and 2 by 0.0167 to obtain the final value of $\textnormal{XMDI}_{73}$.

As a result, the interaction effect between area and temperature is estimated
to be approximately $0.331\times 0.009\times 0.167\approx 0.005$. Similarly,
the interaction effect between species and temperature is approximately 0.032.
All other interaction effects, when rounded to three decimal places, are
negligible and thus are not reported here. Recall that the thickness of each
edge depends on the size of the corresponding interaction effect, while its
grayness depends on the corresponding standardized interaction effect. See
Section 3.1 for details of network diagrams.

| X1 | X2 | X3 | X4 | X5 | X6 | X7
---|---|---|---|---|---|---|---
Standardized Feature Importance | 0.353 | 0.005 | 0.090 | 0.003 | 0.002 | 0.034 | 0.723
Standardized Additive Effects | 0.462 | 0.580 | 0.637 | 0.461 | 0.815 | 0.882 | 0.695
Standardized Interaction Effects
From $\Big{\backslash}$ To | X1 | X2 | X3 | X4 | X5 | X6 | X7
X1 : species | NA | 0.350 | 0.014 | 0.038 | 0.038 | 0.032 | 0.259
X2 : subspecies | 0.001 | NA | 0.003 | 0.029 | 0.006 | 0.005 | 0.002
X3 : area | 0.004 | 0.056 | NA | 0.056 | 0.048 | 0.032 | 0.041
X4 : RMU | 0.000 | 0.016 | 0.002 | NA | 0.026 | 0.017 | 0.001
X5 : amplitude | 0.000 | 0.002 | 0.001 | 0.017 | NA | 0.001 | 0.000
X6 : incubation periods (days) | 0.003 | 0.030 | 0.012 | 0.189 | 0.020 | NA | 0.002
X7 : temperature | 0.531 | 0.282 | 0.331 | 0.210 | 0.047 | 0.032 | NA
Table 6:  Detailed numerical values used in Figure 1, rounded to three decimal
places. Standardized feature importance measures are computed as
$\frac{\textnormal{XMDI}_{m}}{0.0167}$, where $0.0167$ represents the sample
variance of the response variable (sex ratio). Standardized additive effects
are calculated as $\frac{\textnormal{XMDI}_{mm}}{\textnormal{XMDI}_{m}}$ for
each of the $7$ feature groups. For the interaction matrix, values are
obtained as $\frac{\textnormal{XMDI}_{ij}}{\textnormal{XMDI}_{j}}$, where $j$
is the column index of the interaction effect matrix. The term NA denotes “not
applicable.” Figure 6: Detailed view of Figure 1.

####  A.5.2 Additional simulation experiments for Section 6.1

Figure 7: Network diagrams for data simulated from Model (A.1) with $p=100$
and $\lambda=0.1$. The picture here concerns only the first 10 features.
$\lambda=0.1$
---
(I) | (II)
0.95 | 0.97
Table 7: XMDI of Collaborative Trees Ensemble. (I) The rates of successfully
separating significant and insignificant features in Model (A.1), based on the
overall feature importance measures. (II) The rates of
$\min_{j\in\\{2,6,9\\}}\textnormal{XMDI}_{10,j}$ being the maximum interaction
effect among all other pairs. All rates are derived from $100$ simulation
experiments with $p=100$ features.

In this section, we present additional experiments similar to those in Section
5.1 showing that the XMDI can practically distinguish interaction effects
stemming from the same feature group. The response model of simulation
experiment is given as follows.

$\begin{split}Y&=2X_{10}+10\sin{(\pi(X_{2}-0.5)(X_{10}-0.5))}+10\sin{(\pi(X_{6}-0.5)(X_{10}-0.5))}\\\
&\qquad+10\sin{(\pi(X_{9}-0.5)(X_{10}-0.5))}+\varepsilon.\end{split}$ (A.1)

The sample size is $1000$ for each experiment, with $\lambda=0.1$ and $p=100$
for Figure 7 and Table 7; Figure 7 only displays the results of
$X_{1},\cdots,X_{10}$. Here, $p$ denotes the number of features. Other model
assumptions are the same as those for Model (9).

We use $20$ runs of optimizations for hyperparameter tuning (see Section A.1)
and report important hyperparameters: $\textsf{n\\_trees}=K=11$ and
$\textsf{n\\_bins}=5$. This challenging learning example involves three XOR
components. Hence, we consider a sample of size $1000$, which is larger than
the ones of $500$ samples in Section 5.1. In Figure 7, Collaborative Trees
Ensemble identifies four relevant features ($X_{2}$, $X_{6}$, $X_{9}$,
$X_{10}$) based on their XMDIj importance measure. Furthermore, all
interaction effects are visually identified in Figure 7. On the other hand,
the column (II) of Table 7 shows that Collaborative Trees Ensemble
consistently recognizes these three interaction effects as the strongest among
all pairs over 100 repeated experiments. Additionally, column (I) of Table 7
demonstrates that XMDIj effectively differentiates between significant and
insignificant features across all experiments. Here, we emphasize again that
our goal is to verify empirically whether there exists a threshold for
separating significant and insignificant features in each experiment. However,
we do not aim at finding such a threshold.

The results presented in Table 7 complement those of Table 1, demonstrating
the effectiveness of XMDI in distinguishing interaction effects originating
from the same feature. Therefore, the practical performance of XMDI is not
contingent on our technical condition requiring tree heredity in Condition 3
(but the addition of the weak additive component $2X_{10}$ in Model (A.1) does
stabilize the numerical results). These findings support the application of
Collaborative Tree Ensemble in Section 6.1 for the embryo growth dataset
study.

###  A.6 Additional details of prediction evaluation for Section 6.2

| Win Rates
---|---
CTE | 0.7528 | 0.8069 | 0.8909
XGB | 0.4662 | 0.5986 | 0.7011
RF | 0.1612 | 0.2439 | 0.3413
Table 8: For each experiment, an average win rate over $19$ datasets is
calculated. From left to right columns, the lowest, the average, and the
highest win rates across $10$ experiments are respectively reported for each
method.

Following the approach outlined in [20], we utilize 19 selective datasets from
OpenML (https://www.openml.org/), each containing no more than 10,000 samples,
for prediction evaluation based on the goodness-of-fit measure R2 value. The
calculation of each R2 value involves splitting the full sample into training,
validation, and test sets, constituting 48%, 32%, and 20% of the full dataset,
respectively. This division is employed for hyperparameter tuning based on the
training and validation samples, as well as prediction loss calculation based
on the test sample. A detailed description of the hyperparameter tuning
procedure, involving 100 optimization runs for this evaluation task, is
provided in Section A.1.

To summarize the comparison, we calculate adjusted win rates for the resulting
R2 values. For instance, if three methods (CTE, XGBoost, Random Forests) yield
R2 values of $0.13$, $0.135$, and $0.142$, they receive scores of $0$,
$\frac{0.135-0.13}{0.142-0.13}$, and $1$, respectively, rounded to the fourth
decimal if necessary. This performance measure, also adopted in [20] (but
their calculation involves more than three models), ensures a fair comparison
among the methods across datasets. The term “adjusted win rate” is used to
account for the interpolation when more than two models are compared. When
only two models are compared, win rates (1 if a model wins, 0 otherwise) are
calculated based on their resulting R2 scores. For each dataset and method, we
conduct 10 independent repeated experiments. We calculate average win rates
denoted as $wr_{1}^{(k)},\cdots,wr_{10}^{(k)}$ for $10$ experiments and the
$k$th dataset, with $k\in\\{1,\cdots,19\\}$. In Table 8, the reported win
rates for each method represent $\min_{1\leq b\leq 10}wr_{b}$,
$\frac{1}{10}\times\sum_{1\leq b\leq 10}wr_{b}$, and $\max_{1\leq b\leq
10}wr_{b}$ from left to right, where
$wr_{b}=\frac{1}{19}\sum_{k=1}^{19}wr_{b}^{(k)}$. It is important to note that
the entire evaluation procedure, including three sample splits, is repeated
independently 10 times for each dataset.

From Table 8, Collaborative Trees Ensemble consistently outperforms XBGoost
and Random Forests. These benchmark methods, recognized for their strong
performance against modern deep learning models in regression and
classification prediction tasks [20], are surpassed by our tree model. The
results underscore our model’s superior predictive capabilities for
applications with medium sample sizes. Further experiments targeting
regression applications with large sample sizes ($>10,000$) and classification
tasks are reserved for future work.

In Table 9, we present the minimum, average, and maximum R2 scores for four
models across 10 experiments for each dataset, arranged from left to right.
The models included are Collaborative Trees Ensemble, XGBoost, Random Forests,
and least squares regression (LS), with the LS serving as a linear benchmark
model. It is worth noting that the LS is not included in Table 8 to prevent
potential overestimation of the other models due to interpolation. These
results are comparable to Figure 13 of [20], where a snippet of the results
from the sample “house” is displayed in Figure 8.

In Figure 8, they report all results from the first run to the 100th run of
optimization, but they do not give the detailed R2 scores as in our Table 9
since they report many other details. In Figure 8, the green line represents
the R2 scores of XGBoost, while the brown line represents the results of
Random Forests. We refer to [20] for other details. Our results in Table 9
align closely with theirs up to the 100th column. For instance, the average R2
score results for XGBoost are approximately 0.85 for the ”house” sample in
both Table 9 and Figure 8. However, there is an exception with the sample
”delay_zurich_transport,” where XGBoost exhibits higher variation and a lower
average R2 score in Table 9.

All Python source codes used to replicate these results are available upon
request.

Figure 8: A clip from the “house” sample in Figure 13 of [20]. The x-axis is the number of optimizations in hyperparameter tuning, ranging from $1$ to $100$ (each column dot is 10 runs), while the y-axis is the R2 scores.  | cpu_act | pol | elevators
---|---|---|---
CTE | (0.9848, 0.9872, 0.9894) | (0.9840, 0.9879, 0.9918) | (0.8884, 0.8972, 0.9057)
XGB | (0.9851, 0.9873, 0.9889) | (0.9861, 0.9897, 0.9943) | (0.8837, 0.8968, 0.9078)
RF | (0.9793, 0.9828, 0.9849) | (0.9797, 0.9850, 0.9887) | (0.7757, 0.8016, 0.8282)
LS | (0.7029, 0.7170, 0.7523) | (0.4471, 0.4625, 0.4790) | (0.7976, 0.8163, 0.8361)
| Ailerons | houses | house_16H
CTE | (0.8369, 0.8592, 0.8719) | (0.8204, 0.8454, 0.8673) | (0.3910, 0.5332, 0.6559)
XGB | (0.8171, 0.8354, 0.8475) | (0.8040. 0.8352, 0.8566) | (0.3960, 0.5278, 0.6574)
RF | (0.8147, 0.8395, 0.8503) | (0.7833, 0.8117, 0.8279) | (0.4740, 0.5386, 0.6460)
LS | (0.7517, 0.8132, 0.8280) | (0.6070, 0.6357, 0.6747) | (0.1624, 0.2422, 0.3259)
| Brazilian_houses | Bike_Sharing_Demand | nyc-taxi-green-dec-2016
CTE | (0.9776, 0.9952, 0.9998) | (0.6797, 0.6962, 0.7244) | (0.5184, 0.5646, 0.5884)
XGB | (0.9791, 0.9954, 0.9996) | (0.6511, 0.6894, 0.7217) | (0.4954, 0.5292, 0.5620)
RF | (0.9578, 0.9859, 0.9980) | (0.6625, 0.6863, 0.7043) | (0.4631, 0.5371, 0.5702)
LS | (0.4390, 0.7615, 0.8399) | (0.3073, 0.3288, 0.3490) | (0.2636, 0.3067, 0.3522)
| sulfur | medical_charges | MiamiHousing2016
CTE | (0.7812, 0.8480, 0.9049) | (0.9755, 0.9792, 0.9846) | (0.9272, 0.9348, 0.9417)
XGB | (0.8264, 0.8678, 0.8917) | (0.9749, 0.9782, 0.9834) | (0.9206, 0.9305, 0.9350)
RF | (0.7556, 0.8107, 0.8720) | (0.9748, 0.9784, 0.9843) | (0.8947, 0.9113, 0.9249)
LS | (0.3379, 0.3874, 0.4364) | (0.7936, 0.8188, 0.8438) | (0.6803, 0.7134, 0.7320)
| yprop_4_1 | abalone | delay_zurich_transport
CTE | (0.0539, 0.0730, 0.0910) | (0.5136, 0.5427, 0.5782) | (0.0113, 0.0236, 0.0393)
XGB | (0.0182, 0.0599, 0.0955) | (0.4329, 0.5182, 0.5532) | (-0.2871, -0.0324, 0.0178)
RF | (0.0537, 0.0827, 0.1020) | (0.4844, 0.5336, 0.5816) | (0.0104, 0.0213, 0.0338)
LS | (0.0278, 0.0519, 0.0690) | (0.4628, 0.4946, 0.5841) | (0.0013, 0.0059, 0.0115)
| wine_quality | diamonds | house_sales
CTE | (0.4989, 0.5328, 0.5607) | (0.9386, 0.9446, 0.9482) | (0.8803, 0.8920, 0.9064)
XGB | (0.4954, 0.5389, 0.5759) | (0.9391, 0.9428, 0.9466) | (0.8788, 0.8880, 0.9017)
RF | (0.4014, 0.4885, 0.5390) | (0.9391, 0.9441, 0.9467) | (0.8603, 0.8696, 0.8906)
LS | (0.2457, 0.2878, 0.3339) | (0.7999, 0.9122, 0.9384) | (0.7321, 0.7472, 0.7662)
| superconduct |  |
CTE | (0.8963, 0.9065, 0.9152) |  |
XGB | (0.8989, 0.9099, 0.9218) |  |
RF | (0.8935, 0.9048, 0.9158) |  |
LS | (0.7155, 0.7342, 0.7523) |  |
Table 9:  The table displays the minimum, average, and maximum R2 scores
across 10 experiments for each of the 19 selective datasets, arranged from
left to right. Note that the results of the least square method (LS) are
included in this table, but they are excluded from the adjusted win rate
calculation.

##  Appendix B Technical proofs for main paper

In this section, we present the technical proofs of Theorem 1 in Section 4.1,
Corollary 1 in Section 4.1, and Theorem 2 in Section 4.2. These proofs are
given respectively in Sections B.1–B.3.

###  B.1 Proof of Theorem 1

Proof of Theorem 1:  Let us introduce notation for our proofs. First, recall
that the output tree predictor in Algorithm 1 is denoted by
$\sum_{k=1}^{K}T_{k}(\vv{x})$. We will introduce another notation for sample
tree predictors as used in Section 3 and (A.4) below, with subscripts
indicating the update round index. This notation is crucial for presenting our
technical analysis, for which we need to introduce the following set of
notation.

For each $J\subset\\{1,\cdots,M\\}$, let $C_{l,J}$ for $l\in\mathcal{L}(J)$
denote a partition of the feature space (not necessarily a partition of
$\mathbb{R}^{p}$) based on features in $\mathcal{X}(J)$ such that
$\sum_{l\in\mathcal{L}(J)}\boldsymbol{1}_{\boldsymbol{X}\in C_{l,J}}=1$ almost
surely. Specifically, it is required that

$\\{C_{l,J}:l\in\mathcal{L}(J)\\}=\\{\cap_{m\in
J}A_{m}:A_{m}\in\mathcal{C}_{m}\textnormal{ for each }m\in J\\},$ (A.2)

where

$\mathcal{C}_{m}=\begin{cases}\\{\\{\vv{x}\in\mathbb{R}^{p}:x_{j}>0\\},\\{\vv{x}\in\mathbb{R}^{p}:x_{j}\leq
0\\}:j\in\mathcal{X}(m)\\}\hskip 17.07182pt\textnormal{ if
}\texttt{\\#}\mathcal{X}(m)=1,\\\
\\{\\{\vv{x}\in\mathbb{R}^{p}:x_{j}>0\\}:j\in\mathcal{X}(m)\\}\hskip
113.81102pt\textnormal{ if }\texttt{\\#}\mathcal{X}(m)>1.\end{cases}$

As an illustrative example, let $J=\\{l,k\\}\subset\\{1,\cdots,M\\}$,
$\mathcal{X}(l)=\\{j_{1}\\}\subset\\{1,\cdots,p\\}$, and
$\mathcal{X}(k)=\\{j_{2},j_{3}\\}\subset\\{1,\cdots,p\\}$. Then, we define
$\mathcal{L}(J)=\\{1,2,3,4\\}$, and define $C_{l,J}$’s such that
$\\{C_{l,J}:l\in\mathcal{L}(J)\\}=\big{\\{}\\{\vv{x}\in\mathbb{R}^{p}:x_{j_{1}}>0,x_{j_{2}}>0\\},\\{\vv{x}\in\mathbb{R}^{p}:x_{j_{1}}>0,x_{j_{3}}>0\\},\\{\vv{x}\in\mathbb{R}^{p}:x_{j_{1}}\leq
0,x_{j_{2}}>0\\},\\{\vv{x}\in\mathbb{R}^{p}:x_{j_{1}}\leq
0,x_{j_{3}}>0\\}\big{\\}}$. As a result,
$\sum_{l\in\mathcal{L}(J)}\boldsymbol{1}_{\boldsymbol{X}\in C_{l,J}}=1$ almost
surely due to Condition 1. Here, it is worth noting that the precise
definition of $\mathcal{L}(J)$ and $C_{l,J}$ is not crucial as long as (A.2)
is satisfied.

Let $\\{\widehat{J}_{l}\\}_{l=1}^{2K}$ be the sample optimal path of sets
given in Section 4.1. For each set $J$ and $0<s\leq 2K$, we define functions
$R_{s,J}^{\star}:\mathbb{R}^{p}\mapsto\mathbb{R}$,
$\widehat{U}_{s,J}:\mathbb{R}^{p}\mapsto\mathbb{R}$, and
$\widehat{U}_{s}:\mathbb{R}^{p}\mapsto\mathbb{R}$ such that

$\begin{split}\widehat{U}_{s,J}(\vv{x})&=\widehat{U}_{s-1}(\vv{x})+\sum_{l\in\mathcal{L}(J)}\boldsymbol{1}_{\vv{x}\in
C_{l,J}}\frac{\sum_{i\in\\{\boldsymbol{X}_{i}\in
C_{l,J}\\}}\left[Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i})\right]}{\texttt{\\#}\\{i:\boldsymbol{X}_{i}\in
C_{l,J}\\}\vee 1},\\\
\widehat{U}_{s}(\vv{x})&=\widehat{U}_{s,\widehat{J}_{s}}(\vv{x}),\\\
R_{s,J}^{\star}(\boldsymbol{X})&=R_{s-1}^{\star}(\boldsymbol{X})+\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}),\end{split}$
(A.3)

where $R_{l}^{\star}$’s are given as in (7), and $\widehat{U}_{s}(\vv{x})$ is
recursively defined by (A.3) given $\widehat{U}_{0}(\vv{x})=0$. In addition,
we set $\widehat{U}_{0,J}(\vv{x})=R_{0,J}^{\star}(\vv{x})=0$ for each
$\vv{x}\in\mathbb{R}^{p}$ and every set $J$. It is important to note that
$R_{s,\widehat{J}_{s}}^{\star}(\boldsymbol{X})=R_{s}^{\star}(\boldsymbol{X})$.

To clarify, $\widehat{U}_{s}$ in (A.3) represents the sample tree predictor at
the (end of the) $s$th round of updates, obtained from Algorithm 1 with $s\leq
2K$. To see this, let us consider the update at the $s$th round with $s\leq
K$, and that the feature group to be split has a single binary feature with
$\widehat{J}_{s}=\\{q\\}$. Then,
$\\{C_{1,\widehat{J}{s}},C_{2,\widehat{J}{s}}\\}=\\{\\{\vv{x}\in\mathbb{R}^{p}:x_{q}\leq
0\\},\\{\vv{x}\in\mathbb{R}^{p}:x_{q}>0\\}\\}$ represent the two child nodes
of some root node at Step 12 of Algorithm 1. The tree update at Step 12 of
Algorithm 1 is equivalent to the recursive definition in (A.3) when
considering the residuals
$\\{Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i})\\}_{i=1}^{n}$ at the beginning
of the $s$th round. The distinction lies in emphasizing the round index,
denoted by subscript $s$ in $\widehat{U}_{s}(\vv{x})$, rather than the tree
index, denoted by subscript $k$ in $T_{k}(\vv{x})$. Similar arguments can be
applied for general cases with $K<s\leq 2K$, but we omit the details for
simplicity. Note that this notation for predictors has already been introduced
in Section 3. Furthermore, while the expressions of sample tree predictors
presented here may not align with the sample tree predictor in Algorithm 1
when $s>2K$, this discrepancy does not affect our arguments. Our focus is
specifically on the first $2K$ updates involving $K$ collaborative decision
trees in Theorem 1.

Considering the determination of the next feature group to be split in (2),
and the definition of $\widehat{J}_{l}$’s in Section 4.1, it holds that for
$0<s\leq 2K$,

$\widehat{J}_{s}=\operatorname*{arg\,max}_{J\in\Theta_{s}}\left\\{\left[\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right]-\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\right\\},$
(A.4)

where $\Theta_{1}=\cdots=\Theta_{K}=\\{\\{1\\},\cdots,\\{M\\}\\}$, and that

$\begin{split}A_{s}&=\\{\\{l,k\\}:1\leq k\leq
M,k\not=l\in\widehat{J}_{s}\\}\textnormal{ for }0<s\leq K,\\\
\Theta_{s}&=\cup_{l\in Z_{s}}A_{l}\textnormal{ for }K<s\leq 2K,\\\
Z_{s+1}&=Z_{s}\backslash\\{\widehat{k}_{s}\\},\textnormal{ for
}K<s<2K,\end{split}$

where $1\leq\widehat{k}_{s}\leq K$ is such that $\widehat{J}_{s}\in
A_{\widehat{k}_{s}}$ (after $\widehat{J}_{s}$ has been decided) for $K<s\leq
2K$, and $Z_{K+1}=\\{1,\cdots,K\\}$. Ties are broken randomly. The expression
within the curly brackets in (A.4) represents a split score multiplied by $n$.
The specific set restriction for (A.4) is imposed by natural tree heredity
before the initial $2K$ updates, our update priority outlined in (2), and our
assumption that $\widehat{J}_{l}\neq\widehat{J}_{k}$ for each
$\\{l,k\\}\subset\\{1,\cdots,2K\\}$. Recall that $\\{j,j\\}=\\{j\\}$ according
to our definition of sets. The proof of (A.4) is postponed to Section B.1.1
after the end of the proof of Theorem 1.

By that
$\widehat{U}_{s,\widehat{J}_{s}}(\boldsymbol{X})=\widehat{U}_{s}(\boldsymbol{X})$
and
$R_{s,\widehat{J}_{s}}^{\star}(\boldsymbol{X})=R_{s}^{\star}(\boldsymbol{X})$
for each $s\geq 0$, and Lemma 3, it holds that for each $0<s\leq 2K$ and every
set $J$,

$\begin{split}&n^{-1}\left[\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right]-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&=\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}))^{2}]\\\
&\qquad+D_{2,s-1}(\widehat{J}_{s-1})-D_{2,s}(J)+D_{1,s-1}(\widehat{J}_{s-1})-D_{1,s}(J),\end{split}$
(A.5)

where

$\begin{split}D_{1,s}(J)&=\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\right)-\left(n^{-1}\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}\varepsilon_{i}^{2},\\\
D_{2,s}(J)&=\left(n^{-1}\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}\right)-\mathbb{E}(f(\boldsymbol{X})-R_{s,J}^{\star}(\boldsymbol{X}))^{2}.\end{split}$

We will establish upper bounds for $|D_{1,s}(J)|$’s and $|D_{2,s}(J)|$’s for
each $s\geq 0$ and every $J$ with $\texttt{\\#}J\leq 2$, starting with the
expressions of $R_{s,J}^{\star}(\boldsymbol{X})$ and
$R_{s}^{\star}(\boldsymbol{X})$ provided in (A.6) below. Additionally, we will
utilize several other preliminary results, including (A.6)–(A.10) below, to
analyze the upper bounds on $|D_{1,s}(J)|$’s and $|D_{2,s}(J)|$’s.

By (A.3) and (7), it holds that for every set $J$ and each $s>0$,

$\begin{split}R_{s,J}^{\star}(\boldsymbol{X})&=\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)})\\\
&\qquad+\sum_{l=1}^{s-1}\mathbb{E}(f(\boldsymbol{X})-R_{l-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(\widehat{J}_{l})})\\\
&=\left(\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}a_{s,\vv{c},J}\times\boldsymbol{1}_{\boldsymbol{X}_{\mathcal{X}(J)}=\vv{c}}\right)+\sum_{l=1}^{s-1}\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(\widehat{J}_{l})}}a_{l,\vv{c},\widehat{J}_{l}}\times\boldsymbol{1}_{\boldsymbol{X}_{\mathcal{X}(\widehat{J}_{l})}=\vv{c}},\\\
R_{s}^{\star}(\boldsymbol{X})&=\sum_{l=1}^{s}\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(\widehat{J}_{l})}}a_{l,\vv{c},\widehat{J}_{l}}\times\boldsymbol{1}_{\boldsymbol{X}_{\mathcal{X}(\widehat{J}_{l})}=\vv{c}},\end{split}$
(A.6)

where
$a_{l,\vv{c},J}=\mathbb{E}\\{f(\boldsymbol{X})-R_{l-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}=\vv{c}\\}$
for every set $J$, each $l>0$, and each
$\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}$. For the completeness of
notation, we define $a_{0,\vv{c},J}=0$ for all vector $\vv{c}$ and set $J$.
Recall that we have defined
$\mathbb{E}(h(\boldsymbol{X})|\boldsymbol{X}=\vv{c})=0$ when
$\mathbb{P}(\boldsymbol{X}=\vv{c})=0$ for every $\vv{c}\in\\{0,1\\}^{p}$ and
every measurable function $h:\mathbb{R}^{p}\mapsto\mathbb{R}$. Also, a
summation over an empty set is defined to be zero. Furthermore,
$\sum_{l=i}^{i-k}1$ is defined to be zero for every integer $i$ and every
positive integer $k$. By Lemma 2, (A.6), and the definition of $\iota_{s}$ in
(8), for every set $J$ with $\texttt{\\#}J\leq 2$ and each $s\geq 0$,

$\max_{\vv{c}\in\\{0,1\\}^{p}}|R_{s,J}^{\star}(\vv{c})|\leq
3\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}\times p_{\min}^{-1}+\iota_{s-1},$
(A.7)

where we note that $\max_{\vv{c}\in\\{0,1\\}^{p}}|R_{0,J}^{\star}(\vv{c})|=0$
by the definition in (A.3). In addition, since each feature group contains
either a single feature or a group of one-hot indicators, it holds that for
each $l>0$ and each subset $J\subset\\{1,\cdots,M\\}$ with $\texttt{\\#}J\leq
2$,

$\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}\boldsymbol{1}\\{|a_{l,\vv{c},J}|>0\\}\leq
M_{\mathcal{X}}^{2},$ (A.8)

in which $M_{\mathcal{X}}$ is defined in Condition 1.

By the definition of $\iota_{l}$’s in (8), Lemma 2, and (A.8), it holds that
for every $J$ with $\texttt{\\#}J\leq 2$ and each $s\geq 0$,

$\begin{split}\sum_{l=1}^{s-1}\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(\widehat{J}_{l})}}|a_{l,\vv{c},\widehat{J}_{l}}|&\leq
M_{\mathcal{X}}^{2}\iota_{s-1},\\\
\left(\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}|a_{s,\vv{c},J}|\right)+\sum_{l=1}^{s-1}\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(\widehat{J}_{l})}}|a_{l,\vv{c},\widehat{J}_{l}}|&\leq
3M_{\mathcal{X}}^{2}\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}p_{\min}^{-1}+M_{\mathcal{X}}^{2}\iota_{s-1}.\end{split}$
(A.9)

Let us define six events, denoted as $E_{41}$ to $E_{46}$, as follows. Here,
we recall that $\boldsymbol{X}_{iH}$ represents the feature subset in
$H\subset\\{1,\cdots,p\\}$ of the $i$th sample feature vector, which can be
represented as $(X_{ij},j\in H)^{{}^{T}}$. Our main results in this proof hold
on the event $\cap_{l=1}^{6}E_{4l}$, with the probability
$\mathbb{P}((\cap_{l=1}^{6}E_{4l})^{c})=o(1)$ proven in Lemma 5 in Section
C.5, where $A^{c}$ denotes the complementary set of $A$. The intersection
operations $\texttt{\\#}J\leq k$ below are performed over all possible subsets
$J\subset\\{1,\cdots,M\\}$ with $\texttt{\\#}J\leq k$, where we have
$\texttt{\\#}\\{J:J\subset\\{1,\cdots,M\\},\texttt{\\#}J\leq k\\}\leq M^{k}$.
We also note that
$\sum_{\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}1=2^{\texttt{\\#}\mathcal{X}(J)}\leq
2^{k\times M_{\mathcal{X}}}$ when $\texttt{\\#}J\leq k$, and that
$M_{\mathcal{X}}$ given in Condition 1 is assumed to be a finite constant.
Similar notation for union operations is employed in the proof of Lemma 5.

$\begin{split}E_{41}&=\left\\{\left|\sum_{i=1}^{n}\varepsilon_{i}f(\boldsymbol{X}_{i})\right|\leq
t_{1}\right\\},\\\ E_{42}&=\cap_{\texttt{\\#}J\leq
2,\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}\left\\{\left|\sum_{i=1}^{n}\varepsilon_{i}\boldsymbol{1}_{\boldsymbol{X}_{i\mathcal{X}(J)}=\vv{c}}\right|\leq
t_{2}\right\\},\\\ E_{43}&=\cap_{\texttt{\\#}J\leq
2,\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}\left\\{\left|\sum_{i=1}^{n}\left(f(\boldsymbol{X}_{i})\boldsymbol{1}_{\boldsymbol{X}_{i\mathcal{X}(J)}=\vv{c}}-\mathbb{E}(f(\boldsymbol{X})\boldsymbol{1}_{\boldsymbol{X}_{\mathcal{X}(J)}=\vv{c}})\right)\right|\leq
t_{3}\right\\},\\\ E_{44}&=\cap_{\texttt{\\#}J\leq
4,\vv{c}\in\\{0,1\\}^{\texttt{\\#}\mathcal{X}(J)}}\left\\{\left|\sum_{i=1}^{n}\left(\boldsymbol{1}_{\boldsymbol{X}_{i\mathcal{X}(J)}=\vv{c}}-\mathbb{P}(\boldsymbol{X}_{\mathcal{X}(J)}=\vv{c})\right)\right|\leq
t_{4}\right\\},\\\
E_{45}&=\left\\{\left|\sum_{i=1}^{n}\varepsilon_{i}^{2}\right|\leq
t_{5}\right\\},\\\
E_{46}&=\left\\{\left|\sum_{i=1}^{n}\\{(f(\boldsymbol{X}_{i}))^{2}-\mathbb{E}[(f(\boldsymbol{X}))^{2}]\\}\right|\leq
t_{6}\right\\}.\end{split}$ (A.10)

In the following, we deal with establishing an upper bound for $|D_{1,s}(J)|$
in (LABEL:equivalence.recur.4.b). For each $0\leq s\leq 2K$ and any $J$ with
$\texttt{\\#}J\leq 2$, we deduce that

$\begin{split}&\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&=\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})+\varepsilon_{i}-R_{s,J}^{\star}(\boldsymbol{X}_{i})+R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&=\left(\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}\right)+\left(\sum_{i=1}^{n}(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\right)+\sum_{i=1}^{n}\varepsilon_{i}^{2}\\\
&\qquad+2\left(\sum_{i=1}^{n}\varepsilon_{i}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))\right)+2\sum_{i=1}^{n}\varepsilon_{i}(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))\\\
&\qquad+2\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i})),\end{split}$

which, along with Jensen’s inequality, implies that for each $0\leq s\leq 2K$
and every $\texttt{\\#}J\leq 2$,

$\begin{split}|n\times
D_{1,s}(J)|&=\left|\left(\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\right)-\left(\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}\right)-\sum_{i=1}^{n}\varepsilon_{i}^{2}\right|\\\
&\leq\left(\sum_{i=1}^{n}(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\right)+2\left|\sum_{i=1}^{n}\varepsilon_{i}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))\right|\\\
&\qquad+2\sqrt{\sum_{i=1}^{n}\varepsilon_{i}^{2}}\sqrt{\sum_{i=1}^{n}(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}}\\\
&\qquad+2\sqrt{\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}}\sqrt{\sum_{i=1}^{n}(R_{s,J}^{\star}(\boldsymbol{X}_{i})-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}}.\end{split}$
(A.11)

By (A.11) and the regularity conditions assumed by Theorem 1, there exists
$L_{0}>0$ such that for each $n$ with $\frac{2\log{(n)}}{\sqrt{n}}\leq
p_{\min}$, each $0\leq s\leq 2K$, and every set $J$ with $\texttt{\\#}J\leq
2$, it holds that on the event $\cap_{l=1}^{5}E_{4l}$ with
$t_{1}=t_{2}=\log{(n)}n^{\frac{1}{2}+\frac{1}{q_{1}}}$,
$t_{3}=t_{4}=\log{(n)}\sqrt{n}$, and $t_{5}=2n\mathbb{E}\varepsilon^{2}$,

$\begin{split}&|n\times D_{1,s}(J)|\\\
&\leq\left[20p_{\min}^{-2}(\log{n})^{2}M_{\mathcal{X}}^{4}\left(\sum_{l=1}^{s}\iota_{l-1}^{2}\right)+280sM_{f}^{2}p_{\min}^{-4}(\log{n})^{2}n^{\frac{2}{q_{1}}}\right]\\\
&\qquad+2\log{(n)}n^{\frac{1}{2}+\frac{1}{q_{1}}}(3M_{\mathcal{X}}^{2}p_{\min}^{-1}\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}+M_{\mathcal{X}}^{2}\iota_{s-1}+1)\\\
&\qquad+\sqrt{20p_{\min}^{-2}(\log{n})^{2}M_{\mathcal{X}}^{4}\left(\sum_{l=1}^{s}\iota_{l-1}^{2}\right)+280sM_{f}^{2}p_{\min}^{-4}(\log{n})^{2}n^{\frac{2}{q_{1}}}}\times
2\sqrt{2n\mathbb{E}\varepsilon^{2}}\\\
&\qquad+\sqrt{20p_{\min}^{-2}(\log{n})^{2}M_{\mathcal{X}}^{4}\left(\sum_{l=1}^{s}\iota_{l-1}^{2}\right)+280sM_{f}^{2}p_{\min}^{-4}(\log{n})^{2}n^{\frac{2}{q_{1}}}}\\\
&\qquad\qquad\times
2\sqrt{n}\left(\max_{\vv{c}\in\\{0,1\\}^{p}}|f(\vv{c})|+3\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}\times
p_{\min}^{-1}+\iota_{s-1}\right)\\\ &\leq
L_{0}\times\bigg{(}s(\log{(n)})^{2}n^{\frac{2}{q_{1}}}+\left(\sum_{l=1}^{s}\iota_{l-1}^{2}\right)(\log{(n)})^{2}+n^{\frac{1}{2}+\frac{1}{q_{1}}}\log{(n)}(\iota_{s-1}+1)\sqrt{s\vee
1}\\\
&\qquad+(\iota_{s-1}+1)\sqrt{\sum_{l=1}^{s}\iota_{l-1}^{2}}\sqrt{n}\log{(n)}\bigg{)},\end{split}$
(A.12)

where the detailed derivation for the first inequality is given after the end
of the current proof of Theorem 1 (see Section B.1.2), and the the last
inequality is due to subadditivity of the square root function.

Let us simplify the upper bound on the RHS of (LABEL:equivalence.recur.10.b).
When $n^{-\frac{1}{2}+\frac{1}{q_{1}}}\log{(n)}\sqrt{s\vee 1}\leq 1$ and
$\sqrt{\sum_{l=1}^{s}\iota_{l-1}^{2}}(\log{n})n^{-1/2}\leq 1$, under the same
conditions for (LABEL:equivalence.recur.10.b),

$\begin{split}|D_{1,s}(J)|&\leq
L_{0}n^{-\frac{1}{2}+\frac{1}{q_{1}}}\log{(n)}\sqrt{s\vee
1}(\iota_{s-1}+2)+L_{0}(\iota_{s-1}+2)\sqrt{\sum_{l=1}^{s}\iota_{l-1}^{2}}n^{-\frac{1}{2}}\log{(n)}\\\
&=L_{0}n^{-\frac{1}{2}}\log{(n)}\left(n^{\frac{1}{q_{1}}}\sqrt{s\vee
1}+\sqrt{\sum_{l=1}^{s}\iota_{l-1}^{2}}\right)(\iota_{s-1}+2).\end{split}$
(A.13)

Next, proceed to deal with the term $D_{2,s}(J)$’s on the RHS of
(LABEL:equivalence.recur.4.b). We first deduce that

$\begin{split}D_{2,s}(J)&=\left(n^{-1}\sum_{i=1}^{n}(f(\boldsymbol{X}_{i})-R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}\right)-\mathbb{E}(f(\boldsymbol{X})-R_{s,J}^{\star}(\boldsymbol{X}))^{2}\\\
&=n^{-1}\sum_{i=1}^{n}[(f(\boldsymbol{X}_{i}))^{2}-\mathbb{E}(f(\boldsymbol{X}))^{2}]\\\
&\qquad-2n^{-1}\sum_{i=1}^{n}[f(\boldsymbol{X}_{i})R_{s,J}^{\star}(\boldsymbol{X}_{i})-\mathbb{E}(f(\boldsymbol{X})R_{s,J}^{\star}(\boldsymbol{X}))]\\\
&\qquad+n^{-1}\sum_{i=1}^{n}[(R_{s,J}^{\star}(\boldsymbol{X}_{i}))^{2}-\mathbb{E}(R_{s,J}^{\star}(\boldsymbol{X}))^{2}].\end{split}$
(A.14)

By (A.14) and the assumed regularity conditions, there exists some
sufficiently large $L_{0}>0$ such that on the event $E_{43}\cap E_{44}\cap
E_{46}$ defined in (A.10) with $t_{3}=t_{4}=t_{6}=\log{(n)}\sqrt{n}$, it holds
that for each $0\leq s\leq 2K$ and every set $J$ with $\texttt{\\#}J\leq 2$,

$\begin{split}|D_{2,s}(J)|&\leq\log{(n)}n^{-\frac{1}{2}}+2\log{(n)}n^{-\frac{1}{2}}(3M_{\mathcal{X}}^{2}p_{\min}^{-1}\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}+M_{\mathcal{X}}^{2}\iota_{s-1})\\\
&\qquad+\log{(n)}n^{-\frac{1}{2}}(3M_{\mathcal{X}}^{2}p_{\min}^{-1}\sqrt{\mathbb{E}[(f(\boldsymbol{X}))^{2}]}+M_{\mathcal{X}}^{2}\iota_{s-1})^{2}\\\
&\leq L_{0}\log{(n)}n^{-\frac{1}{2}}(\iota_{s-1}^{2}+1).\end{split}$ (A.15)

where the detailed derivation for the first inequality on the RHS of (A.15) is
given after the end of the current proof of the Theorem 1 (see Section B.1.3),
and the second inequality is due to simple calculations and an application of
the quadratic mean inequality [21].

To summarize the results from (LABEL:equivalence.recur.4.b), (A.13), and
(A.15), under the assumed conditions, on the event $\cap_{l=1}^{6}E_{4l}$, it
holds that for each $0\leq s\leq 2K$, and every set $J$ with
$\texttt{\\#}J\leq 2$,

$\begin{split}|D_{1,s}(J)|&\leq
L_{0}n^{-\frac{1}{2}}(\log{n})\left(n^{\frac{1}{q_{1}}}\sqrt{s\vee
1}+\sqrt{\sum_{l=1}^{s}\iota_{l-1}^{2}}\right)(\iota_{s-1}+2),\\\
|D_{2,s}(J)|&\leq
L_{0}(\log{n})n^{-\frac{1}{2}}(\iota_{s-1}^{2}+1).\end{split}$

Particularly, for each $0<s\leq 2K$ and every set $J$ with $\texttt{\\#}J\leq
2$,

$|D_{2,s-1}(\widehat{J}_{s-1})|+|D_{2,s}(J)|+|D_{1,s-1}(\widehat{J}_{s-1})|+|D_{1,s}(J)|\leq
d_{n},$ (A.16)

where
$d_{n}=L_{1}n^{-\frac{1}{2}}(\log{n})(n^{\frac{1}{q_{1}}}\sqrt{2K}+\sqrt{\sum_{l=1}^{2K}\iota_{l-1}^{2}})(\iota_{2K-1}^{2}+\iota_{2K-1}+3)$
for sufficiently large $L_{1}\geq 1$. Additionally,
$\mathbb{P}((\cap_{l=1}^{6}E_{4l})^{c})=o(1)$ by Lemma 5. We note that the
assumption $d_{n}\leq 1$ implies the conditions
$\sqrt{\sum_{l=1}^{2K}\iota_{l-1}^{2}}(\log{n})n^{-1/2}\leq 1$ and
$n^{-\frac{1}{2}+\frac{1}{q_{1}}}(\log{n})\sqrt{2K}\leq 1$. These conditions,
along with the requirement $\frac{2\log{n}}{\sqrt{n}}\leq p_{\min}$, are
required for (A.13) to hold.

The results (LABEL:equivalence.recur.4.b) and (A.16) implies that for each
$0<s\leq 2K$ and every set $J$ with $\texttt{\\#}J\leq 2$,

$\begin{split}&\Big{|}\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&\qquad-\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}))^{2}]\Big{|}\\\
&\leq d_{n}.\end{split}$ (A.17)

In light of (LABEL:equivalence.recur.4.c), we now show that
$\widehat{J}_{1},\cdots,\widehat{J}_{K}$ include all significant component in
$S_{1}^{*}$ such that $S_{1}^{*}\subset\cup_{l=1}^{K}\widehat{J}_{l}$. By
(LABEL:equivalence.recur.4.c), for $0<s\leq K$, and those insignificant
components $m\in\\{1,\cdots,M\\}\backslash S_{1}^{*}$,

$\begin{split}&\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,\\{m\\}}(\boldsymbol{X}_{i}))^{2}\\\
&\leq\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(m)}))^{2}]+d_{n}.\end{split}$
(A.18)

Since we have assumed that Algorithm 1 makes the initial $2K$ selections
without any repetitions, we do not need to consider the cases with
$\widehat{J}_{s}\subset\cup_{l=1}^{s-1}\widehat{J}_{l}\cap S_{1}^{*}$ when
deciding $\widehat{J}_{s}$ for $0<s\leq K$. Meanwhile, by
(LABEL:equivalence.recur.4.c), for $0<s\leq K$, and those significant
components that have not been selected $m\in
S_{1}^{*}\backslash(\cup_{l=1}^{s-1}\widehat{J}_{l})$,

$\begin{split}&\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,\\{m\\}}(\boldsymbol{X}_{i}))^{2}\\\
&\geq\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(m)}))^{2}]-d_{n}.\end{split}$
(A.19)

From (LABEL:bias.control.1) and (LABEL:bias.control.2), it is seen that the
analysis of
$\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(m)}))^{2}]$
is crucial for our desired results. We can show that for each
$m\in\\{1,\cdots,M\\}\backslash(\cup_{l=1}^{s-1}\widehat{J}_{l})$ and each
$0<s\leq K$,

$\begin{split}&\left|\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(m)}))^{2}]-\textnormal{Var}(g_{m}(\boldsymbol{X}))\right|\leq\frac{21\delta_{0}}{(p_{\min})^{2}}\mathbb{E}[(f(\boldsymbol{X}))^{2}],\end{split}$
(A.20)

whose proof is deferred to the end of the current proof of Theorem 1 (see
Section B.1.4).

By (A.4), (LABEL:bias.control.1)–(LABEL:bias.control.3.b), Condition 2, and
other regularity conditions assumed by Theorem 1, it holds that on the event
$\cap_{l=1}^{6}E_{4l}$ with $d_{n}\leq 1$ and $\frac{2\log{(n)}}{\sqrt{n}}\leq
p_{\min}$,

$S_{1}^{\star}\subset\cup_{1\leq l\leq K}\widehat{J}_{l}.$ (A.21)

We have finished the proof of the results concerning $0<s\leq K$, and now we
proceed to deal with the case with $K<s\leq 2K$. By
(LABEL:equivalence.recur.4.c), for $K<s\leq 2K$, and those insignificant
components $J\subset\\{1,\cdots,M\\}$ with $\texttt{\\#}J=2$ and $J\not\in
S_{2}^{*}$,

$\begin{split}&\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&\leq\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}))^{2}]+d_{n}.\end{split}$
(A.22)

Since we have assumed that Algorithm 1 makes the initial $2K$ selections
without any repetitions, we do not need to consider those has already been
selected when deciding $\widehat{J}_{s}$. Meanwhile, by
(LABEL:equivalence.recur.4.c), for $K<s\leq 2K$, and those significant
components that have not been selected $J\in
S_{2}^{*}\backslash\\{\widehat{J}_{K+1},\cdots,\widehat{J}_{s-1}\\}$,

$\begin{split}&\left(n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s-1}(\boldsymbol{X}_{i}))^{2}\right)-n^{-1}\sum_{i=1}^{n}(Y_{i}-\widehat{U}_{s,J}(\boldsymbol{X}_{i}))^{2}\\\
&\geq\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}))^{2}]-d_{n}.\end{split}$
(A.23)

Let us analyze the term
$\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)})$
on the RHS of (LABEL:bias.control.1.2k) and (LABEL:bias.control.2.2k) as
follows. With the conditions required for (A.21) and Condition 3, it holds
that on the event $\cap_{l=1}^{6}E_{4l}$, for each $K<s\leq
K+\texttt{\\#}S_{2}^{*}$, each $J\subset\\{1,\cdots,M\\}$ with
$\texttt{\\#}J=2$ and
$J\not\in\\{\widehat{J}_{K+1},\cdots,\widehat{J}_{s-1}\\}$,

$\begin{split}&\left|\mathbb{E}[(\mathbb{E}(f(\boldsymbol{X})-R_{s-1}^{\star}(\boldsymbol{X})|\boldsymbol{X}_{\mathcal{X}(J)}))^{2}]-\textnormal{Var}\left[g_{J}(\boldsymbol{X})-\left(\sum_{j\in
J}g_{\\{j\\}}(\boldsymbol{X})\right)\right]\right|\\\ &\leq