Fr\'echet random forests for metric space valued regression with non euclidean predictors

TL;DR

This paper introduces Fréchet trees and forests, extending random forests to handle data in general metric spaces, including complex data types like curves, images, and shapes, with theoretical guarantees and practical applications.

Contribution

It develops a novel splitting and prediction framework for random forests in metric spaces, enabling analysis of heterogeneous data types with theoretical consistency results.

Findings

Effective handling of heterogeneous data including images and curves.

Theoretical consistency of the Fréchet regressogram predictor.

Successful application to air quality data.

Abstract

Random forests are a statistical learning method widely used in many areas of scientific research because of its ability to learn complex relationships between input and output variables and also its capacity to handle high-dimensional data. However, current random forest approaches are not flexible enough to handle heterogeneous data such as curves, images and shapes. In this paper, we introduce Fr\'echet trees and Fr\'echet random forests, which allow to handle data for which input and output variables take values in general metric spaces. To this end, a new way of splitting the nodes of trees is introduced and the prediction procedures of trees and forests are generalized. Then, random forests out-of-bag error and variable importance score are naturally adapted. A consistency theorem for Fr\'echet regressogram predictor using data-driven partitions is given and applied to Fr\'echet purely uniformly random trees. The method is studied through several simulation scenarios on heterogeneous data combining longitudinal, image and scalar data. Finally, one real dataset about air quality is used to illustrate the use of the proposed method in practice.

Tables

Table 1. Table 1: Random draws in the MNIST dataset of the output images from the realizations G i 1 , G i 2 subscript superscript 𝐺 1 𝑖 subscript superscript 𝐺 2 𝑖 G^{1}_{i},\ G^{2}_{i} and β i subscript 𝛽 𝑖 \beta_{i} used to simulate the input curves.

	$G_i^1=0,G_i^2=0$	$G_i^1=1,G_i^2=0$	$G_i^1=0,G_i^2=1$	$G_i^1=1,G_i^2=1$
${\beta }_i>1$	${ℳ}_i^1$	${ℳ}_i^3$	${ℳ}_i^5$	${ℳ}_i^7$
${\beta }_i\le 1$	${ℳ}_i^2$	${ℳ}_i^4$	${ℳ}_i^6$	${ℳ}_i^8$

Equations

$$\begin{array}[]{ccccc}\mathcal{F}_{n}&:&\mathcal{Z}&\longmapsto&\mathbb{R}^{+}\\ &&z&\longmapsto&\frac{1}{n}\overset{n}{\underset{i=1}{\sum}}d^{2}(z,z_{i})\\ \end{array}$$

$${\overline{z}}_n\in \underset{z\in \mathcal{Z}}{\mathrm{argmin}}{\mathcal{F}}_n(z)$$

$${\mathcal{V}}_n={\mathcal{F}}_n({\overline{z}}_n)=\frac{1}{n}\sum _{i=1}^n{d}^2(z_i,{\overline{z}}_n)$$

$$H_{n,j}(A,c_{j,l},c_{j,r})={\mathcal{V}}_n(A)-\frac{N_n\left(A_{j,r}\right)}{N_n\left(A\right)}{\mathcal{V}}_n(A_{j,r})-\frac{N_n\left(A_{j,l}\right)}{N_n\left(A\right)}{\mathcal{V}}_n(A_{j,l})$$

$$\mathcal{V}_{n}(A)=\frac{1}{N_{n}(A)}\sum_{i:X_{i}\in A}d_{\mathcal{Y}}^{2}(Y_{i},\overline{Y}_{A})\qquad\mbox{(resp. for $A_{j,l}$ and $A_{j,r}$).}$$

$$\overline{Y}_{A}=\operatorname*{arg\,min}_{y\in\mathcal{Y}}\sum_{i:X_{i}\in A}d_{\mathcal{Y}}^{2}(y,Y_{i})\qquad\mbox{(resp. for $\overline{Y}_{A_{j,l}}$ and $\overline{Y}_{A_{j,r}}$).}$$

$$j_n^{\ast }=\underset{j\in \{1,\dots ,p\}}{\mathrm{argmax}}{H}_{n,j}.$$

$${\overline{Y}}_t=\underset{y\in \mathcal{Y}}{\mathrm{argmin}}\sum _{i:X_i\in t}{d}_{\mathcal{Y}}^2\left(y,Y_i\right)$$

$$\mathrm{err}(T_n(x))=d_{\mathcal{Y}}^2(T_n(x),y).$$

$$\hat{y}=\underset{z\in \mathcal{Y}}{\mathrm{argmin}}\sum _{l=1}^q{d}_{\mathcal{Y}}^2(z,T_n(x,{\Theta }_l,{\Theta }_l^{\prime })).$$

$$\mathrm{errOOB}=\frac{1}{n}\sum _{i=1}^n{d}_{\mathcal{Y}}^2(Y_i,{\hat{Y}}_i^{OOB}).$$

$$VI(X^{(j)})=\frac{1}{q}\sum _{l=1}^q\left(\mathrm{err}{\widetilde{\mathrm{OOB}}}_l^j-errOO{B}_l\right).$$

$$L\left(y,\varphi (x)\right)=d^2\left(y,\varphi (x)\right)y\in \mathcal{Y},x\in \mathcal{X}$$

$$R(\varphi )=\mathbb{E}\left[L(Y,\varphi \left(X\right)\right]=\mathbb{E}\left[d^2(Y,\varphi \left(X\right)\right]$$

$${\varphi }^{\ast }\in \underset{\varphi :\mathcal{X}⟶\mathcal{Y}}{\mathrm{argmin}}R(\varphi )$$

$$R(\varphi )={\mathbb{E}}_X\left({\mathbb{E}}_Y\left[d^2(Y,X)|X\right]\right)$$

$$r\left(x,\varphi (x)\right)={\mathbb{E}}_Y\left[L(Y,\varphi (X))|X=x\right]={\mathbb{E}}_Y\left[d^2(Y,\varphi (X))|X=x\right]$$

$${\varphi }^{\ast }(x)\in \underset{y\in \mathcal{Y}}{\mathrm{argmin}}r\left(x,y\right),\text{mbox}wherer\left(x,y\right)={\mathbb{E}}_Y\left[d^2(Y,y)|X=x\right].$$

$${\mathcal{A}}_n:=\{{\pi }_n\left(z_1,\dots ,z_n\right),\left(z_1,\dots ,z_n\right)\in {\mathcal{Z}}^n\}$$

$$\left\{A_1\cap B,A_2\cap B,\dots ,A_n\cap B\right\}$$

$${\Delta }_n^{\ast }\left(\mathcal{A}\right)=\underset{x_1^n\in {\mathbb{R}}^{d.n}}{\max }\Delta \left(\mathcal{A},x_1^n\right)$$

$$\mathbb{P}\left(\underset{\pi \in \mathcal{A}}{\sup }\sum _{A\in \pi }\left|\mu (A)-{\mu }_n(A)\right|>ϵ\right)\le 4{\Delta }_n^{\ast }\left(\mathcal{A}\right)2^{\mathcal{C}\left(\mathcal{A}\right)}\exp \left(-n{ϵ}^2/32\right)$$

$$T_n\left(x\right)=\underset{y\in \mathcal{Y}}{\mathrm{argmin}}\frac{1}{n}\sum _{i=1}^n{d}^2(Y_i,y)1\{X_i\in {\pi }_n[x]\}$$

$$N(ϵ,\mathcal{Y},d)\le {\left(\frac{2\mathrm{diam}(\mathcal{Y})}{ϵ}\right)}^{\mathrm{ddim}(\mathcal{Y})}$$

$$\underset{n\to \infty }{\lim }\left|r\left(x,T_n(x)\right)-\underset{y\in \mathcal{Y}}{\min }r\left(x,y\right)\right|=0,\text{mbox}a.s.$$

$$\underset{n\to \infty }{\lim }R(T_n)-R({\varphi }^{\ast })=0,\text{mbox}a.s$$

$$FPUR{T}_n(x)=\underset{y\in \mathcal{Y}}{\mathrm{argmin}}\frac{1}{n}\sum _{j=0}^{k_n}\sum _{i=1}^n{d}^2(y,Y_i)1\{U_{(j)}\le x\le U_{(j+1)}\}\forall x\in [0,1]$$

$$\underset{n\to \infty }{\lim }R(FPUR{T}_n)-R({\varphi }^{\ast })=0,\text{mbox}a.s$$

$$\frac{\mathcal{C}({\Pi }_n)}{n}=\frac{k_n+1}{n}\underset{n⟶\infty }{\to }0$$

$$\frac{\log {\Delta }_n^{\ast }({\Pi }_n)}{n}\le \frac{k_n\log n}{n}\underset{n\to \infty }{⟶}0$$

Full text

K****eywords Random forests  $\cdot$ Nonparametric regression  $\cdot$ Metric spaces regression  $\cdot$ Longitudinal data  $\cdot$ heterogeneous data

1 Introduction

Random Forests [Breiman, 2001] are one of the state-of-the-art machine learning methods. It owes its success to very good predictive performance coupled with very few parameters to tune. Moreover, as a tree-based method, it is able to handle regression and classification (2-class or multi-class) problems in a consistent manner and deals with quantitative or qualitative input variables. Finally, its non-parametric nature allows to proceed high-dimensional data where the number of input variables is very large in regards of statistical units.

The general principle of a tree predictor is to recursively partition the input space. Starting from the root node which contains all learning observations, it repeatedly splits each node into two or more child nodes until a stopping rule is reached. Let us focus, for the sake of clarity, on the case where all input variables are quantitative. For most of tree predictors, splits are binary and consist in an input variable $X^{j}$ and a threshold $s$, leading to two child nodes containing observations that verify $\{X^{j}\leq s\}$ and $\{X^{j}>s\}$ respectively [Breiman et al., 1984]. The splitting variable as well as the threshold are most of the time sought to minimize an heterogeneity criterion on child nodes. The main idea is to partition the input space into more and more homogeneous regions in terms of the output variable.

One limitation of the previously described splitting strategy is that all input variables must live in an ordered space (the method must decide if an observation of the splitting variable is less or larger than the threshold). Yet, with complex data structures, inputs can belong to unordered spaces. As an illustrative example, the real datasets to be analyzed in this paper are from HIV vaccine trials made of input and output variables which are discretely sampled curves representing repeated measurements over time. This example typically corresponds to observations collected in longitudinal studies. In such settings, a main objective is to predict for a given individual the output trajectory given the knowledge of inputs trajectories. If this objective is reached at the trajectory level, then the notion of order between explanatory variables is lost. However, ignoring the fact that measurements are repeated observations over time generally leads to an important loss of information for prediction. Thus, one way of analyzing this kind of data is to generalize the notion of split in unordered metric spaces. Recently, random forests have been adapted to the general metric space framework but in the special case where neither the representation of the data nor the distances between data points are available [Haghiri et al., 2018]. In the present work, we take into account the distances between any items of the space.

Hence, we consider the framework of a learning sample $\mathcal{L}_{n}=\{(X_{1},Y_{1}),\ldots,(X_{n},X_{n})\}$ made of i.i.d. observations of a pair of random variable $(X,Y)\in\mathcal{X}\times\mathcal{Y}$, where $\mathcal{X}$ is a product of $p$ metric spaces (which can be unordered) $(\mathcal{X}_{1},d_{1})\times\cdots\times(\mathcal{X}_{p},d_{p})$, and where $\mathcal{Y}$ is also a metric space with distance $d_{\mathcal{Y}}$. The core idea of this work is to generalize the notion of split that only uses the distance of each metric space. Furthermore, as the notion of mean in the output space $\mathcal{Y}$ is also needed to allocate predictions to terminal nodes of a tree, the Fréchet mean that generalizes the mean in general metric spaces [Fréchet, 1906] is used. This justifies the names Fréchet trees and Fréchet random forests hereafter. Once the notion of split is defined, the building of a maximal tree and the pruning of that tree to obtain an optimal tree are extended. Finally, with this generalization of CART trees, Fréchet random forests are derived in a rather standard way: a forest predictor is an aggregation of a collection of randomized trees. In our framework, the aggregation step therefore consists in taking the Fréchet mean of individual tree predictions.

The use of the Fréchet mean has now become a standard tool for statistical inference from manifold-valued data. For example, it is the key notion allowing to perform PCA for non-Euclidean data such as functional data on Riemanian manifolds (see e.g. [Dai and Muller, 2018], [Fletcher et al., 2004], [Sommer et al., 2010]) or histograms [Cazelles et al., 2018], and to analyze ensemble of complex objects with their shape, such as ECG curves [Bigot, 2013] or phylogenetic trees [Nye et al., 2017]. New innovative regression methods have also emerged for the framework of a metric space valued output variable with Euclidean predictors [Petersen and Müller, 2019]. The methods proposed in this paper allow to perform nonparametric regression between predictors taking their values in different metric spaces and a metric space valued output.

In this paper, we first present the Fréchet tree predictor in Section 2 before introducing Fréchet random forests in Section 3. We introduce an extremely randomized version of the Fréchet random forests method in Section 3.3. Section 4 is dedicated to the analysis of the consistency of Fréchet regressogram estimators using data-driven partitions with output lying in metric space. We report numerical experiments using simulated longitudinal data to compare our approach with competitive methods, then we analyze two scenarios of heterogeneous data simulations involving curves, images and scalars in Section 5. An application of Fréchet random forests for statistical inference from longitudinal data is presented in Section 6. Finally, we discuss in Section 7 potential extensions of this work. All the numerical experiments of this paper are reproducible from our R package FrechForest111https://github.com/Lcapitaine/FrechForest.

2 Fréchet Trees

2.1 Fréchet means and Fréchet variance

The notions of mean and variance are central to the construction of regression trees [Breiman et al., 1984]. We introduce in this section the notions of Fréchet empirical mean and Fréchet empirical variance [Fréchet, 1948], which are the natural generalization of mean and variance in metric spaces. Let $(z_{1},\ldots,z_{n})$ a sample from a metric space $\left(\mathcal{Z},d\right)$, the empirical Fréchet function is given by

[TABLE]

the function $\mathcal{F}_{n}(z)$ measures the average squared distance between $z\in\mathcal{Z}$ and $z_{1},\ldots,z_{n}$. We define the empirical Fréchet means $\overline{z}_{n}$ of the sample $(z_{1},\ldots,z_{n})$ as any minimizer of the empirical Fréchet function, i.e.

[TABLE]

Note that the Fréchet mean can be non unique. The empirical Fréchet variance $\mathcal{V}_{n}$ of the sample $(z_{1},\ldots,z_{n})$ is then given by

[TABLE]

Note that even if the empirical Fréchet mean may not be a unique element of the metric space, the Fréchet variance is unique. Throughout the paper, Fréchet’s mean and Fréchet variance will always refer to Fréchet empirical mean and Fréchet’s empirical variance. For the sake of simplicity, we assume in the rest of the paper that the Fréchet mean is unique.

2.2 Splitting rule

One key ingredient in the building of a decision tree is the way its nodes are split [Breiman et al., 1984]. Splitting a node $t$ of a tree according to some variable $X^{(j)}$ amounts to find a way of grouping observations of this node into two subsets constituting the child nodes. This grouping is usually performed to maximize the differences between the two resulting child nodes according to the output variable. However, if variable $X^{(j)}$ is strongly related to the output variable $Y$, then it is expected that for two observations with "close" $X^{(j)}$ values in $(\mathcal{X}_{j},d_{j})$, associated outputs will be "close" in $(\mathcal{Y},d_{\mathcal{Y}})$. From this idea, we introduce a way of splitting nodes in general metric spaces.  Let $\left(\mathcal{Z},d\right)$ be a metric space, a split is any couple of distinct elements $\left(c_{1},c_{2}\right)$ of $\mathcal{Z}$. We define the partition $\mathcal{P}=\left\{P_{1},P_{2}\right\}$ associated with elements $\left(c_{1},c_{2}\right)$ by $P_{1}=\left\{z\in\mathcal{Z},d\left(z,c_{1}\right)\leq d\left(z,c_{2}\right)\right\}$ and $P_{2}=\left\{z\in\mathcal{Z},d\left(z,c_{2}\right)<d\left(z,c_{1}\right)\right\}$.

Let $A$ be a subset of the input space $\mathcal{X}$ and for any $j=1,\ldots,p$, let $A_{j}=\{x^{(j)},x=\left(x^{(1)},\ldots,x^{(p)}\right)\in A\}$ denotes the set of the $j$-th coordinates of the components of A. Let $\left(c_{j,l},c_{j,r}\right)$ be a split on $(A_{j},d_{j})$, denote $A_{j,r}$ and $A_{j,l}$ the right and left child nodes (i.e. the associated partition) obtained from the split $\left(c_{j,l},c_{j,r}\right)$.

The quality of the split $\left(c_{j,l},c_{j,r}\right)$ is then defined by the following measure of Fréchet variance decrease:

[TABLE]

where $N_{n}(A)$ is the number of observations of the learning set $\mathcal{L}_{n}$ belonging to $A$ and $\mathcal{V}_{n}(A),\ \mathcal{V}_{n}(A_{j,l})$ and $\mathcal{V}_{n}(A_{j,r})$ are the empirical Fréchet variances of outputs in $A$, $A_{j,r}$ and $A_{j,l}$ i.e.

[TABLE]

$\overline{Y}_{A}$, $\overline{Y}_{A_{j,l}}$ and $\overline{Y}_{A_{j,r}}$ are the Fréchet means of outputs associated to observations belonging to nodes $A$, $A_{j,l}$ and $A_{j,r}$ i.e.

[TABLE]

It is worth noting that the decrease in Fréchet variance for each possible split is compared with the output space metric, which makes it possible to compare splits made on input variables from different metric spaces. At last, the split variable $j_{n}^{*}$, chosen for splitting the node is the one that maximizes $H_{n,j}$, that is

[TABLE]

It is easy to show that $H_{n,j_{n}^{*}}\geq 0$ for all $n$, thanks to the use of the Fréchet mean, which means that each split leads to a decrease of Fréchet variance.

To determine the successive splits $(c_{j,l},c_{j,r})$, the user defines, in a preliminary step, a split function i.e. a way to find the two representatives $(c_{j,l},c_{j,r})$. More precisely, a split function is an application which associates a couple $\left(c_{1},c_{2}\right)\in\mathcal{Z}$ to any sample $\{h_{1},\ldots,h_{n}\}$ from a general metric space $\left(\mathcal{Z},d\right)$. For example, the $2$-means algorithm ($k$-means with $k=2$) can be used to determine the representatives. Note that for each metric space $\left(\mathcal{X}_{j},d_{j}\right)$ we can use a different split function.

2.3 Tree building

Starting from the root node (associated with the whole input space $\mathcal{X}$), nodes are recursively split in order to give a partition of the input space $\mathcal{X}$. A node $t$ of the tree is not split if it is pure, that is if the Fréchet variance of this node is null. As a first step in the building process, the tree is developed until all nodes are pure, leading to the so-called maximal tree. Then, the pruning algorithm of CART [Breiman et al., 1984] is applied, with the use of the Fréchet variance instead of the standard empirical variance. At the end of this step, a sequence of nested sub-trees of the maximal tree is obtained. Next, the sub-tree associated to the lowest prediction error (estimated by cross-validation) is selected as the final tree predictor. The way a Fréchet tree predicts new inputs is detailed in the next section.

2.4 Prediction

Let $T_{n}$ be a Fréchet tree, we note $\widetilde{T}_{n}$ the set of leaves (i.e., terminal nodes) of $T_{n}$. For each leaf $t\in\widetilde{T}_{n}$, the Fréchet mean of the outputs of observations belonging to $t$ is associated to $t$. Then the prediction of the output variable associated with any $x\in\mathcal{X}$ is given by $\widehat{y}=T_{n}(x)=\sum_{t\in\widetilde{T}_{n}}\overline{Y}_{t}\mathbb{1}_{x\in t}$, where $\mathbb{1}_{S}$ denotes the indicator function of a set $S$ and $\overline{Y}_{t}$ is the Fréchet mean of outputs in $t$

[TABLE]

In order to determine to which leaf belongs an observation $x$, it is dropped down the tree as follows. Starting from the root node, the associated split variable $X^{(j_{1})}$ is considered, together with its two child nodes $A_{j_{1},l}$ and $A_{j_{1},r}$, as well as the corresponding representatives $c_{j_{1},l}$ and $c_{j_{1},l}$. To decide in which child node $x$ must fall, its $d_{j_{1}}$-distance with $c_{j_{1},l}$ and $c_{j_{1},l}$ must be computed and $x$ goes to $A_{j_{1},l}$ if $d_{j_{1}}(x^{(j_{1})},c_{j_{1},l})<d_{j_{1}}(x^{(j_{1})},c_{j_{1},r})$ and to $A_{j_{1},r}$ otherwise. This process is then repeated until $x$ falls into a leaf. The error made by $T_{n}$ on $x$ is defined as:

[TABLE]

3 Fréchet random forests

3.1 An aggregation of Fréchet trees

A Fréchet random forest is derived as standard random forests [Breiman, 2001]: it consists in an aggregation of a collection of randomized Fréchet trees. Here, the same random perturbations as standard random forests [Breiman, 2001] are used. Let $l\in\{1,\ldots,q\}$ and consider the $l$-th tree. First, it is built on a bootstrap sample of the learning sample $\mathcal{L}_{n}^{\Theta_{l}}$ ($n$ observations drawn with replacement among $\mathcal{L}_{n}$), and secondly, the search for the optimized split for each node is restricted to a subset of $mtry$ variables randomly drawn among the $p$ input variables (this random subset is denoted $\Theta_{l}^{\prime}$ hereafter). Hence, the $l$-th tree is denoted $T_{n}\left(.,\Theta_{l},\Theta_{l}^{\prime}\right)$ and can be viewed as a doubly-randomized Fréchet tree. Once all randomized trees are built, the Fréchet mean is again used to aggregated them. Thus, for any $x\in\mathcal{X}$ the prediction made by the Fréchet random forest is:

[TABLE]

3.2 OOB error and variable importance scores

Fréchet random forests inherit from standard random forest quantities: OOB (Out-Of-Bag) error and variable importance scores. The OOB error provides a direct estimation of the prediction error of the method and proceeds as follows. The predicted output value, $\widehat{Y}_{i}^{OOB}$, of the $i$-th observation $(X_{i},Y_{i})\in\mathcal{L}_{n}$, is obtained by aggregating only trees built on bootstrap samples that do not contain $(X_{i},Y_{i})$. The OOB error is then computed as the average squared distance between those predictions and the $Y_{i}$:

[TABLE]

Variable importance (VI) provides information on the use of input variables in the learning task that can be used e.g. to perform variable selection. For $j\in\{1,\ldots,p\}$, variable importance of input variable $X^{(j)}$, $VI(X^{(j)})$, is computed as follows. For the $l$-th bootstrap sample $\mathcal{L}_{n}^{\Theta_{l}}$, let us define the associated $\mathrm{OOB}_{l}$ sample of all observations that were not picked in $\mathcal{L}_{n}^{\Theta_{l}}$. First, $\mathrm{errOOB}_{l}$, the error made by tree $T_{n}\left(.,\Theta_{l},\Theta_{l}^{\prime}\right)$ on $\mathrm{OOB}_{l}$ is computed. Then, the values of $X^{(j)}$ in the $\mathrm{OOB}_{l}$ sample are randomly permuted, to get a disturbed sample $\widetilde{\mathrm{OOB}}_{l}^{j}$, and the error, $\mathrm{err}\widetilde{\mathrm{OOB}}_{l}^{j}$, made by $T_{n}\left(.,\Theta_{l},\Theta_{l}^{\prime}\right)$ on $\widetilde{\mathrm{OOB}}_{l}^{j}$ is calculated. Finally, VI of $X^{(j)}$ is defined as:

[TABLE]

3.3 Extremely randomized Fréchet random forests

The construction of a Fréchet tree is conditioned by: i) the existence of the Fréchet mean for the output space $\left(\mathcal{Y},d_{\mathcal{Y}}\right)$; ii) the use of a calculable split function for each input space. As mentioned in Section 2.2, in practice the 2-means algorithm can be used as the split function. However, it may not be applicable on all input spaces, for example on input spaces where the Fréchet mean does not exist. In order to have a split function applicable on all input metric spaces, we use the split function introduced by [Geurts et al., 2006] for regression and classification trees in $\mathbb{R}^{p}$: let ntry be an integer between 1 and $n(n-1)/2$, we randomly draw ntry different splits i.e. ntry different couples of representatives, then we calculate the reduction of the Fréchet variance associated to each of these splits for the response variable and finally we select the split which maximizes the reduction of the Fréchet variance on the response variable. An extremely randomized Fréchet tree (ERFT) is any tree built with this random split function. An aggregation of extremely randomized Fréchet trees is called an extremely randomized Fréchet random forest (ERFRF). Note that when $\texttt{ntry}=n(n-1)/2$ the node split is no longer random. This splitting strategy has two advantages: it is applicable for any type of input and by taking a low value of ntry, it allows to drastically reduce calculation times while having excellent prediction capabilities (see Section 5.4.1).

4 Theory

In this section we study the consistency of Fréchet regressogram using data-driven partitions. First, we recall the notions of specific risk and global risk in a general framework before recalling the notion of Fréchet function. Then we remind the notion of family of partitions on $\mathbb{R}^{p}$. Finally we give the definition of Fréchet regressogram using data-driven partition and a result of its consistency in the case where the input space is $\mathbb{R}^{p}$ and the output space is a metric space.

4.1 Problem

In this section we present some notations in the general framework where $\mathcal{X}$ is any separable space and $(\mathcal{Y},d)$ is a sparable metric space. Consider the pair of random variables $\left(X,Y\right)\in\mathcal{X}\times\left(\mathcal{Y},d\right)$. The task is to learn a mapping $\phi:\mathcal{X}\longrightarrow\mathcal{Y}$.

For any mapping $\phi:\mathcal{X}\longrightarrow\mathcal{Y}$ the loss function $L$ is given by

[TABLE]

The global risk associated with the mapping $\phi$ is defined by

[TABLE]

The bayes optimal mapping $\phi^{*}$ is any minimizer of the global risk function i.e.

[TABLE]

When $\mathcal{X}$ and $\mathcal{Y}$ are separable, according to [Blackwell and Maitra, 1984] the global risk can be factorized as

[TABLE]

We define the point risk function of $\phi$ by

[TABLE]

The Bayes optimal point-risk mapping $\phi^{*}$ is defined by

[TABLE]

This mapping introduced in [Petersen and Müller, 2019] is called Fréchet regression function.

4.2 Family of partitions

Let $\mathcal{X}=\mathbb{R}^{p}$, denote $\mathcal{Z}=\mathbb{R}^{p}\times\mathcal{Y}$ and let $\pi_{n}$ be a partitioning rule of $\mathbb{R}^{p}$ i.e a function that associates a measurable partition of $\mathbb{R}^{p}$ to any vector $\left(z_{1},\ldots,z_{n}\right)\in\mathcal{Z}^{n}$. We note $\mathcal{A}_{n}$ the family of all the partitions we can obtain with $\pi_{n}$:

[TABLE]

We denote $\mathcal{C}\left(\mathcal{A}_{n}\right)=\underset{\pi\in\mathcal{A}_{n}}{\sup}|\pi|$ the maximal number of cells for the partitions family $\mathcal{A}_{n}$. Finally, let $\mathcal{A}$ be a family of partitions, let $x_{1},\ldots,x_{n}$ $n$ points of $\mathbb{R}^{p}$ and let $B=\left\{x_{1},\ldots,x_{n}\right\}$. We note $\Delta\left(\mathcal{A},x_{1}^{n}\right)$ the number of distinct partitions

[TABLE]

induced by the partitions $\{A_{1},\ldots,A_{n}\}\in\mathcal{A}$. The growing function of the partitions family $\mathcal{A}$ is defined by

[TABLE]

Let $\left(X_{1},\ldots,X_{n}\right)$ a sample made of independent observations with the same distribution as $X$. Denote $\mu$ the distribution of $X$ and $\mu_{n}$ the empirical distribution of the sample $\left(X_{1},\ldots,X_{n}\right)$. The following Lemma can be found in [Lugosi and Nobel, 1996, lemma 1].

Lemma 1

Let $\mathcal{A}$ be any collection of partitions of $\mathbb{R}^{p}$. For every $n\geq 1$ and every $\epsilon>0$,

[TABLE]

4.3 Fréchet regressogram

Let $\mathcal{L}_{n}=\left\{(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\right\}$ be a learning sample made of independent observations with same distribution as $\left(X,Y\right).$ Let $\pi_{n}$ a partitioning rule, we define the Fréchet regressogram estimator by

[TABLE]

where $\pi_{n}[x]$ denotes the unique cell containing $x$. The goal is then to show that under certain assumptions on the metric space $(\mathcal{Y},d)$, on the distribution of $(X,Y)$ and on the partitioning rule, this estimator is consistent for the punctual risk as well as for the global risk.

We recall the definitions of doubling dimension and covering numbers given in [Gottlieb et al., 2016].

Definition 1** (Doubling dimension)**

Let $(\mathcal{Y},d)$ be a metric space, let $\lambda_{\mathcal{Y}}>0$ be the smallest positive integer such that every ball in $\mathcal{Y}$ can be covered by $\lambda_{\mathcal{Y}}$ balls of half its radius. The doubling dimension of $(\mathcal{Y},d)$ is then defined as $\operatorname*{ddim}(\mathcal{Y}):=\log_{2}(\lambda_{\mathcal{Y}})$.

Definition 2** (Covering numbers)**

The $\epsilon$-covering number $\mathcal{N}(\epsilon,\mathcal{Y},d)$ of a metric space $(\mathcal{Y},d)$ is defined as the smallest number of balls of radius $\epsilon$ that suffices to cover $\mathcal{Y}$.

The diameter of a metric space $(\mathcal{Y},d)$, denoted $\operatorname*{diam}(\mathcal{Y})$ is defined by $\operatorname*{diam}(\mathcal{Y})=\underset{y_{1},y_{1}\in\mathcal{Y}}{\sup}d(y_{1},y_{2})$. When both the diameter and doubling dimension of the metric space $(\mathcal{Y},d)$ are finite, according to [Gottlieb et al., 2016], the following lemma allows to bound the $\epsilon-$covering number.

Lemma 2

Let $\left(\mathcal{Y},d\right)$ be a metric space with finite diameter $\operatorname*{diam}(\mathcal{Y})<\infty$ and finite doubling dimension $\operatorname*{ddim}(\mathcal{Y})<\infty$. Then, for every $0<\epsilon\leq\operatorname*{diam}(\mathcal{Y})$

[TABLE]

We now state the main result of our analysis.

Theorem 1

Let $\left(\mathcal{Y},d\right)$ with finite diameter $\operatorname*{diam}(\mathcal{Y})$ and finite doubling dimension $\operatorname*{ddim}\left(\mathcal{Y}\right)$. Let $\pi_{n}$ be a partitioning rule on $\mathbb{R}^{p}$, $\Pi_{n}$ be the family of partitions of $\mathbb{R}^{p}$ obtained from $\pi_{n}$ and $\mathcal{V}_{n}[x]=\mathbb{E}(\operatorname*{Vol}(\pi_{n}[x]))$ be the expected volume of the cell containing $x$. Assume that the following properties hold:

1. P1.

We assume that $\left(X,Y\right)$ has uniformly continuous and bounded density $\rho$ and the marginal $\rho_{X}$ verifies $0<\rho_{\min}\leq\rho_{X}$ 2. P2.

$\frac{\mathcal{C}(\Pi_{n})}{n}\rightarrow 0$** 3. P3.

$\frac{\log(\Delta_{n}^{*}(\Pi_{n}))}{n}\rightarrow 0$** 4. P4.

$\frac{\log\mathcal{V}_{n}[x]}{n}\rightarrow 0$** 5. P5.

$\frac{1}{\mathcal{V}_{n}[x]}=o(\frac{n}{\log n})$** 6. P6.

$\operatorname*{diam}(\pi_{n}[x])\rightarrow 0$* almost surely*

then

[TABLE]

Furthermore,

[TABLE]

Proof: The proof can be found in Appendix A.

4.4 Fréchet purely uniformly random trees

In this (sub)section the input space considered is $\mathcal{X}=[0,1]$. As several theoretical works on regression trees, we consider a simplified version of Fréchet trees. Hence, we study a variant of the purely random trees introduced in [Genuer, 2012], denoted Fréchet purely random tree.

Definition 3** (Fréchet purely uniformly random tree)**

Let $\mathcal{L}_{n}=\{(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\}$ be a learning sample of i.i.d measurements in $[0,1]\times(\mathcal{Y},d)$. Let $k_{n}$ be a positive integer and $U_{1},\ldots,U_{k_{n}}$ be $k_{n}$ i.i.d uniformly drawn random variables on $[0,1]$. Denote $U_{(1)},\ldots,U_{(k_{n})}$ the order statistics, the Fréchet purely random tree predictor $FPURT_{n}$ is given by

[TABLE]

with $U_{(0)}=0$ and $U_{(k_{n}+1)}=1$

Corollary 1

Let $\left(\mathcal{Y},d\right)$ with finite diameter $\operatorname*{diam}(\mathcal{Y})$ and finite doubling dimension $\operatorname*{ddim}\left(\mathcal{Y}\right)$. Let $k_{n}$ be an integer depending on $n$. Assume consider the following assumptions:

1. A1.

We assume that $\left(X,Y\right)$ has uniformly continuous and bounded density $\rho$ and the marginal $\rho_{X}$ verifies $0<\rho_{\min}\leq\rho_{X}$ 2. A2.

$k_{n}\rightarrow\infty$* as $n\rightarrow\infty$ and $k_{n}=o(n/\log n)$*

if the assumptions A1 and A2 hold then the Fréchet purely uniformly random tree estimator is consistent for the global risk i.e

[TABLE]

Proof: Let $\pi_{n}$ be the partitioning rule used to build $FPURT_{n}$ and let $\Pi_{n}$ the family of partitions associated with $\pi_{n}$. The interval $[0,1]$ is partitioned into $k_{n}+1$ intervals, then $\mathcal{C}(\Pi_{n})=k_{n}+1$ which implies

[TABLE]

It is easy to show that $\Delta_{n}^{*}(\Pi_{n})\leq n^{k_{n}}$, then we deduce from $k_{n}=o(n/\log n)$ that

[TABLE]

From [Arlot and Genuer, 2014] (page 34-36) we have that the expected volume (diameter in dimension one) of the interval containing $x$ is:

[TABLE]

Hence, $\mathcal{V}_{n}[x]\leq\frac{2}{k_{n}+1}$, then

[TABLE]

Finally, for $x\in\{0,1\}$

[TABLE]

and for every $0<x<1$

[TABLE]

We demonstrated that the properties P1-P5 of Theorem 14 are verified. We thus conclude that the the $FPURT_{n}$ estimator is consistent point-wise consistent as well as consistent for the global risk.

5 Simulation study

In this section, we study the behavior of Fréchet random forests through three simulation scenarios.

5.1 First scenario, longitudinal data

The first scenario, inspired by our real data applications, deals with the analysis of longitudinal data where inputs and outputs are curves. We simulate $n=100,\ 200,\ 400$ and $1000$ observations of $p=6$ input variables according to the following model for any $i=1,\ldots,n$ and for any $j\in\{1,\ldots,6\}$:

[TABLE]

where $X^{(j)}_{i}(t)$ is the observation of the $j$th input variable at time $t$ for the $i$th curve (individual/patient); $t$ browses a regular subdivision of $[0,1]$ with a step size of 0.05, $G^{j}_{i}$ and $G^{\prime j}_{i}\sim B(0.5)$, $\beta_{i}$ and $\beta_{i}^{\prime}\sim\mathcal{N}\left(1,0.3\right)$, $W_{i}^{1}(t)$ is a Gaussian white noise with standard deviation $0.02$ and $f_{j,1}$ and $f_{j,2}$ are defined as follows:

[TABLE]

The terms $G^{j}_{i}$ and $G^{\prime j}_{i}$ allow to randomly affect typical temporal behaviors, defined by $f_{j,1}$ and $f_{j,2}$ functions, to observations. The $\beta_{i}$ and $\beta^{\prime}_{i}$ are dilatation/shrinkage terms of $f_{j,1}$ or $f_{j,2}$, while $W_{i}^{1}(t)$ corresponds to an additive noise. As illustrated in Figure 1, for each input variable, the observed trajectories are variations of the typical temporal behavior functions. The observations are divided into two groups of trajectories.

Output variable $Y$ is simulated in a similar way. The pair $(G_{i}^{1},G_{i}^{2})$ is used to determine a trajectory for the output variable, this is the primary link between $X_{i}$ and $Y_{i}$

[TABLE]

where $Y_{i}(t)$ is the $i$th output curve measured at time $t$; $t$ browses the same subdivision as in (18), $\beta_{i}$ are the same coefficients used in (18), $W_{i}^{2}\left(t\right)$ is a Gaussian white noise with standard deviation $0.05$ and $g_{j,k}$ are given by:

[TABLE]

The response curves are distributed according to four different trajectory shapes, one for each pair of possible trajectory shapes for the first two input curve variables $X^{(1)}$ and $X^{(2)}$. Of note, the variables $X^{(3)}$ and $X^{(4)}$ are simulated using the same temporal functions as variables $X^{(1)}$ and $X^{(2)}$; however, the trajectories of variables $X^{(3)}$ and $X^{(4)}$ are simulated from $G^{\prime j}$ and not from $G^{j}$ and thus have no relation with the output variable $Y$.

5.2 Second scenario, predict curves with images, scalars and curves

In this scenario, we want to predict output curves from inputs that are curves, scalars and images to illustrate the flexibility of the Fréchet RF method, in particular its ability to learn about different types of inputs and outputs. The input curve variables are simulated according to the model (18) of the first scenario with $\beta_{i}$ and $\beta_{i}^{\prime}$ drawn according to $\mathcal{N}(1,1)$ in order to have large variations of the curves around their average temporal behavior. Similarly, the output curves are simulated according to the model (19) of the first scenario. Let $\left(\mathcal{M}^{1}_{i}\right)_{i}$ and $\left(\mathcal{M}^{2}_{i}\right)_{i}$ two sequences of handwritten images of numbers 1 (for $\mathcal{M}^{1}_{i}$) and 2 (for $\mathcal{M}^{2}_{i}$) randomly drawn from the MNIST dataset [LeCun and Cortes, 2010]. We simulate two input image variables $I^{(1)}$ and $I^{(2)}$ according to the following model:

[TABLE]

where $G^{j}_{i}$ are the same draws as those used to simulate the input and output curves in model (18) and model (19). Finally, consider the two real input variables $R^{(1)}_{i}=\beta_{i}$ and $R^{(2)}_{i}=\beta_{i}^{\prime}$, where the $\beta_{i}$ and $\beta_{i}^{\prime}$ are the same as those used to simulate the input and output curves. The first variable $R^{(1)}$ determines the intensity of the contraction/expansion of the $X^{(1)}$ and $X^{(2)}$ response curves. It is important to note that the link between the output curves and the input variables is entirely contained in the pairs $(G^{1}_{i},G_{i}^{2})$ which determine the general shape of the output curve as well as the $\beta_{i}$ which determine the compression/expansion of the output curves. The pairs $(G^{1}_{i},G_{i}^{2})$ as well as the $\beta_{i}$ are used to simulate the first two curve input variables $X^{(1)}$ and $X^{(2)}$. However the two image variables are constructed only from the pairs $(G^{1}_{i},G_{i}^{2})$ and the scalar variables are $\beta_{i}$ and $\beta_{i}^{\prime}$.

5.3 Third scenario, predict images with curves, a toy example

The purpose of this scenario is to illustrate the ability of the Fréchet RF method to predict images from input curves. We simulate a dataset of $n=500$ observations, the input curve variables are simulated according to the model (18) of the first scenario. As in the second scenario, the output images are taken from the MNIST dataset [LeCun and Cortes, 2010]. For any $k=1,\ldots,8$ and for any $i=1,\ldots,n$ we note $\mathcal{M}_{i}^{k}$ the random draw of the handwritten $k$ digit in the MNIST dataset for the $i$th observation. Let the pair $(G_{i}^{1},G^{2}_{i})$ used to attribute their shape to the curves of the first two input variables for the $i$th observation and $\beta_{i}$ the expansion/contraction parameter of these same curves, then the output images are drawn according to the combinations summarized in the Table 1.

As in the first scenario, the output images depend only on the first two input variables, the link between the images and the curves is entirely contained in the pairs $(G_{i}^{1},G_{i}^{2})$ as well as in the $\beta_{i}$. This means that the handwritten number images then depend both on the shape of the curves of the first two input variables and their amplitude.

5.4 Results

5.4.1 First scenario

First, we need to determine a metric for each input and output space. In the case of longitudinal data i.e. when repeated measurements of quantitative variables are available over time. The observations for $p$ input and one output variables can thus be represented by time-dependent curves. In this case, the $i$-th observation $X_{i}$ is a curve from $\mathcal{I}_{1}\times\cdots\times\mathcal{I}_{p}\subset\mathbb{R}^{p}_{+}$ to $\mathbb{R}^{p}$ (where $\mathcal{I}_{1}=[0,1]$ and $\mathcal{I}_{2}=[0,1]$ in the first scenario), and $Y_{i}$ is a curve from $\mathcal{J}\subset\mathbb{R}_{+}$ to $\mathbb{R}$. We choose to equip the resulting curves spaces with the Fréchet distance $d_{\mathcal{F}}$ introduced in [Fréchet, 1906] defined for two real-valued curves $f$ and $g$ with support in $\mathcal{I}\subset\mathbb{R}_{+}$ as

[TABLE]

where $\alpha$ and $\beta$ are any re-parameterizations of $\mathcal{I}$. The definition is the same in the discrete case (polygonal curves), except that $t$ takes values on $\mathcal{I}$ by intervals, see [Alt and Godeau, 1995] for a full description of Fréchet distance for discretely sampled curves. This distance is a natural measure of similarity between the shapes of curves and has been widely used in various applications such as signature authentication [Zheng et al., 2008], path classification [Genolini et al., 2016] and speech recognition [Kwong et al., 1998]. Note that, unlike several classical distances, the calculation of the Fréchet distance does not require the same number of measurements, nor the same observations times on the two trajectories. Once we have determined the metrics used for the different spaces we need to define the split function used to cut on the input spaces. Finally, the 2-means algorithm for longitudinal data using Fréchet distance and Fréchet mean introduced in [Genolini et al., 2016] is chosen on each input space to determine the different competing splits. This split function called kmlShape is an adaptation of the $k$-means method tailored to one-dimensional curves. It allows to find groups of trajectories based on their shapes (which are usually not found by conventional methods, e.g. based on Euclidean distance).

Fréchet trees and Fréchet random forests were compared on simulated datasets to standard CART trees [Breiman et al., 1984] and standard random forests [Breiman, 2001] as well as standard existing methods for longitudinal data analysis such as linear mixed effects model (LMEM) with a random intercept and a random effect on time and the boosting functional regression method FDboost [Brockhaus et al., 2017] with optimized number of iterations.

The prediction errors (mean squared error) of all the methods are estimated on several sample sizes $n=100,\ 200,\ 400$ and 1000 using for each sample size, $100$ datasets simulated according to models (18) and (19). For each simulated dataset $\mathcal{L}_{n}$, we randomly divide $\mathcal{L}_{n}$ into a training set (with $0.8n$ observations) and a test set (made of the remaining $0.2n$ observations). The Fréchet distance is used on the curved input and output spaces to build Fréchet trees and Fréchet random forests, however in order not to advantage our method, prediction errors are calculated with the usual $L^{2}$ Euclidean distance (time by time) which benefits to the standard approaches like CART trees, RF and FDboost.

The number of randomly drawn variables mtry at each node has usually a strong impact on random forests performance: if mtry is too small, individual trees would give too poor predictions, and if mtry is too high, the collection of trees could be not diverse enough ([Díaz-Uriarte and Alvarez De Andres, 2006];[Genuer et al., 2008]). As illustrated in Figure 2 the prediction error (MSE) of the Fréchet random forest decreases as the value of the mtry increases. In all our experiments in the first scenario, we chose mtry=5. We set the number of trees $q$ to $250$ (justified by the fact that, in this experiment, the OOB error stabilizes as soon as $100$ trees are included in the forest). The standard random forest was composed of 500 trees and the mtry parameter was otpimized to 2. The number of iterations for FDboost is selected between 1 and 500 through the internal procedure of the package.

As illustrated in Figure 3 for any sample size, FDboost, Fréchet tree, and Fréchet random forests clearly outperform the standard LMEM, CART and RF methods. Not surprisingly, the transition from a Fréchet tree to a Fréchet RF greatly improves predictive capacity by reducing prediction error and error variance. For instance, when $n=100$, the estimated MSE obtained with a Fréchet tree is 0.047 while the one obtained with a Fréchet RF is 0.028 which is a 40% decrease in prediction error. Even though FDboost (our principal competitor) shows very good performances, Fréchet tree and Fréchet RF are the methods that obtain the lowest prediction errors for all sample sizes. More precisely, for small dataset ($n=100$) FDboost obtains an estimated MSE of 0.05 while Fréchet tree and Fréchet RF obtain respectively 0.047 and 0.028 while for large dataset ($n=1000$) the estimated MSE of FDboost is 0.031 and the Fréchet tree and Fréchet RF estimated MSE are respectively 0.012 and 0.006. Finally, note that the prediction error of the FDboost, Fréchet tree and Fréchet RF methods decreases as the sample size $n$ increases which is not the case with other methods that keep a stable prediction error. Additionally, this decrease is much larger with the Fréchet tree and Fréchet RF methods than with the FDboost method. Moreover, the error prediction of the Fréchet RF converges to zero as $n$ tends to infinity.

The extremely randomized version of Fréchet random forests introduced in Section 3.3 has some advantages over the Fréchet RF method. In particular, they are easy to implement, can be used for any type of data and reduce calculation times. In order to verify this claim we calculate the prediction error obtained by extremely randomized Fréchet forests (ERFRF) for different values of ntry on 100 data sets of size $n=100$ simulated according to the first scenario. As shown in Figure 4, the prediction error of the ERFRF method decreases as the value of the ntry increases. When ntry is large enough (here ntry=3), the error obtained by ERFRF is similar to that obtained by Fréchet RF. Moreover, the execution time of an ERFRF is much lower than that of a Fréchet RF. For example, the build time of an Fréchet RF is 281 seconds while the build time of an ERFRF with ntry=3 is 191 seconds which is 30% lower. Similar results are obtained on larger datasets (not shown here).

As mentioned in the presentation of the Fréchet distance at the beginning of this section, using the Fréchet distance allows to calculate the distance between two curves measured at different times. Thus, having missing observation times for some curves does not prevent the construction of the trees, as long as not all observation times are missing for a given curve. In order to study the robustness of Fréchet RF to missing observations, we simulate new datasets with $n=100$ individuals according to models (18) and (19) by randomly removing 10%, 20% and 30% of the observation times for each curve. It is important to note that the removed observation times are different for each curve. For example, the observations removed for the first variable of the first individual will not necessarily be the same as those removed for the second or third variable or even the output curve of the same individual. It is then impossible to use the standard LMEM, RF and FDboost methods (it is always possible to use the CART method by removing the missing observations for the output curves). As shown in Figure 5, the prediction error obtained by Fréchet RF increases slightly as the percentage of missing observations increases. Moreover, the prediction error obtained with Fréchet RF on simulated data sets with 30% missing data remains competitive with that obtained by the FDboost method on datasets without missing observations.

It is rather common in biostatistical applications and more particularly in clinical trials to have a response variable observed after the measurement times of the input variables, this is the case of our application on DALIA-I vaccine trial (see Section 6). In order to study the stability of Fréchet RF method to time shifts, we transform the output curves by shifting them: i) by the same time shift of 1 for all the curves, i.e., the output curves are observed on windows $[1,2]$ instead of $[0,1]$ (keeping the same shapes); ii) by randomly shifting each of them according to a uniform $\mathcal{U}([0,1])$, making the windows of observation of the output curves all different in this case (see Figure 6 for the simulated dynamics according to the time shifts). As illustrated in Figure 7, the constant time shift for the response curves has no influence on the Fréchet RF prediction error. When the offsets are randomly drawn for each output curve, the prediction error increases slightly to an average error of 0.039, however note that this error is still 20% smaller than the average error of 0.05 obtained by FDBoost on simulated data without time shifts.

Finally, Figure 8 gives the importance scores of variables calculated with the Fréchet RF method on 4 datasets of size $n=100$ simulated according to models (18) and (19):

With no time shifts on the output curves and no missing observation times. 2. 2.

With random time shifts according to a uniform $\mathcal{U}([0,1])$ on the output curves but with no missing measurement times. 3. 3.

With no time shifts but with 30% missing observation times. 4. 4.

With 30% missing observations and time shifts on the output curves.

This graph shows that neither time shifts nor missing observation times have an impact on the importance of the variables. In fact, the first two variables (those related to the output variable) are always the ones with the highest importance scores. The other four variables (unrelated to the output variable) have extremely low importance scores compared to the first two variables.

As a conclusion, we illustrate the superiority on longitudinal data (in terms of prediction error) of the Fréchet trees and Fréchet RF methods compared to the standard LMEM, CART, RF methods as well as the longitudinal boosting method FDboost. In addition, we illustrate the great robustness of the method to missing data and time shifts, both in terms of prediction error and the importance of the variables. Lastly, we show that the extremely randomized variant ERFRF can obtain a prediction error similar to that of Fréchet RF while having lower computation times, making it a method of choice for analyzing very large datasets.

5.4.2 Second scenario

The Fréchet distance is used on curve spaces while the standard Euclidean distance is used on scalar spaces and image variables. Since there is no comparison with other methods in this scenario, the OOB error will be used as a measure of the performance of the Fréchet RF. Throughout this section we study the ERFRF method, the version implemented in our package FrechForest that can handle images and shapes as inputs.

We study the OOB error obtained by ERFRF according to the types of input variables (images, curves or scalars) on 100 datasets of size $n=$100 simulated according to the second scenario. We consider the following models:

Only scalar variables $R^{(1)}$ and $R^{(2)}$ are used to predict output curves. 2. 2.

Only curve variables $X^{(1)},\ldots,X^{(6)}$ are used to predict output curves. 3. 3.

Image variables $I^{(1)}$ and $I^{(2)}$ and scalar variables are used. 4. 4.

all variables *i.e. *curves, scalars and images are used to predict output curves.

Note that case 2 corresponds to the first simulation scenario. Figure 9 shows an example of an extremly randomised Fréchet tree of depth 2 (only the first three splits are shown here) for each model above. When the models incorporate different types of inputs, in the case of models 3) and 4), the constructed trees are mixed in the sense that they can alternate the split spaces. For example, in the case of model 4), Figure 9 shows an example of a tree with the first three splits in the three different types of input spaces: curves, scalars and images. The parameters mtry and ntry are selected for each model to minimize the OOB error of the ERFRF, and are therefore different for each model.

As shown in Figure 10, the highest OOB error is obtained when only scalar variables are used. When the image variables are added to the scalar variables, the OOB error is the same as the one obtained on the model using only the input curves. This was expected since the input curve variables provide the same information as the image and scalar variables combined. More precisely, the input curve variables provide both information on the shape of the output curves as well as on their amplitude, whereas the information on the shape is only provided by the images and the ones on the amplitude is only provided by the scalars. Individually, the input variables of images or scalars provide only part of the information that is provided by the input variables of curves. Finally, when the image and scalar variables are added to the curve variables, the OOB error of the ERFRF decreases. This is explained by the fact that in some cases, when the contraction or dilation of the input curves is too large, the dilated or contracted curves may have a very different shape than their initial shape and thus lose the information they brought due to their shape. Thus the addition of image variables allows to always have access to information on the shapes of the output curves. Finally, the results of these simulations emphasizes the main strength of the ERFRF method, which is to handle heterogeneous data, i.e. input and output variables of different natures.

5.4.3 Third scenario

The Fréchet distance is used on the curve spaces, i.e. on the 6 input variables. The distance used on the output space is the standard Euclidean distance. A Fréchet RF is constructed with $q=500$ trees (justified by the fact that the OOB error of the Fréchet RF becomes stable as long as 350 trees compose the forest). Similarly, the mtry parameter is set to 5. As shown in Figure 11, OOB predictions of output images always give the correct written digit. However, ghosting can be seen on some digit predictions. This is due to the simulation scheme itself. The input curves only give information about the written number, and do not provide any information about its individual characteristics such as the width of the number, its height, the presence or not of a loop (for writing a 2 for example). More precisely, there is within the same group of numbers (for example the set of numbers 4 drawn) a variability in the written numbers that is not explained by the input curves. By introducing variables that provide information on the fine characteristics of each written number (such as its height, width, etc.) we would get even more accurate predictions. Moreover, it is noticeable that this phenomenon of ghosting is not present for numbers that have a very low variability in their writing such as the number 1. In order to highlight this point a Fréchet RF is constructed on the same simulated dataset and the images are replaced by factors indicating what the written number is. When the outputs are images, the percentage of explained variance is 20%, which was expected since there is a large variability between the same numbers that is not explained by the input curves. When the outputs are factors expressing the written numbers, the percentage of explained variance is 98%. It is therefore clear that the link between the output images and the input curves relates only to the digits and not to its individual characteristics. So even if the explained variance percentage is only 20% for the images, the Fréchet RF (almost) always predicts the right digit.

6 Application to real data

6.1 DALIA vaccine trial

DALIA is a therapeutic vaccine trial including 17 HIV-infected patients who received an HIV vaccine candidate before stopping their antiretroviral treatment. For a full description of the DALIA vaccine trial we refer to [Lévy et al., 2014]. At each harvest time before stopping their treatment, 5,399 gene transcripts which significantly vary over time during the vaccination phase were selected by [Hejblum et al., 2015] among more than 32000. The plasma HIV viral load (which was log-transformed) for every patient was measured at each harvest time after the antiretroviral treatment interruption (called HAART interruption). In this application the measurement times of the inputs (gene transcripts) differ from the ones of the output (HIV viral load). The objective is to be able to predict the HIV viral load dynamics after antiretroviral treatment interruption for a patient given the evolution of his/her gene expression during the vaccination phase [Thiébaut et al., 2019]. Figure 12 illustrates the design of the DALIA vaccine trial and the dynamics of the viral replication after antiretroviral treatment interruption with a large between-individuals variability. The analysis with Fréchet random forest was performed on the 17 patients. The $mtry$ parameter was fixed to $1500$ and the number of trees, $q$, was set to $500$. The OOB error of the Fréchet random forest converged and stabilized for almost $100$ trees composing the forest. Figure 13 illustrates both the OOB predictions and the predictions on the learning samples (fits) of the evolution of the viral load after the HAART interruption for 4 patients of the vaccine trial.

The predictions of the Fréchet forest on the learning sample were close to the observed viral load curves. Moreover, despite a very small number of individuals, the OOB predictions obtained with this forest are quite close in shape to the true curves.

Among the 100 most important variables, many belong to the groups of genes (modules) that were selected in [Thiébaut et al., 2019] because i) their dynamics was influenced by the vaccine ii) their abundance after vaccination were associated with the maximum of the observed viral load. For instance, five genes from the inflammation module 3.2 and three genes from the T cell module 4.1 were selected with the current approach. Both groups of genes were extremely relevant in the context of such vaccine that generates strong T cell response.

Thus, the Fréchet random forests method applied on the complex example of the DALIA vaccine trial is extremely effective both for its capacity to predict the output variable as well as for its ability to find relevant genes in order to explain the evolution of the viral load after the treatment interruption. Previous analyses performed in [Thiébaut et al., 2019] were done by looking at the association between one time of measurement for the input gene expressions and only one characteristic of the viral load dynamics: its maximum value. The Fréchet random forest allows a direct analysis of the whole longitudinal information available. It should be noted that standard CART trees and random forests methods cannot be used on such an application. Indeed, both the number and the observation times of the input and output variables were different.

6.2 LIGHT vaccine trial

LIGHT is a therapeutic vaccine trial including 97 HIV-infected patients. Using the data available in this trial, the objective of the present analysis was to assess the capacity of predicting the abundance of CD4 T cells using gene expression data as measured by RNA sequencing in whole blood, dealing with high-dimensional longitudinal data with extremely unbalanced trajectories. The dataset is composed by 1150 input variables of genes abundance. Those 1150 input genes were preselected using dearseq [Gauthier et al., 2019], a method for differential expression analysis. These are found to be differentially expressed overtime among more than 17000 genes. A targeted control rate for the False Discovery Rate was used at a nominal level of 5%.

Furthermore, there were 234 observations in the dataset: one to four measurements were available for each patient. In Section 5.4.1 we showed that the Fréchet RF method was very robust to unbalanced experimental designs. However, in the LIGHT dataset there is a significant proportion of patients with only one observation time, which is an extreme case that can cause some problems in predicting trajectories with three or four measurement times. We then consider the output as scalars, i.e. all the output observations are independent. We used 4 different RF models:

A standard RF. 2. 2.

A Fréchet RF considering the inputs as scalars. 3. 3.

A Fréchet RF considering the inputs as curves. 4. 4.

A Fréchet RF considering both the inputs as curves and as scalars.

As in our simulations, the Fréchet distance is used for curve spaces while the Euclidean distance is used for scalar spaces. For each random forest model the mtry parameters as well as the number of trees composing the forest are optimized to minimize the OOB error. The prediction error estimated on OOB samples for the different RF models are compared. The standard RF method as well as the Fréchet RF method with inputs and outputs treated as scalars obtain the highest OOB errors with 22.65 and 22.61 respectively. When the input is treated as curves, the Fréchet RF method obtains an OOB error of 17.59 which represents a decrease of more than 20% of the OOB error compared to the first two models. Unsurprisingly, and in line with our simulations, considering repeated observations from the same patient as curves and partitioning the input space according to the shape of these curves greatly improves predictive performance. When the inputs are treated as both curves and scalars, the prediction error obtained from the Fréchet RF is 17.36. Thus, adding the same input variables considered as scalars does not improve the OOB error already obtained by the Fréchet RF on the curve variables only. Genes as scalars (i.e. when all observations are independent) do not therefore provide additional information compared to genes as curves, because integrating these genes in the form of scalars does not improve the prediction error.

Finally, the importance scores of the variables are considered. Figure 14 shows the importance scores of the 200 most important input genes calculated with the RF Fréchet built on the genes considered as curves as well as scalars. The genes as curves globally obtain much higher importance scores than those obtained by the genes as scalars. A closer analysis of the importance scores of genes according to their type (curve or scalar) shows that the importance score for a gene as curve is systematically higher than that obtained by the same gene as a scalar. Moreover, of the 200 genes in curve form that obtain the highest importance scores, less than 47% are common to the 200 most important genes for genes in scalar form. To go further, only 42% of the 200 genes in the form of curves that obtain the highest importance scores are common to the 200 most important genes found with the standard RF method. Thus, in the LIGHT vaccine trial, considering the observations from the same patient as a trajectory not only greatly improves the prediction error but also allows us to find genes that were not previously found with the standard RF method. From the importance scores of the variables we selected the 100 most important genes for the genes considered as curves. Among these 100 genes, many of them (e.g. CD8A, CD40, TLR3, TRBV1, TRBV18…) were completely related to T cell pathways as expected because they are associated to the abundance of CD4+ T cells. Hence, Th1 pathway, Th2 pathway, T cell exhaustion pathways were among the most enriched pathways among the 100 selected genes (Ingenuity Pathway Analysis). Interestingly, the genes and pathways associated to the analysis considering inputs as scalar were less relevant with for instance the absence of the genes associated to communication between immune cells (CCL4, CCL5, CCL7) that were selected when taking into account the curves. Finally, the RF Fréchet considering the repeated measures as curves provided very relevant gene selection in this application.

7 Discussion

Two new tree-based methods, Fréchet trees and Fréchet random forests, for general metric spaces-valued data were introduced. Let us emphasize that the proposed methods are very general. Indeed, input variables can thus all be of different kinds, each one having its own metric, and the kind of the output variable can also be a different one.

The example of learning curve shapes in the context of longitudinal data was presented to illustrate the capacity of the methods to learn from data in unordered metric spaces. A simulation study in this framework demonstrated the superiority of Fréchet trees and forests over the existing classical methods, both in terms of prediction error as well as robustness and flexibility. An important aspect highlighted in our study is the great robustness of Fréchet trees and Fréchet random forests. Indeed, our simulations illustrated the ability to handle missing data as well as different observation times for the different variables, which is common in longitudinal datasets. Two other simulation scenarios demonstrated the capacity of the methods to simultaneously handle data of different natures such as curves, images, scalars, factors, shapes etc. This great flexibility allows the construction of more efficient predictors while being able to compare the information provided by each of these variables of different natures thanks to the importance score. In this article, we analyzed two high-dimensional longitudinal datasets from vaccine trials. For the first time and due to the very high flexibility of Fréchet RF, it was possible to associate the entire evolution of the transcriptome during the vaccination phase with the entire evolution of the immune response after interruption of treatment on DALIA-I dataset. We have thus highlighted the groups of genes that best explain the different immune response after interruption of treatment. Finally, within the framework of LIGHT vaccine trial on which the classical methods stumbled, we highlighted the superiority of the Fréchet RF method. We illustrated that regression on curve shapes could greatly improve the prediction error while selecting new variables.

However, there are two main limitations to Fréchet trees and forests: the first one is that the Fréchet mean has to exist in the output space [Le Gouic and Loubes, 2017] and has to be fairly approximated. The second concerns the computation time. Indeed, even if the proposed approaches have been fully coded in the R package FrechForest for the trajectories case, Fréchet random forests can still be computationally intensive. This problem can be alleviated by the fact that, as all forests methods, they are easily parallelized (the different trees can be built in parallel).

Finally, we are working on the rates of convergence of Fréchet trees and Fréchet random forests. In parallel, we are working on the efficient implementation of metrics adapted to the image problem, such as the Wassertein distance [Vallender, 1974], in order to apply the Fréchet RF method to large brain imaging databases.

Appendix A Proof of Theorem 14

First, we demonstrate the point-wise consistency given by (13). We introduce the following quantity

[TABLE]

From 11 we have $T_{n}(x)=\underset{y\in\mathcal{Y}}{\operatorname*{arg\,min}}\ r_{n}\left(x,y\right)$. First, we use the following classical upper bound in $M$-estimation:

[TABLE]

We are going to decompose the above supremum in several terms that we are going to appropriately upperbound to obtain their decay to zero under the assumptions Theorem 14. Consider a $\delta$ covering of $\mathcal{Y}$ with centers $\{y_{\alpha}\}_{\alpha=1}^{Q}$ where $Q=N(\delta,\mathcal{Y},d)$. Thus, for every $y\in\mathcal{Y}$, there is $\alpha=\alpha_{y}\in\{1,\ldots,Q\}$ such as $d\left(y,y_{\alpha}\right)<\delta$. We introduce the following quantity

[TABLE]

Then, the following decomposition is used

[TABLE]

We are now going to derive upper bounds for each of the four terms above that do not depend on $y$. Let us start with the term $(i)$ of (25), we introduce the following event

[TABLE]

We can upper bound the probability of the complementary of the event $\mathcal{E}_{n}$ (denoted $\mathcal{E}_{n}^{c}$) as

[TABLE]

Then, we upper bound the last probability using Lemma 1

[TABLE]

On the event $\mathcal{E}_{n}$, one has that $\left|\frac{\frac{1}{n}\sum_{i=1}^{n}\mathbb{1}\{X_{i}\in\pi_{n}[x]\}}{\mathbb{P}\left(X\in\pi_{n}[x]\right)}\right|<\frac{3}{2}$ which implies that

[TABLE]

On the complementary of the event $\mathcal{E}_{n}$, we use similar arguments and the upper bound $\left|\frac{\frac{1}{n}\sum_{i=1}^{n}\mathbb{1}\{X_{i}\in\pi_{n}[x]\}}{\mathbb{P}\left(X\in\pi_{n}[x]\right)}\right|\leq\frac{1}{\mathbb{P}\left(X\in\pi_{n}[x]\right)}$ to derive that

[TABLE]

Therefore, we finally obtain that

[TABLE]

Now we consider the term $(ii)$ in (25). To this end, we propose to bound the following probability

[TABLE]

For a fixed $\alpha\in\{1,\ldots,Q\}$, we define $W_{i}=\frac{d^{2}(Y_{i},y_{\alpha})\mathbb{1}\{X_{i}\in\pi_{n}[x]\}}{\mathbb{P}\left(X\in\pi_{n}[x]\right)}$, and we thus have that

[TABLE]

which can be controlled thanks to Bernstein’s inequality by finding upper bounds on $|W_{i}|$ and $\operatorname*{var}(W_{i})$. To this end, we first derive a lower bound on $\mathbb{P}\left(X\in\pi_{n}[x]\right)$.

[TABLE]

Therefore, we obtain that

[TABLE]

Moreover,

[TABLE]

Then, by Bernstein’s inequality, we have for every $\alpha\in\{1,\ldots,Q\}$

[TABLE]

For $\epsilon<1$ we have

[TABLE]

We deduce from Equation (A) and Lemma 2

[TABLE]

Let us now bound the term $(iii)$ in (25) as follows

[TABLE]

We combine inequalities (29), (32) and (33) such that with probability, we have: $1-2\exp\left(\operatorname*{ddim}(\mathcal{Y})\log\frac{2\operatorname*{diam}(\mathcal{Y})}{\delta}-Cn\epsilon^{2}\mathcal{V}_{n}[x]\right)$.

[TABLE]

Thanks to the assumptions $\frac{\mathcal{C}(\Pi_{n})}{n}\rightarrow 0$, $\frac{\log(\Delta_{n}^{*}(\Pi_{n}))}{n}\rightarrow 0$ and $\frac{\log\mathcal{V}_{n}[x]}{n}\rightarrow 0$ of Theorem 14, the term $\frac{\Delta_{n}^{*}(\Pi_{n})2^{\mathcal{C}(\Pi_{n})}\exp-n/128}{\rho_{\min}\mathcal{V}_{n}[x]}$ appearing in the right hand side of the Inequality (36) converges to zero. Hence, there is a constant $D$ such that

[TABLE]

for every $n$. Thus we deduce the following inequality that holds with probability $1-2\exp\left(\operatorname*{ddim}(\mathcal{Y})\log\frac{2\operatorname*{diam}(\mathcal{Y})}{\delta}-Cn\epsilon^{2}\mathcal{V}_{n}[x]\right)$

[TABLE]

with $B=\operatorname*{diam}(\mathcal{Y})(5+8D)$.

Let $s>0$ and $\delta=n^{-s}$, for $s$ large enough $B\delta$ is bounded by $\epsilon$. Thus, for $s$ large enough we deduce that

[TABLE]

Under the assumption on $\frac{1}{\mathcal{V}_{n}[x]}=o\left(\frac{n}{\log n}\right)$, the probability upper bound on the right hand side in Inequality (38) becomes summable over $n$. We thus conclude the almost sure convergence of $\underset{y\in\mathcal{Y}}{\sup}\left|r_{n}(x,y)-r^{E}(x,y)\right|$ towards zero by the Borel-Cantelli Lemma.

Finally, we analyze the term $(iv)$ in (25) $\left|r\left(x,y\right)-r^{E}\left(x,y\right)\right|$ For fixed $x_{0}\in\mathbb{R}^{p}$ and $y_{0}\in\mathcal{Y}$, we have

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

Since $\rho$ is uniformly continuous, for every $(x_{0},y)\in\mathbb{R}^{p}\times\mathcal{Y}$, $\forall\epsilon>0,\exists\delta^{1}_{\epsilon}>0$ such that $||x_{0}-x||\leq\delta^{1}_{\epsilon}\Rightarrow|\rho(x_{0},y)-\rho(x,y)|\leq\epsilon$. Thus, there exists $\delta^{1}_{\epsilon}>0$ such that

[TABLE]

Using the same argument of continuity on the density $\rho_{X}$, for all $\epsilon$, there is $\delta_{\epsilon}^{2}$ such that

[TABLE]

We define $\delta_{\epsilon}=\min(\delta_{\epsilon}^{1},\delta_{\epsilon}^{2})$. We will apply the dominated convergence theorem to conclude. To this end, we remark that for every sequence of functions $(f_{n})_{n}$, $(g_{n})_{n}$ and for every functions $f$ and $g$ we have

[TABLE]

We take

[TABLE]

[TABLE]

We deduce the following upper bound

[TABLE]

with the same arguments we also get using (43)

[TABLE]

From the assumptions of Theorem 14, we have that $\operatorname*{diam}(\pi_{n}[x_{0}])$ converges towards zero almost surely. Hence, with probability 1, for every $\delta_{\epsilon}$, there is $N_{\epsilon}>0$ such that for every $n\geq N_{\epsilon}$, $\operatorname*{diam}(\pi_{n}[x_{0}])\leq\delta_{\epsilon}/2$. Thus, for every $n\geq N_{\epsilon}$, $\operatorname*{Vol}(\pi_{n}(x_{0})\backslash B(x_{0},\delta_{\epsilon}))=0$ almost surely. Then from (A) and (45) we deduce that for every $n\geq N_{\epsilon}$ the following inequalities hold almost surely

[TABLE]

Finally, we deduce from (46)

[TABLE]

Moreover

[TABLE]

Using the dominated convergence theorem we thus get

[TABLE]

Finally, we demonstrate the weak consistency given by (14). The proof uses the arguments from [Hein, 2009]. Under the assumptions of Theorem 14, for every $x\in\mathbb{R}^{p}$, one has that, $\lim\limits_{n\rightarrow\infty}r\left(x,T_{n}(x)\right)=r\left(x,\phi^{*}(x)\right)$ almost surely. Now, remark that

[TABLE]

As $\operatorname*{diam}\left(\mathcal{Y}\right)<\infty$, we have that $\mathbb{E}\left(r\left(X,T_{n}(X)\right)\right)<+\infty$ and $\mathbb{E}\left(r\left(X,\phi^{*}(X)\right)\right)<+\infty$. Therefore, an extension of the dominated convergence theorem given in [Glick, 1974] allows to conclude.