<figure><img src="image_1.png" /></figure>

HAL
open science

# Quantile oriented sensitivity analysis with random forest based on pinball loss

Ri Wang, Véronique Maume-Deschamps, Clémentine Prieur

► To cite this version:

Ri Wang, Véronique Maume-Deschamps, Clémentine Prieur. Quantile oriented sensitivity analysis with random forest based on pinball loss. 2024. hal-04606380

HAL Id: hal-04606380

https://inria.hal.science/hal-04606380

Preprint submitted on 10 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantile oriented sensitivity analysis with random forest based on pinball loss

Ri Wang, Véronique Maume-Deschamps, Clémentine Prieur

► To cite this version:

Ri Wang, Véronique Maume-Deschamps, Clémentine Prieur. Quantile oriented sensitivity analysis with random forest based on pinball loss. 2024. hal-04606380v2

HAL Id: hal-04606380

https://inria.hal.science/hal-04606380v2

Preprint submitted on 30 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantile oriented sensitivity analysis with random
forest based on pinball loss

Ri Wang<sup>(1)</sup>, Véronique Maume-Deschamps<sup>(1)</sup> and Clémentine Prieur<sup>(2)</sup>

<sup>(1)</sup>Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon,
Université Jean Monnet, ICJ UMR5208, 69622 Villeurbanne, France.

<sup>(2)</sup>Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France

Abstract

Global Sensitivity Analysis (GSA) is an important tool to better
understand the behavior of black box models. Among the numerous
methods for GSA, variance-based approaches have received much at-
tention. Only a few papers focus on Quantile Oriented Sensitivity
Analysis (QOSA), which can help in analysing the behavior of the
response at different quantile levels. Moreover, existing QOSA estima-
tion methods have flaws: bias when input variables are dependent, loss
of accuracy and efficiency as input space dimension increases. In this
paper, we propose a new estimation procedure of QOSA indices based
on the notion of projected random forest, with the initial random for-
est built from a criterion designed for quantiles: the pinball loss also
known as quantile loss.

**Key words:** Quantile Oriented Sensitivity Analysis, random forest, conditional distribution, conditional quantile, projected forest.

# 1 Introduction

Sensitivity Analysis (SA), as defined in Saltelli et al. (2004), is the study of
how the uncertainty in the output of a system can be divided and allocated
to different sources of uncertainty in its inputs. It is an invaluable tool to
understand the behavior of numerical models, determining the most con-
tributing input variables and ascertaining interaction effects within model
(see, e.g., Da Veiga et al. (2021) for a recent review). Global Sensitivity
Analysis (GSA) focuses on the impact of inputs on some output of interest,
when varied over their whole domain (Iooss and Lemaître (2015)). Variance-
based methods are well-established and widely used for GSA. In this context,
Sobol' indices introduced in Sobol' (1993) are very popular.

In the present paper, we focus on Quantile Oriented Sensitivity Analysis
(QOSA), a special case of indices based on contrast functions considered in
Fort et al. (2016). QOSA indices are designed to capture the impact of inputs
# Quantile oriented sensitivity analysis with random forest based on pinball loss

Ri Wang<sup>(1)</sup>, Véronique Maume-Deschamps<sup>(1)</sup> and Clémentine Prieur<sup>(2)</sup>

<sup>(1)</sup>Université Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, 69622 Villeurbanne, France.

<sup>(2)</sup>Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France

## Abstract

Global Sensitivity Analysis (GSA) is an important tool to better understand the behavior of black box models. Among the numerous methods for GSA, variance-based approaches have received much attention. Only a few papers focus on Quantile Oriented Sensitivity Analysis (QOSA), which can help in analysing the behavior of the response at different quantile levels. Moreover, existing QOSA estimation methods have flaws: bias when input variables are dependent, loss of accuracy and efficiency as input space dimension increases. In this paper, we propose a new estimation procedure of QOSA indices based on the notion of projected random forest, with the initial random forest built from a criterion designed for quantiles: the pinball loss also known as quantile loss.

**Key words:** Quantile Oriented Sensitivity Analysis, random forest, conditional distribution, conditional quantile, projected forest.

## 1 Introduction

Sensitivity Analysis (SA), as defined in Saltelli et al. (2004), is the study of how the uncertainty in the output of a system can be divided and allocated to different sources of uncertainty in its inputs. It is an invaluable tool to understand the behavior of numerical models, determining the most contributing input variables and ascertaining interaction effects within model (see, e.g., Da Veiga et al. (2021) for a recent review). Global Sensitivity Analysis (GSA) focuses on the impact of inputs on some output of interest, when varied over their whole domain (Iooss and Lemaître (2015)). Variance-based methods are well-established and widely used for GSA. In this context, Sobol' indices introduced in Sobol' (1993) are very popular.

In the present paper, we focus on Quantile Oriented Sensitivity Analysis (QOSA), a special case of indices based on contrast functions considered in Fort et al. (2016). QOSA indices are designed to capture the impact of inputs
at a given response quantile level. Browne et al. (2017); Maume-Deschamps and Niang (2018) proposed a kernel based estimator of first-order QOSA indices. Elie-Dit-Cosaque and Maume-Deschamps (2024) introduced an estimator based on random forest. Both approaches require tuning parameters and behave badly as dimension increases. In particular, computation burden of the approach in Elie-Dit-Cosaque and Maume-Deschamps (2024) increases dramatically with input dimension as their approach requires to build a different forest for each QOSA index, in other words, the estimation of the full set of first-order QOSA indices requires to build *p* different forests, with *p* the input space dimension.

In this paper, we introduce a new estimation procedure for QOSA indices. This estimation procedure is mainly based on three ingredients. First we leverage the notion of projected random forest (see Bénard et al. (2022a)) to estimate, for a given quantile level, QOSA indices of any order with only one initial forest. Secondly, to build the initial forest, we use a criterion based on the pinball loss, that Bhat et al. (2015) introduced to build a single quantile oriented decision tree. Finally, in order to avoid overfitting, we use out-of-bag (OOB) samples. Consistency results are proven.

The paper is organized as follows. In Section 2, we recall the definition of QOSA indices and background knowledge on random forests. Section 3 presents the projected random forest idea and consistency results for conditional distribution functions and conditional quantiles. Section 4 is devoted to consistency results for QOSA index estimators. Numerical study and application are shown in Section 5. Some supplemental materials are also provided in Appendix A (technical proofs).

# 2 General framework and tools

In this section, we recall the definition of quantile oriented sensitivity indices, the so-called QOSA indices, of any order. Then, we recall basics on random forests and introduce a new criterion to build quantile oriented random forests.

## 2.1 First-order QOSA indices

We are considering the usual framework in sensitivity analysis, that is $Y = m(\mathbf{X})$ with $Y$ a scalar response, $\mathbf{X} = (X_1, ..., X_p) \in \mathcal{X} \subseteq \mathbb{R}^p$ a $p$-dimensional input vector. One could also consider $Y = m(\mathbf{X}) + \varepsilon$ with $\varepsilon$ a random noise independent on $\mathbf{X}$, this setting is closer to statistical frameworks. The function $m(\cdot)$ is a deterministic model which could be unknown and / or computationally heavy (see Da Veiga et al. (2021) for a more detailed presentation of GSA). Consider input-output samples $\mathcal{D}_n = \{Y^i, \mathbf{X}^i\}_{i=1}^n, \mathbf{X}^i = (X_1^i, X_2^i, ..., X_p^i)$. For any subset $U \subseteq \{1, 2, ..., p\}$ and $\mathbf{x} \in \mathbb{R}^p$, define $\mathbf{x}_U =$
at a given response quantile level. Browne et al. (2017); Maume-Deschamps and Niang (2018) proposed a kernel based estimator of first-order QOSA indices. Elie-Dit-Cosaque and Maume-Deschamps (2024) introduced an estimator based on random forest. Both approaches require tuning parameters and behave badly as dimension increases. In particular, computation burden of the approach in Elie-Dit-Cosaque and Maume-Deschamps (2024) increases dramatically with input dimension as their approach requires to build a different forest for each QOSA index, in other words, the estimation of the full set of first-order QOSA indices requires to build *p* different forests, with *p* the input space dimension.

In this paper, we introduce a new estimation procedure for QOSA indices. This estimation procedure is mainly based on three ingredients. First we leverage the notion of projected random forest (see Bénard et al. (2022a)) to estimate, for a given quantile level, QOSA indices of any order with only one initial forest. Secondly, to build the initial forest, we use a criterion based on the pinball loss, that Bhat et al. (2015) introduced to build a single quantile oriented decision tree. Finally, in order to avoid overfitting, we use out-of-bag (OOB) samples. Consistency results are proven.

The paper is organized as follows. In Section 2, we recall the definition of QOSA indices and background knowledge on random forests. Section 3 presents the Projected Random Forest algorithm and consistency results for the estimation of the conditional cumulative distribution functions and the conditional quantiles. Section 4 is devoted to consistency results for the estimation of QOSA indices. Numerical study and application to a real dataset are shown in Section 5. Some supplemental materials are also provided in Appendix A (technical proofs).

# 2 General framework and tools

In this section, we recall the definition of quantile oriented sensitivity indices, the so-called QOSA indices, of any order. Then, we recall basics on random forests and introduce a new criterion to build quantile oriented random forests.

## 2.1 First-order QOSA indices

We are considering the usual framework in sensitivity analysis, that is $Y = m(\mathbf{X})$ with $Y$ a scalar response, $\mathbf{X} = (X_1, ..., X_p) \in \mathcal{X} \subseteq \mathbb{R}^p$ a $p$-dimensional input vector. One could also consider $Y = m(\mathbf{X}) + \varepsilon$ with $\varepsilon$ a random noise independent on $\mathbf{X}$, this setting is closer to statistical frameworks. The function $m(\cdot)$ is a deterministic model which could be unknown and / or computationally heavy (see Da Veiga et al. (2021) for a more detailed presentation of GSA). Consider input-output samples $\mathcal{D}_n = \{Y^i, \mathbf{X}^i\}_{i=1}^n, \mathbf{X}^i = (X_1^i, X_2^i, ..., X_p^i)$. For any subset $U \subseteq \{1, 2, ..., p\}$ and $\mathbf{x} \in \mathbb{R}^p$, define $\mathbf{x}_U =$
(x_j, j \in U). We shall also consider $\mathcal{X}_U$ the trace of $\mathcal{X}$ in $\mathbb{R}^{|U|}$: $\mathcal{X}_U = \{x \in \mathbb{R}^{|U|} / \exists \tilde{x} \in \mathcal{X}, \tilde{x}_U = x\}$. Our target is to quantify the impact of $X_U$ on the output Y at different quantile levels. Let $\alpha \in (0, 1)$, the QOSA index with respect to the j-th feature at level $\alpha$, is defined as:

$$ S_j^\alpha = \frac{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)] - \mathbb{E}[\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta) | X_j]]}{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)]} \quad (2.1) \\ = \frac{\mathbb{E}[\psi_\alpha(Y, q^\alpha(Y))] - \mathbb{E}[\psi_\alpha(Y, q^\alpha(Y | X_j))]}{\mathbb{E}[\psi_\alpha(Y, q^\alpha(Y))]} = 1 - \frac{\mathbb{E}[\psi_\alpha(Y, q^\alpha(Y | X_j))]}{\mathbb{E}[\psi_\alpha(Y, q^\alpha(Y))]} $$

with $\psi_\alpha(y, \theta) = (y - \theta)(\alpha - \mathbb{I}_{\{y \le \theta\}})$ the pinball function (also known as check function, introduced in Koenker and Hallock (2001)). We denote by $q^\alpha(Y)$ the $\alpha$-level quantile of Y and $q^\alpha(Y | X_j)$ the conditional $\alpha$-level quantile of Y given $X_j$. QOSA indices are a particular case of sensitivity indices based on a contrast function, first introduced in Fort et al. (2016).

## 2.2 Higher order QOSA indices

Following Fort et al. (2016), we define for any $U \subseteq \{1, \dots, p\}$, higher order QOSA indices as:

$$ S_U^\alpha = \frac{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)] - \mathbb{E}[\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta) | X_U]]}{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)]} \quad (2.2) $$

In the particular case of second-order QOSA indices we get:

$$ S_{ij}^\alpha = \frac{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)] - \mathbb{E}[\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta) | X_i, X_j]]}{\min_{\theta \in \mathbb{R}} \mathbb{E}[\psi_\alpha(Y, \theta)]} $$

With $S_U^\alpha$, we measure the contribution on the output Y of the set of features indexed by U, from a quantile perspective. The aim of this paper is to estimate QOSA indices of any order by building a quantile oriented random forest and using its projection for subset U.

Remark that in Kala (2019), another definition for higher order QOSA indices is given. This last definition is similar to higher order Sobol' indices but since in for QOSA indices there is no ANOVA like decomposition, we think the above definition is more relevant.

## 2.3 Random forests

The idea of random forest proposed by Breiman (2001) is to build many randomized CART trees (Breiman et al. (1984)) at first and then average the estimation from different trees as the final estimation. Building different trees from the same dataset requires randomness in the tree growing process. Breiman (2001) proposed to introduce two sources of randomness in the tree construction. The first one is a subsampling step preliminary to each
tree construction. Like the bootstrap strategy from Efron (1992), Breiman (2001) randomly samples n observations with replacement from the original dataset $\mathcal{D}_n$ to construct a new dataset $\mathcal{D}_n^* = ((Y^{*i}, X^{*i}), i = 1, \dots, n)$, then grows the tree with $\mathcal{D}_n^*$. The second randomness source appears during the growing tree process. Instead of optimizing the splitting criterion over all p variables, only part of the features is considered. Each tree is built node by node in a greedy fashion. For each node, the CART splitting criterion is optimized over a subset of features, $\mathcal{M}_{try} \subseteq \{1, 2, \dots, p\}$, randomly selected, with the constraint that each feature has a positive probability to be chosen. For each node $A \subseteq \mathcal{X}$, the best split selects a variable $X_{jA}$ in $\mathcal{M}_{try}$ and a threshold $z_A$ to maximize the CART splitting criterion, whose definition is recalled in (2.3). This splitting criterion measures the decrease of output variance between the parent node and the child nodes. The CART splitting criterion is defined as:

$$L_A^n(j, z) = \frac{1}{N_n^*(A)} \sum_{i=1}^n (Y^i - \bar{Y}_A)^2 \mathbb{I}\{\mathbf{X}^i \in A\} - \\ \frac{1}{N_n^*(A)} \sum_{i=1}^n \left(Y^i - \bar{Y}_{A_L} \mathbb{I}\{X_j^i \le z\} - \bar{Y}_{A_R} \mathbb{I}\{X_j^i > z\}\right)^2 \mathbb{I}\{\mathbf{X}^i \in A\}, \quad (2.3)$$

where for any $j \in \mathcal{M}_{try}$, for any $z \in \mathbb{R}$, the left child node $A_L$ and the right child node $A_R$ are defined as $A_L = \{\mathbf{x} \in A, x_j \le z\}$, $A_R = \{\mathbf{x} \in A, x_j > z\}$, $N_n^*(A) = \#\{i = 1, \dots, n / \mathbf{X}^{*i} \in A\}$ and $\bar{Y}_A$ is the empirical mean of $Y$ on $A$ on the bootstrap sample: $\bar{Y}_A = \frac{1}{N_n^*(A)} \sum_{i=1}^n Y^{*i} \mathbb{I}\{\mathbf{X}^{*i} \in A\}$. To build each tree of a forest composed with B trees, we choose randomly for each tree a bootstrap sample $\mathcal{D}_n^*(b)$. This random choice is modeled by a random vector of indices $\Theta_b^1$. We also select randomly splitting candidate features in each cell of the tree. This random selection is modeled by a random vector of indices $\Theta_b^2$. The vectors $\Theta_b = (\Theta_b^1, \Theta_b^2)$, $b = 1, \dots, B$, are sampled independently from each other and independently from the initial sample $\mathcal{D}_n$. Then, for a new query point $\mathbf{x} \in \mathcal{X}$, $m(\mathbf{x})$ is estimated from the b-th tree as:

$$m_n(\mathbf{x}; \Theta_b, \mathcal{D}_n) = \frac{1}{N_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)} \sum_{i=1}^{n} B_i(\Theta_b^1, \mathcal{D}_n) \mathbb{I}\{\mathbf{X}^i \in A_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)\} Y^i \quad (2.4)$$

with $B_i(\Theta_b^1, \mathcal{D}_n)$ the number of occurrences of $\mathbf{X}^i$ in $\mathcal{D}_n^*(b)$, $A_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)$ the leaf in which $\mathbf{x}$ falls and $N_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)$ the number of elements of $\mathcal{D}_n^*(b)$ in $A_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)$. The estimators from different trees are aggregated as:

$$m_n^B(\mathbf{x}; \Theta_1, \dots, \Theta_b, \mathcal{D}_n) = \frac{1}{B} \sum_{b=1}^{B} m_n(\mathbf{x}; \Theta_b, \mathcal{D}_n). \quad (2.5)$$
space $\mathcal{X} \subseteq \mathbb{R}^p$. We focus here on classification and regression trees (CART) introduced in Breiman et al. (1984). The main principle of decision trees is to recursively partition the input space with splits of the form $X_j \le z$, where $z$ is a real threshold, using the training data $\mathcal{D}_n$. The idea of random forests proposed by Breiman (2001) is to build $B$ randomized trees, with $B \in \mathbb{N}^*$, and then to average the estimates from different trees as the final estimate. Building different trees from the same dataset requires randomness in the tree growing process. Breiman (2001) proposed to introduce two sources of randomness in the tree construction. Let us label each tree by an integer $b \in \{1, \dots, B\}$. The first source of randomness is a subsampling step preliminary to each tree construction. Following the bootstrap strategy from Efron (1992), Breiman (2001) randomly samples $n$ observations with replacement from the original dataset $\mathcal{D}_n$ to construct a new dataset $\mathcal{D}_n^*(b) = ((Y^{*i}, X^{*i}), i = 1, \dots, n)$, then grows the tree with $\mathcal{D}_n^*(b)$. This way $\mathcal{D}_n^*(1), \dots, \mathcal{D}_n^*(B)$ are obtained from $B$ independent realizations $\Theta_1^1, \dots, \Theta_B^1$ of a $n$-dimensional random vector of indices $\Theta^1$ in $\{1, \dots, n\}^n$. The second source of randomness appears during the growing tree process. Instead of optimizing the splitting criterion over all $p$ variables, only part of the features is considered. Each tree is built node by node in a greedy fashion. For each node $A \subseteq \mathcal{X}$, a subset of features $U \subseteq \{1, 2, \dots, p\}$ of cardinality $m$ is randomly selected, with the constraint that each feature has a positive probability to be chosen. Then the CART splitting criterion is optimized over this subset. More precisely, for each node $A \subseteq \mathcal{X}$, the best split selects a variable $X_{j_A}$ in $U$ and a threshold $z_A$ in $\mathbb{R}$ to maximize the CART splitting criterion, whose definition is recalled in (2.4). This splitting criterion measures the decrease of output variance between the parent node and the child nodes. The CART splitting criterion is defined as:

$$L_A^n(j, z) = \frac{1}{N_n^*(A)} \sum_{i=1}^n (Y^i - \bar{Y}_A)^2 \mathbb{I}\{\mathbf{X}^i \in A\} - \frac{1}{N_n^*(A)} \sum_{i=1}^n (Y^i - \bar{Y}_{A_L}) \mathbb{I}\{\mathbf{X}^i \le z\} - \bar{Y}_{A_R} \mathbb{I}\{\mathbf{X}^i > z\} \mathbb{I}\{\mathbf{X}^i \in A\}, \quad (2.4)$$

where for any $j \in U$, for any $z \in \mathbb{R}$, the left child node $A_L$ and the right child node $A_R$ are defined as $A_L = \{\mathbf{x} \in A, x_j \le z\}$, $A_R = \{\mathbf{x} \in A, x_j > z\}$, $N_n^*(A) = \#\{i = 1, \dots, n / \mathbf{X}^{*i} \in A\}$ and $\bar{Y}_A$ is the empirical mean of $Y$ on $A$ computed on the bootstrap sample: $\bar{Y}_A = \frac{1}{N_n^*(A)} \sum_{i=1}^n Y^{*i} \mathbb{I}\{\mathbf{X}^{*i} \in A\}$. For each tree in the forest, let $N_b^{\text{node}}$ denote its number of nodes. The random selection of candidate features is modeled by $B$ independent random vectors $\Theta_b^2 := (\Theta_b^2(1), \dots, \Theta_b^2(N_b^{\text{node}}))$ where $\Theta_b^2(1), \dots, \Theta_b^2(N_b^{\text{node}})$ are independently sampled on the set of vectors of dimension $m$ with components uniformly sampled without replacement on $\{1, \dots, p\}$. Also, $\Theta_1^2, \dots, \Theta_B^2$ and
Equation (2.5) rewrites as:

$$m_n^B(x; \Theta_1, \dots, \Theta_b, \mathcal{D}_n) = \sum_{i=1}^n w_i(x) Y^i, \qquad (2.6)$$

with

$$w_i(\mathbf{x}) = \frac{1}{B} \sum_{b=1}^{B} \frac{B_i(\Theta_b^1, \mathcal{D}_n) \mathbb{I}_{\{\mathbf{X}^i \in A_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)\}}}{N_n(\mathbf{x}; \Theta_b, \mathcal{D}_n)}. \qquad (2.7)$$

Meinshausen and Ridgeway (2006) proposed to replace Y with $\mathbb{I}_{\{Y \le y\}}$ in (2.6) in order to estimate the conditional distribution:

$$\hat{F}_{Y|\mathbf{X}=\mathbf{x}}(y) = \sum_{i=1}^{n} w_i(\mathbf{x}) \mathbb{I}_{\{Y^i \le y\}}.$$

Finally, the conditional quantile function can be estimated by the generalized inverse of $\hat{F}_{Y|\mathbf{X}=\mathbf{x}}(\cdot)$.

However CART-split criterion is sensitive to changes in conditional mean of Y given $\mathbf{X}$ and is not designed for changes in conditional quantiles. It is one of the reasons why alternative splitting criteria, more focused on quantiles, have been introduced. Bhat et al. (2015) uses the pinball loss to grow one quantile regression tree. They also propose an online update algorithm to improve the search efficiency. Athey et al. (2019) proposes a general framework, gradient tree, which could be applied for different regression tasks. In the quantile regression case, their splitting criterion is similar to Gini impurity for classification problem. Cévid et al. (2022) consider the Maximum Mean Discrepancy (MMD), a distributional metric between left and right child node. They propose a fast approximation algorithm. Because the pinball loss is the quantile loss function, we generalize Bhat et al. (2015)'s splitting criterion to build a random forest. More precisely, we would like to maximize the quantity:

$$E[\psi_{\theta}(Y)|X \in A] - P[X \in A_L | X \in A] E[\psi_{\alpha}(Y, \theta_L)|X \in A_L] \\ - P[X \in A_R | X \in A] E[\psi_{\alpha}(Y, \theta_R)|X \in A_R], \qquad (2.8)$$

with $\psi_{\alpha}(y, \theta) = (y - \theta)(\alpha - \mathbb{I}_{\{y < \theta\}})$, $\alpha \in (0, 1)$, $\theta_L$ and $\theta_R$ the $\alpha$-level quantile in left and right child node respectively. In practice, we optimize the empirical form of (2.8). Note that the first term in (2.8) can be ignored in the optimization because it does not depend on $A_L$ and $A_R$. Thus we shall minimize the empirical counterpart of the last two terms in (2.8):

$$\frac{n_L}{n_P} \frac{1}{n_L} \sum_{X_i \in A_L} \psi_{\alpha}(Y^i, \hat{\theta}_L) + \frac{n_R}{n_P} \frac{1}{n_R} \sum_{X_i \in A_R} \psi_{\alpha}(Y^i, \hat{\theta}_R), \qquad (2.9)$$

where $n_L$ and $n_R$ are the sample size in left and right child nodes respectively, $\hat{\theta}_L$ and $\hat{\theta}_R$ are the quantile estimation in each child node.
$\Theta_1^1, \dots, \Theta_B^1$ are independent from each other and independent from the initial sample $\mathcal{D}_n$. Then, for a new query point $\boldsymbol{x} \in \mathcal{X}$, $m(\boldsymbol{x})$ is estimated from the b-th tree as:

$$m_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n) = \frac{1}{N_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)} \sum_{i=1}^{n} B_i(\Theta_b^1, \mathcal{D}_n) \mathbb{I}_{\{\boldsymbol{X}^i \in A_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)\}} Y^i \quad (2.5)$$

with $\Theta_b = (\Theta_b^1, \Theta_b^2)$, $B_i(\Theta_b^1, \mathcal{D}_n)$ the number of occurrences of $\boldsymbol{X}^i$ in $\mathcal{D}_n^*(b)$, $A_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)$ the leaf in which $\boldsymbol{x}$ falls and $N_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)$ the number of elements of $\mathcal{D}_n^*(b)$ in $A_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)$. The estimators from different trees are aggregated as:

$$m_n^B(\boldsymbol{x}; \Theta_1, \dots, \Theta_b, \mathcal{D}_n) = \frac{1}{B} \sum_{b=1}^{B} m_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n). \quad (2.6)$$

Equation (2.6) rewrites as:

$$m_n^B(\boldsymbol{x}; \Theta_1, \dots, \Theta_b, \mathcal{D}_n) = \sum_{i=1}^{n} w_i(\boldsymbol{x}) Y^i, \quad (2.7)$$

with

$$w_i(\boldsymbol{x}) = \frac{1}{B} \sum_{b=1}^{B} \frac{B_i(\Theta_b^1, \mathcal{D}_n) \mathbb{I}_{\{\boldsymbol{X}^i \in A_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)\}}}{N_n(\boldsymbol{x}; \Theta_b, \mathcal{D}_n)}. \quad (2.8)$$

Meinshausen and Ridgeway (2006) proposed to replace $Y$ with $\mathbb{I}_{\{Y \le y\}}$ in (2.7) in order to estimate the conditional distribution:

$$\hat{F}_{Y|\boldsymbol{X}=\boldsymbol{x}}(y) = \sum_{i=1}^{n} w_i(\boldsymbol{x}) \mathbb{I}_{\{Y^i \le y\}}.$$

Finally, the conditional quantile function can be estimated by the generalized inverse of $\hat{F}_{Y|\boldsymbol{X}=\boldsymbol{x}}(\cdot)$. Recall that for $F$ a cumulative distribution function, the generalized inverse of $F$ is defined as $F^{-}(\alpha) = \inf\{x : F(x) \ge \alpha\}$, for any $\alpha \in [0, 1]$.

**Beyond CART splitting criterion** The major drawback of CART splitting criterion is that, although sensitive to changes in conditional mean of $Y$ given $\boldsymbol{X}$, it may ignore changes in conditional quantiles. It is one of the reasons why alternative splitting criteria, more focused on quantiles, have been introduced. Bhat et al. (2015) uses the pinball loss to grow one quantile regression tree. They also propose an online update algorithm to improve the search efficiency. Athey et al. (2019) proposes a general framework, gradient tree, which could be applied for different regression tasks. In the quantile regression case, their splitting criterion is similar to Gini impurity
In the next two sections we explain how the projection algorithm introduced in Bénard et al. (2022a) can be applied on our forest to estimate the conditional distribution function of the output and conditional quantiles.

To conclude this section, let us show a simple tree structure and introduce some notation for the rest of the paper. As shown in Fig.1, there are

<figure><img src="image_3.png" /><figcaption>Figure 1: Decision tree structure</figcaption></figure>

two kinds of nodes in a tree. Terminal leaf nodes are denoted by $A_i$. The other nodes are labeled by a pair of indices $(l, k)$, where $l$ is the tree level and $k$ is the index of each node at $l$-level. We denote $i(l, k) \in \{1, 2, ..., p\}$ and $a_{l,k} \in \mathbb{R}$ the index of the feature and the corresponding threshold selected by optimizing the splitting criterion (2.9) at node $(l, k)$.

# 3 Projected random forests

In this section, we present the notion of projected forest in Bénard et al. (2022a) which uses an idea from Lundberg and Lee (2017). The key idea is to build one random forest and then to project it adequately to estimate the indices $S_U^\alpha$. This reduces computational burden in high dimension. The algorithm is available for all $U \subseteq \{1, 2, ..., p\}$. The details can be found in Algo. 1. The projected idea is as in Bénard et al. (2022a). The conditional quantile estimation is then similar to the random forest estimation using the projected leaves concept. So, we only display the core part of projecting in Algo. 1. The illustration in one dimension is shown in Fig. 2.

## 3.1 Estimation of conditional cumulative distribution functions

We shall need the following properties on tree construction. These construction hypotheses are also done in Athey et al. (2019) and Bénard et al. (2022a).

(P1) [Data sampling] Rather than using bootstrap sampling with replacement, we use subsampling with growing size $s_n$ such that $s_n/n = \kappa$,
for classification problem. Cévid et al. (2022) consider the Maximum Mean Discrepancy (MMD), a distributional metric between left and right child node. They propose a fast approximation algorithm. Because the pinball loss is the quantile loss function, we generalize Bhat et al. (2015)'s splitting criterion to build a random forest. More precisely, we would like to maximize the quantity:

$$E[\psi_{\alpha}(Y, \theta)|X \in A] - P[X \in A_L | X \in A]E[\psi_{\alpha}(Y, \theta_L)|X \in A_L] - P[X \in A_R | X \in A]E[\psi_{\alpha}(Y, \theta_R)|X \in A_R], \quad (2.9)$$

with $\psi_{\alpha}(y, \theta) = (y - \theta)(\alpha - \mathbb{I}_{\{y<\theta\}})$, $\alpha \in (0, 1)$, $\theta_L$ and $\theta_R$ the $\alpha$-level quantile in left and right child nodes, respectively. In practice, we optimize the empirical form of (2.9). Note that the first term in (2.9) can be ignored in the optimization because it does not depend on $A_L$ and $A_R$. Thus we shall minimize the empirical counterpart of the last two terms in (2.9):

$$\frac{n_L}{n_P} \frac{1}{n_L} \sum_{X_i \in A_L} \psi_{\alpha}(Y^i, \hat{\theta}_L) + \frac{n_R}{n_P} \frac{1}{n_R} \sum_{X_i \in A_R} \psi_{\alpha}(Y^i, \hat{\theta}_R), \quad (2.10)$$

where $n_L$ and $n_R$ are the sample size in left and right child nodes respectively, $\hat{\theta}_L$ and $\hat{\theta}_R$ are the quantile estimation in each child node.

Useful notation: There are two kinds of nodes in a tree. Terminal leaf nodes are denoted by $A_i$. The other nodes are labeled by a pair of indices $(l, k)$, where $l$ is the tree level and $k$ is the index of each node at $l$-level. We denote $i(l, k) \in \{1, 2, ..., p\}$ and $a_{l,k} \in \mathbb{R}$ the index of the feature and the corresponding threshold selected by optimizing the splitting criterion in (2.10) at node $(l, k)$. This notation is illustrated on Fig. 1 below.

<figure><img src="image_2.png" /><figcaption>Figure 1: Decision tree structure</figcaption></figure>

## 3 Projected random forests

The Projected Random Forest algorithm was introduced in Bénard et al. (2022a). The algorithm is stated in Section 3. Then we leverage it to esti-
<caption>**Algorithm 1:** Projected Random Forest (Bénard et al. (2022a))</caption>

**Inputs** : A random forest, $\mathcal{T}_1, \mathcal{T}_2, ..., \mathcal{T}_B$, fit with $\mathcal{D}_n$,
A variable subset $U \subseteq \{1, 2, ..., p\}$,
A query point $\boldsymbol{x} = (x_1, x_2, ..., x_p)$.

**Outputs**: Projected leafs: $\mathcal{A}(\boldsymbol{x}|\Theta_1, U), \mathcal{A}(\boldsymbol{x}|\Theta_2, U), ..., \mathcal{A}(\boldsymbol{x}|\Theta_B, U)$

1 **for** all trees in the forest **do**
| **/* Step 1: initialize variables */**
| 2 initialize $\mathit{nodes\_level}$ as a list of nodes containing only the root node;
| 3 initialize $\mathit{nodes\_child}$ as an empty list of child nodes;
| 4 initialize $\mathit{samples}$ as the list of observation indices of the full training data of the tree;
| 5 **for** all levels in the tree **do**
|   | **/* Step 2: drop $x_U$ to the next tree level with the relevant training samples */**
|   | 6 **for** all nodes in $\mathit{nodes\_level}$ **do**
|   |   | **if** the node splits on a variable in $U$, **then**
|   |   |   | compute whether $x_U$ falls in the left or right child node;
|   |   |   | append the child node to $\mathit{nodes\_child}$;
|   |   |   | set $\mathit{samples\_child}$ as the observations in $\mathit{samples}$ which satisfy the split
|   |   |   | **else**
|   |   |   | append both the left and right children nodes to $\mathit{nodes\_child}$;
|   |   |   | **end**
|   |   |   | **if** the size of $\mathit{samples\_child}$ is lower then
|   |   |   |   | **min\_node\_size** **then**
|   |   |   |   |   | break the loop through the tree levels;
|   |   |   | **else**
|   |   |   |   | set $\mathit{samples} = \mathit{samples\_child}$;
|   |   |   | **end**
|   |   |   | **end**
|   |   |   | **end**
|   |   |   | set $\mathit{nodes\_level} = \mathit{nodes\_child}$;
|   |   |   | **end**
|   |   |   | $\mathcal{A}(\boldsymbol{x}|\Theta_b, U) = \mathit{samples}$;
|   |   |   | **end**
mate, for each subset $U \subseteq \{1, \dots, p\}$, the conditional cumulation distribution function of the output $Y$ (see Section 3.2) and its conditional quantile (see Section 3.3). These are essential tools for the estimation of QOSA indices in Section 4.

## 3.1 Projected Random Forest algorithm

The Projected Random Forest algorithm takes as inputs the initial forest and a given subset $U \subseteq \{1, \dots, p\}$. Then, its general principle is to project the partition of each tree of the forest on the subspace spanned by the variables in $U$, as illustrated in Fig. 2 below.

<figure><img src="image_3.png" /><figcaption>Figure 2: Example in dimension $p = 2$ and with $U = \{1\}$. The partition of $[0, 1]^2$ given by a tree in the original forest (top left) is projected on the subspace spanned by $X_1$ (top right). The tree is plotted (bottom left) and the evolution of the projected partition associated to point $x$ in Algorithm 1 is provided (bottom right).</figcaption></figure>

The Projected Random Forest algorithm, as introduced in Bénard et al. (2022a), is stated below as Algorithm 1.
<figure><img src="image_4.png" /><figcaption>Figure 2: Projected algorithm illustration for $p = 2, U = \{1\}$. Left panel represents the structure of the original tree and partition in $\mathbb{R}^2$. Each dash line represents a node in the tree. Right panel is the projected leaf, which could be tracked by level.</figcaption></figure>
<caption>**Algorithm 1:** Projected Random Forest (Bénard et al. (2022a))</caption>

**Inputs** : A random forest, $\mathcal{T}_1, \mathcal{T}_2, ..., \mathcal{T}_B$, fit with $\mathcal{D}_n$,
A variable subset $U \subseteq \{1, 2, ..., p\}$,
A query point $\boldsymbol{x} = (x_1, x_2, ..., x_p)$.

**Outputs**: Projected leafs: $\mathcal{A}(\boldsymbol{x}|\Theta_1, U), \mathcal{A}(\boldsymbol{x}|\Theta_2, U), ..., \mathcal{A}(\boldsymbol{x}|\Theta_B, U)$

1 **for** all trees in the forest **do**
| **/* Step 1: initialize variables */**
| 2 initialize $\text{nodes\_level}$ as a list of nodes containing only the root node;
| 3 initialize $\text{nodes\_child}$ as an empty list of child nodes;
| 4 initialize $\text{samples}$ as the list of observation indices of the full training data of the tree;
| 5 **for** all levels in the tree **do**
| | **/* Step 2: drop $x_U$ to the next tree level with the relevant training samples */**
| | 6 **for** all nodes in $\text{nodes\_level}$ **do**
| | | if the node splits on a variable in $U$, then
| | | 7 | compute whether $x_U$ falls in the left or right child node;
| | | 8 | append the child node to $\text{nodes\_child}$;
| | | 9 | set $\text{samples\_child}$ as the observations in $\text{samples}$ which satisfy the split
| | | 10 | else
| | | 11 | append both the left and right children nodes to $\text{nodes\_child}$;
| | 12 | end
| | 13 | if the size of $\text{samples\_child}$ is lower then $\min_{\text{node\_size}}$
| | 14 | | then
| | | 15 | break the loop through the tree levels;
| | 16 | else
| | 17 | | set $\text{samples} = \text{samples\_child}$;
| | 18 | end
| | 19 | end
| | 20 | set $\text{nodes\_level} = \text{nodes\_child}$;
| | 21 | end
| | 22 | $\mathcal{A}(\boldsymbol{x}|\Theta_b, U) = \text{samples}$;
| | 23 | end
where $\kappa \in (0, 1)$ is a fixed number. For convenience, we reuse the notation $D_n^*(b)$ for subsamples.

(P2) [γ-regularity] Each split leaves at least a fraction $\gamma \in (0, 0.5)$ of the available training sample on each side.

(P3) [Random-split] At each tree node, the number mtry of candidate variables drawn to optimize the split is set to mtry = 1 with a small probability $\pi > 0$ and for mtry = 1 and any $j = 1, \dots, p$ the probability that the split occurs along feature $X_j$ is set to $\pi/p$.

Let us consider the following estimator of the conditional distribution function $y \mapsto F_Y(y|\mathbf{X}_U) = \mathbb{P}(Y \le y|\mathbf{X}_U)$: for $\mathbf{x} \in \mathcal{X}_U$, $U \subseteq \{1, \dots, p\}$

$$w_i(\mathbf{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}_{\{\mathbf{X}^i \in A_n(\mathbf{x}; D_n^*(b), U)\}}}{N(\mathbf{x}; D_n^*(b), D_n^*(b), U)},$$

$$\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \sum_{i=1}^{n} w_i(\mathbf{x}|U) \mathbb{I}_{(Y^i \le y)},$$

where $A_n(\mathbf{x}; D_n^*(b), U)$ is the projected leaf containing $\mathbf{x}$ and $N(\mathbf{x}; D_n^*(b), D_n^*(b), U)$ is the number of elements in $D_n^*(b)$ which belong to $A_n(\mathbf{x}; D_n^*(b), U)$. In what follows, we will refer to $\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x})$ as the standard estimator. Also, for any $j \in \{1, \dots, n\}$, we denote by $\hat{F}_{s_n,n-1}(y|\mathbf{X}_U^j)$ the standard estimator built with $D_n \setminus (\mathbf{X}_U^j, Y^j)$ and with sample size $a_n$.

**Assumption 1.** The covariate vector $X = (X_1, \dots, X_p)$ admits a density over $\mathcal{X}$ bounded from below by a strictly positive constant.

**Assumption 2** (sample size in a projected leaf). For fixed $\beta > 1$ and $C > 0$, for any $U \subseteq \{1, 2, \dots, p\}$, and $\mathbf{x}$ in $\mathcal{X}_U$, $N(\mathbf{x}; D_n^*(\Theta), D_n^*(\Theta), U) \ge C\sqrt{n}(\ln n)^\beta$.

The number of elements in a leaf is an hyper-parameter in the building of the forest, so it can be controlled.

**Assumption 3.** Assumption on $B$ (number of trees), $B = O(n^\beta)$, for some fixed $\beta > 0$.

**Assumption 4.** For any $U \subseteq \{1, 2, \dots, p\}$ and $\mathbf{x}$ in $\mathcal{X}_U$, the conditional distribution function $F_Y(\cdot|\mathbf{X}_U = \mathbf{x})$ is continuous and increasing, with the conditional density function $f_Y(\cdot|\mathbf{X}_U = \mathbf{x})$ continuous. For any $y \in \mathbb{R}$, $F_Y(y|\mathbf{X}_U = \cdot)$ is continuous.

Note that properties (P1) to (P3) as far as Assumptions 1 to 4 are standard assumptions in the literature (see, e.g., Bénard et al. (2022b); Elie-Dit-Cosaque and Maume-Deschamps (2022b)).
## 3.2 Conditional cumulative distribution function estimation

In this section, we leverage Algorithm 1 to estimate conditional cumulative distribution functions. We first state the following assumptions on tree construction. These hypotheses are also done in Athey et al. (2019) and Bénard et al. (2022a).

(P1) [Data sampling] Rather than using bootstrap sampling with replacement, we use subsampling with growing size $s_n$ such that $s_n/n = \kappa$, where $\kappa \in (0, 1)$ is a fixed number. For convenience, we reuse the notation $D_n^*(b)$ for subsamples.

(P2) [γ-regularity] Each split leaves at least a fraction $\gamma \in (0, 0.5)$ of the available training sample on each side.

(P3) [Random-split] At each tree node, the number mtry of candidate variables drawn to optimize the split is set to mtry = 1 with a small probability $\pi > 0$ and for mtry = 1 and any $j = 1, ..., p$ the probability that the split occurs along feature $X_j$ is set to $\pi/p$.

Then we consider the following estimator of the conditional cumulative distribution function $y \mapsto F_Y(y|\mathbf{X}_U) = \mathbb{P}(Y \le y|\mathbf{X}_U)$: for $\mathbf{x} \in \mathcal{X}_U, U \subseteq \{1, ..., p\}$,

$$ \hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \sum_{i=1}^{n} w_i(\mathbf{x}|U) \mathbb{I}_{(Y^i \le y)} $$

with

$$ w_i(\mathbf{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}_{\{\mathbf{X}^i \in A_n(\mathbf{x}; D_n^*(b), U)\}}}{N(\mathbf{x}; D_n^*(b), D_n^*(b), U)} $$

where $A_n(\mathbf{x}; D_n^*(b), U)$ is the projected leaf containing $\mathbf{x}$ and $N(\mathbf{x}; D_n^*(b), D_n^*(b), U)$ the number of elements in $D_n^*(b)$ belonging to $A_n(\mathbf{x}; D_n^*(b), U)$. In what follows, we will refer to $\hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x})$ as the standard estimator. Also, for any $j \in \{1, ..., n\}$, we denote by $\hat{F}_{s_n,n-1}(y|\mathbf{X}_U^j)$ the standard estimator built with $D_n \setminus (\mathbf{X}_U^j, Y^j)$ and with sample size $s_n$.

We also state a set of additional assumptions.

**Assumption 1.** The covariate vector $\mathbf{X} = (X_1, ..., X_p)$ admits a density over $\mathcal{X}$ bounded from below by a strictly positive constant.

**Assumption 2** (sample size in a projected leaf). For fixed $\beta > 1$ and $C > 0$, for any $U \subseteq \{1, 2, ..., p\}$, and $\mathbf{x}$ in $\mathcal{X}_U$, $N(\mathbf{x}; D_n^*(\Theta), D_n^*(\Theta), U) \ge C\sqrt{n}(\ln n)^\beta$.

Assumption 2 concerns the number of elements in a leaf, this hyper-parameter can be controlled in the algorithm used to grow a forest.
One of the main ingredient in the consistency proof of QOSA indices is Theorem 3.3 below. It uses as a technical tool, a dummy estimator constructed with an additional sample, whose consistency is stated in Lemma 3.1. Consider an additional sample $\mathcal{D}_n^\diamond = ((Y^{\diamond i}, \mathbf{X}^{\diamond i}), i = 1, \dots, n)$, independent of $\mathcal{D}_n$ with the $(Y^{\diamond i}, \mathbf{X}^{\diamond i})$'s independent and distributed as $(Y, \mathbf{X})$. The dummy estimator is defined as

$$ \hat{F}^{\diamond}(y | \mathbf{X}_U = \mathbf{x}) = \sum_{i=1}^{n} w_i^{\diamond}(\mathbf{x}|U) \mathbb{I}_{(Y^{\diamond i} \le y)}, $$

where

$$ w_i^{\diamond}(\mathbf{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}_{\{\mathbf{X}^{\diamond i} \in A_n(\mathbf{x}; \mathcal{D}_n^*(b), U)\}}}{N(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^\diamond, U)}, \quad (3.1) $$

$A_n(\mathbf{x}; \mathcal{D}_n^*(b), U)$ is as above and $N(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^\diamond, U)$ is the number of elements in $\mathcal{D}_n^\diamond$ which belong to $A_n(\mathbf{x}; \mathcal{D}_n^*(b), U)$.

The following two lemmas are key ingredients in the proof of Theorem 3.3.

**Lemma 3.1.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, \dots, p\}$,

$$ \forall \mathbf{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, |\hat{F}^{\diamond}(y|\mathbf{X}_U = \mathbf{x}) - F(y|\mathbf{X}_U = \mathbf{x})| \xrightarrow{n \to \infty} 0. $$

**Lemma 3.2.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, \dots, p\}$,

$$ \forall \mathbf{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, |\hat{F}^{\diamond}(y|\mathbf{X}_U = \mathbf{x}) - \hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x})| \xrightarrow{n \to \infty} 0. $$

The proof of Lemmas 3.1 and 3.2 is the same as the one of Proposition 6.1 and Lemma 6.2 in Elie-Dit-Cosaque and Maume-Deschamps (2022b) using Corollary A.2 and Lemma A.3 in the Appendix.

**Theorem 3.3.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, \dots, p\}$,

$$ \forall \mathbf{x} \in \mathcal{X}_U, \sup_{y \in \mathbb{R}} |\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) - F(y|\mathbf{X}_U = \mathbf{x})| \xrightarrow{n \to \infty} 0. $$

*Proof.* The proof idea is similar to the one of Theorem 4.3 in Elie-Dit-Cosaque and Maume-Deschamps (2022b). The consistency result for the standard estimator is obtained by combining Lemmas 3.1 and 3.2 above:

$$ \forall \mathbf{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, \hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) \xrightarrow{n \to \infty} F(y|\mathbf{X}_U = \mathbf{x}). $$

Then, Dini's second theorem, leads to the uniform a.s. convergence under Assumption 4. □
**Assumption 3.** *Assumption on B (number of trees), B = O(n<sup>β</sup>), for some fixed β > 0.*

**Assumption 4.** For any $U \subseteq \{1, 2, ..., p\}$ and $\boldsymbol{x}$ in $\mathcal{X}_U$, the conditional distribution function $F_Y(\cdot|\boldsymbol{X}_U = \boldsymbol{x})$ is continuous and increasing, with the conditional density function $f_Y(\cdot|\boldsymbol{X}_U = \boldsymbol{x})$ continuous. For any $y \in \mathbb{R}$, $F_Y(y|\boldsymbol{X}_U = \cdot)$ is continuous.

Assumptions 1 to 4, as long as Properties (P1) to (P3), are standard assumptions in the literature (see, e.g., Bénard et al. (2022b); Elie-Dit-Cosaque and Maume-Deschamps (2022b)).

One of the main ingredient in the proof consistency of QOSA indices is the consistency of the conditional cumulative distribution function estimator. This result is stated in Theorem 3.3 below. It uses as a technical tool, a dummy estimator constructed with an additional sample, whose consistency is stated in Lemma 3.1. Consider an additional sample $\mathcal{D}_n^\diamond = ((Y^{\diamond i}, \mathbf{X}^{\diamond i}), i = 1, \dots, n)$, independent of $\mathcal{D}_n$ with the $(Y^{\diamond i}, \mathbf{X}^{\diamond i})$'s independent and distributed as $(Y, \mathbf{X})$. The dummy estimator is defined as

$$ \hat{F}^{\diamond}(y|\mathbf{X}_U = \boldsymbol{x}) = \sum_{i=1}^{n} w_i^{\diamond}(\boldsymbol{x}|U) \mathbb{I}_{(Y^{\diamond i} \le y)} $$

with

$$ w_i^{\diamond}(\boldsymbol{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}_{\{\boldsymbol{X}^{\diamond i} \in A_n(\boldsymbol{x}; \mathcal{D}_n^*(b), U)\}}}{N(\boldsymbol{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^\diamond, U)} \quad (3.1) $$

where $A_n(\boldsymbol{x}; \mathcal{D}_n^*(b), U)$ is as above and $N(\boldsymbol{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^\diamond, U)$ is the number of elements in $\mathcal{D}_n^\diamond$ belonging to $A_n(\boldsymbol{x}; \mathcal{D}_n^*(b), U)$. Lemmas 3.1 and 3.2 are key ingredients in the proof of Theorem 3.3.

**Lemma 3.1.** *Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, ..., p\}$,*

$$ \forall \boldsymbol{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, |\hat{F}^{\diamond}(y|\mathbf{X}_U = \boldsymbol{x}) - F(y|\mathbf{X}_U = \boldsymbol{x})| \xrightarrow{n \to \infty} 0. $$

**Lemma 3.2.** *Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, ..., p\}$,*

$$ \forall \boldsymbol{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, |\hat{F}^{\diamond}(y|\mathbf{X}_U = \boldsymbol{x}) - \hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \boldsymbol{x})| \xrightarrow{n \to \infty} 0. $$

The proof of Lemmas 3.1 and 3.2 is the same as the one of Proposition 6.1 and Lemma 6.2 in Elie-Dit-Cosaque and Maume-Deschamps (2022b), using Corollary A.2 and Lemma A.3 in the Appendix.

**Theorem 3.3.** *Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $U \subseteq \{1, 2, ..., p\}$,*

$$ \forall \boldsymbol{x} \in \mathcal{X}_U, \sup_{y \in \mathbb{R}} |\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \boldsymbol{x}) - F(y|\mathbf{X}_U = \boldsymbol{x})| \xrightarrow{n \to \infty} 0. $$
**Remark:** For $U = \{1, 2, ..., p\}$, the projected random forest will return the original forest. So the consistency result of the original random forest is included in Theorem 3.3.

Next, we introduce the OOB estimator (Out Of Bag estimator). Recall that for each tree a subsample of size $s_n$ is chosen, with respect to $\Theta_b$. So that, all indices $i \in \{1, ..., n\}$ will not participate in the tree construction. Let us consider the random set $\Lambda_{n,i}$,

$$ \Lambda_{n,i} = \{b \in \{1, 2, ..., B\} : (Y^i, \mathbf{X}^i) \notin \mathcal{D}_n^*(\Theta_b)\}, $$

it is the set of trees built without using the i-th element of $\mathcal{D}_n$. Then, for all $j = 1, ..., n$, for all $U \subseteq \{1, ..., p\}$, the OOB estimator is defined as

$$ \hat{F}^{OOB}(y|\mathbf{X}_U^j) = \sum_{i=1}^{n} w_i^{OOB}(\mathbf{X}^j|U) \mathbb{I}_{\{Y^i \le y\}}, $$

where $w_i^{OOB}(\mathbf{X}^j|U) = 0$ if $\Lambda_{n,i} = \emptyset$ and otherwise,

$$ w_i^{OOB}(\mathbf{X}^j|U) = \frac{\mathbb{I}_{\{|\Lambda_{n,j}|>0\}}}{|\Lambda_{n,j}|} \sum_{b \in \Lambda_{n,j}} \frac{\mathbb{I}_{\{\mathbf{X}^i \in A(\mathbf{X}^j; \mathcal{D}_n^*(\Theta), U)\}}}{N(\mathbf{X}^j; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)}. $$

The $L^2$ consistency of the OOB estimator is stated in Theorem 3.5 below. It uses Lemma 3.4 which shows that the OOB error can be controlled by the standard one and which proof is postponed to Appendix A.

**Lemma 3.4.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $B \in \mathbb{N}^*, j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \\ \le \frac{2}{1-s_n/n} \sup_y \mathbb{E}[(\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] + O((s_n/n)^B). $$

**Theorem 3.5.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then $\forall j \in \{1, ..., n\}, \forall U \subseteq \{1, 2, ..., p\}$,

$$ \lim_{n \to \infty} \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \to 0. $$

*Proof.* The proof follows directly from Lemma 3.4 which provides a relation between the OOB error and the standard error and Theorem 3.3. Remark that $\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j)$ and $F(y|\mathbf{X}_U^j)$ are distribution functions, so that their values are in $[0, 1]$ and the a.s. consistency from Theorem 3.3 leads to $L^2$ consistency. $\square$
Proof. The proof idea is similar to the one of Theorem 4.3 in Elie-Dit-Cosaque and Maume-Deschamps (2022b). The consistency result for the standard estimator is obtained by combining Lemmas 3.1 and 3.2 above:

$$ \forall \boldsymbol{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, \hat{F}_{B,s_n,n}(y|\boldsymbol{X}_U = \boldsymbol{x}) \xrightarrow[n \to \infty]{\text{a.s.}} F(y|\boldsymbol{X}_U = \boldsymbol{x}). $$

Then, Dini's second theorem, leads to the uniform a.s. convergence under Assumption 4. □

**Remark:** For $U = \{1, 2, ..., p\}$, the projected random forest will return the original forest. So the consistency result of the original random forest is included in Theorem 3.3.

Next, in order to avoid overfitting, we introduce the OOB estimator (Out Of Bag estimator). Recall that for each tree a subsample of size $s_n$ is chosen, with respect to $\Theta_b$. As a consequence, not all indices $i \in \{1, ..., n\}$ will participate in the tree construction. Let us consider the random set $\Lambda_{n,i}$ defined as

$$ \Lambda_{n,i} = \{b \in \{1, 2, ..., B\} : (Y^i, \mathbf{X}^i) \notin \mathcal{D}_n^*(\Theta_b)\}. $$

It corresponds to the set of trees built without using the i-th element of $\mathcal{D}_n$. Then, for all $j = 1, ..., n$, for all $U \subseteq \{1, ..., p\}$, the OOB estimator is defined as

$$ \hat{F}^{OOB}(y|\mathbf{X}_U^j) = \sum_{i=1}^{n} w_i^{OOB}(\mathbf{X}^j|U) \mathbb{I}_{\{Y^i \le y\}} $$

where $w_i^{OOB}(\mathbf{X}^j|U) = 0$ if $\Lambda_{n,i} = \emptyset$ and otherwise,

$$ w_i^{OOB}(\mathbf{X}^j|U) = \frac{\mathbb{I}_{\{|\Lambda_{n,j}|>0\}}}{|\Lambda_{n,j}|} \sum_{b \in \Lambda_{n,j}} \frac{\mathbb{I}_{\{\mathbf{X}^i \in A(\mathbf{X}^j; \mathcal{D}_n^*(\Theta), U)\}}}{N(\mathbf{X}^j; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)}. $$

The $L^2$ consistency of the OOB estimator is stated in Theorem 3.5 below. It uses Lemma 3.4 which shows that the OOB error can be controlled by the standard one and which proof is postponed to Appendix A.

**Lemma 3.4.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $B \in \mathbb{N}^*, j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \\ \le \frac{2}{1-s_n/n} \sup_y \mathbb{E}[(\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] + O((s_n/n)^B). $$

**Theorem 3.5.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then $\forall j \in \{1, ..., n\}, \forall U \subseteq \{1, 2, ..., p\}$,

$$ \lim_{n \to \infty} \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \to 0. $$
## 3.2 Conditional quantile estimation

In what follows, we are interested in the estimation of the conditional quantiles at level $\alpha \in (0, 1)$ fixed:

$$q^{\alpha}(Y|\mathbf{X}_U) = \inf_{y \in \mathbb{R}} \{y, F_Y(y|\mathbf{X}_U) \geq \alpha\}.$$

We shall use the OOB estimator of the conditional distribution function to get an estimation of $q^{\alpha}(Y|\mathbf{X}_U)$:

$$\hat{Q}_{n}^{OOB}(\alpha|\mathbf{X}_{U}=\mathbf{x})=\inf_{y \in \mathbb{R}}\{y, \hat{F}^{OOB}(y|\mathbf{X}_{U}=\mathbf{x}) \geq \alpha\}.$$

It is easily noticed that

$$\hat{Q}_{n}^{OOB}(\alpha|\mathbf{X}_{U}=\mathbf{x})=\inf_{k=1, \ldots, n}\{Y^{k}, \hat{F}^{OOB}(Y^{k}|\mathbf{X}_{U}=\mathbf{x}) \geq \alpha\}.$$

The standard estimator $\hat{Q}_{B,s_n,n}(\alpha|\mathbf{X}_U = \mathbf{x})$ of the conditional quantile is defined in the same way replacing $\hat{F}^{OOB}$ with $\hat{F}_{B,s_n,n}$.

In the following, in order to lighten the notation, for $\mathbf{x} \in \mathcal{X}_U$ we shall write $Q^{\alpha}(\mathbf{x})$ for $Q^{\alpha}(Y|\mathbf{X}_U = \mathbf{x})$ and $\hat{Q}_n(\alpha|\mathbf{x})$, for an estimator of $Q^{\alpha}(Y|\mathbf{X}_U = \mathbf{x})$.

**Theorem 3.6.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then standard and OOB estimators of the conditional quantile are consistent. That is, if $\mathbf{x} \in \mathcal{X}_U$,

$$\forall \alpha \in (0,1), \hat{Q}_{B,s_n,n}(\alpha|\mathbf{X}_U = \mathbf{x}) \xrightarrow{n \to \infty} Q(\alpha|\mathbf{X}_U = \mathbf{x}), \quad (3.2)$$

$$\forall j \in \{1, \dots, n\}, \forall \alpha \in (0,1), \hat{Q}_{n}^{OOB}(\alpha|\mathbf{X}_{U}^{j}) \xrightarrow{n \to \infty} Q(\alpha|\mathbf{X}_{U}^{j}). \quad (3.3)$$

*Proof.* The proof of (3.2) follows from the result in Theorem 3.3. The one of (3.3) follows from the result in Theorem 3.5 which implies, $\forall y$, $\hat{F}^{OOB}(y|\mathbf{X}_U^j) \xrightarrow{n \to \infty} F(y|\mathbf{X}_U^j)$. $\square$

## 4 QOSA index estimation

Let $U \subseteq \{1, \dots, p\}$, the first-order QOSA index estimation is obtained from (2.2) by plugin, namely $\hat{S}_U^\alpha = 1 - \hat{O}/\hat{P}$, with $\hat{O} = \frac{1}{n} \sum_{i=1}^n \psi_\alpha(Y^i, \hat{Q}_n^{OOB}(\alpha|\mathbf{X}_U^i))$ and $\hat{P} = \frac{1}{n} \sum_{i=1}^n \psi_\alpha(Y^i, \hat{q}^\alpha(Y))$.

**Assumption 5.** (Y is bounded.) $\exists C > 0$, s.t. $|Y| \le C$ almost surely.
Proof. The proof follows directly from Lemma 3.4 which provides a relation between the OOB error and the standard error and Theorem 3.3. Note that $\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j)$ and $F(y|\mathbf{X}_U^j)$ are distribution functions, so that their values are in $[0, 1]$ and the a.s. consistency from Theorem 3.3 leads to $L^2$ consistency. $\square$

### 3.3 Conditional quantile estimation

In this section, we are interested in the estimation of the conditional quantiles at fixed level $\alpha \in (0, 1)$:

$$q^{\alpha}(Y|\mathbf{X}_U) = \inf_{y \in \mathbb{R}} \{y, F_Y(y|\mathbf{X}_U) \ge \alpha\}.$$

We shall use the OOB estimator of the conditional cumulative distribution function to estimate $q^{\alpha}(Y|\mathbf{X}_U)$:

$$\hat{Q}_n^{OOB}(\alpha|\mathbf{X}_U = \mathbf{x}) = \inf_{y \in \mathbb{R}} \{y, \hat{F}^{OOB}(y|\mathbf{X}_U = \mathbf{x}) \ge \alpha\}.$$

It is easily noticed that

$$\hat{Q}_n^{OOB}(\alpha|\mathbf{X}_U = \mathbf{x}) = \inf_{k=1,...,n} \{Y^k, \hat{F}^{OOB}(Y^k|\mathbf{X}_U = \mathbf{x}) \ge \alpha\}.$$

The standard estimator $\hat{Q}_{B,s_n,n}(\alpha|\mathbf{X}_U = \mathbf{x})$ of the conditional quantile is defined in the same way, replacing $\hat{F}^{OOB}$ with $\hat{F}_{B,s_n,n}$. In the following, in order to lighten the notation, for $\mathbf{x} \in \mathcal{X}_U$ we shall write $Q(\alpha|\mathbf{x})$ for $Q^{\alpha}(Y|\mathbf{X}_U = \mathbf{x})$ and $\hat{Q}_n(\alpha|\mathbf{x})$, for an estimator of $Q^{\alpha}(Y|\mathbf{X}_U = \mathbf{x})$.

**Theorem 3.6.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then standard and OOB estimators of the conditional quantile are consistent. That is, if $\mathbf{x} \in \mathcal{X}_U$,

$$\forall \alpha \in (0, 1), \hat{Q}_{B,s_n,n}(\alpha|\mathbf{X}_U = \mathbf{x}) \xrightarrow{n \to \infty} Q(\alpha|\mathbf{X}_U = \mathbf{x}), \quad (3.2)$$

$$\forall j \in \{1, \dots, n\}, \forall \alpha \in (0, 1), \hat{Q}_n^{OOB}(\alpha|\mathbf{X}_U^j) \xrightarrow{n \to \infty} Q(\alpha|\mathbf{X}_U^j). \quad (3.3)$$

Proof. The proof of (3.2) follows from the result in Theorem 3.3. The one of (3.3) follows from the result in Theorem 3.5 which implies, $\forall y$, $\hat{F}^{OOB}(y|\mathbf{X}_U^j) \xrightarrow{n \to \infty} F(y|\mathbf{X}_U^j)$. $\square$

## 4 QOSA index estimation

Let $U \subseteq \{1, \dots, p\}$, the first-order QOSA index estimator is obtained from (2.2) by plugin, namely $\hat{S}_U^{\alpha} = 1 - \frac{\hat{O}}{\hat{P}}$, with $\hat{O} = \frac{1}{n} \sum_{i=1}^n \psi_{\alpha}(Y^i, \hat{Q}_n^{OOB}(\alpha|\mathbf{X}_U^i))$ and $\hat{P} = \frac{1}{n} \sum_{i=1}^n \psi_{\alpha}(Y^i, \hat{q}^{\alpha}(Y))$.
**Theorem 4.1.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 5 are verified, then for all $j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \lim_{n \to \infty} \mathbb{E}[(\hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^j) - Q(\alpha | \mathbf{X}_U^j))^2] \to 0. $$

*Proof.* Assumption 5 and (3.3) in Theorem 3.6 lead to the announced result

<figure><img src="image_5.png" /></figure>

**Theorem 4.2.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 5 are verified, then the plugin QOSA estimator is consistent in probability:

$$ \hat{S}_U^\alpha \xrightarrow{p} S_U^\alpha. $$

*Proof.*

$$
\begin{align*}
\hat{S}_U^\alpha &= 1 - \frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \\
&= 1 - \frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \\
&\quad + \underbrace{\frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Q(\alpha | \mathbf{X}_U^i) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}})}_{A} \\
&= 1 - \frac{1}{\hat{P}} \frac{1}{n} \underbrace{\sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}})}_{D} \\
&\quad + \underbrace{\frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}} - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}})}_{B} + \frac{1}{\hat{P}} A.
\end{align*}
$$

The law of large numbers and the almost sure consistency of $\hat{q}_\alpha(Y)$ lead to, $\hat{P} \xrightarrow{a.s.} \mathbb{E}[\psi_\alpha(Y, q^\alpha(Y))]$ and $D \xrightarrow{a.s.} \mathbb{E}[\psi_\alpha(Y, q^\alpha(Y | \mathbf{X}_U))]$. The consistency of $\hat{S}_U^\alpha$ is verified if both A and B terms go to 0 in probability. Now,

$$
\begin{align}
\mathbb{E}[|A|] &\le \mathbb{E} \left[ \frac{1}{n} \sum_{i=1}^{n} \left| (Q(\alpha | \mathbf{X}_U^i) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \right| \right] \nonumber \\
&= \mathbb{E} \left[ \left| (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1))(\alpha - \mathbb{I}_{\{Y^1 \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}}) \right| \right] \nonumber \\
&\le \mathbb{E} |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|. \tag{4.1}
\end{align}
$$
**Assumption 5.** (Y is bounded.) $\exists C > 0$, s.t. $|Y| \le C$ almost surely.

**Theorem 4.1.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 5 are verified, then for all $j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \lim_{n \to \infty} \mathbb{E}[(\hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^j) - Q(\alpha | \mathbf{X}_U^j))^2] \to 0. $$

*Proof.* Assumption 5 and (3.3) in Theorem 3.6 lead to the announced result □

**Theorem 4.2.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 5 are verified, then the plugin QOSA estimator is consistent in probability:

$$ \hat{S}_U^\alpha \xrightarrow{p} S_U^\alpha. $$

*Proof.*

$$
\begin{align*}
\hat{S}_U^\alpha &= 1 - \frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \\
&= 1 - \frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \\
&\qquad + \underbrace{\frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Q(\alpha | \mathbf{X}_U^i) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}})}_{A} \\
&= 1 - \underbrace{\frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}})}_{D} \\
&\qquad + \underbrace{\frac{1}{\hat{P}} \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}} - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) + \frac{1}{\hat{P}} A}_{B}.
\end{align*}
$$

The law of large numbers and the almost sure consistency of $\hat{q}_\alpha(Y)$ lead to, $\hat{P} \xrightarrow{a.s.} \mathbb{E}[\psi_\alpha(Y, q^\alpha(Y))]$ and $D \xrightarrow{a.s.} \mathbb{E}[\psi_\alpha(Y, q^\alpha(Y | \mathbf{X}_U))]$. The consistency of $\hat{S}_U^\alpha$ is verified if both A and B terms go to 0 in probability. Now,

$$
\begin{align}
\mathbb{E}[|A|] &\le \mathbb{E} \left[ \frac{1}{n} \sum_{i=1}^{n} \left| (Q(\alpha | \mathbf{X}_U^i) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i))(\alpha - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \right| \right] \nonumber \\
&= \mathbb{E} \left[ \left| (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1))(\alpha - \mathbb{I}_{\{Y \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}}) \right| \right] \nonumber \\
&\le \mathbb{E} |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|. \tag{4.1}
\end{align}
$$
The second line is due to equidistribution of $(Y^i, \mathbf{X}^i)$. Using Eq. (3.3) in Theorem 3.6 and Assumption 5, the upper bound in (4.1) tends to zero so that by Markov's inequality, A tends to zero in probability. We next consider the B term.

$$
\begin{align*}
\mathbb{E}[|B|] &= \mathbb{E} \left[ \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}} - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \right] \\
&\le \mathbb{E} \left[ |(Y - Q(\alpha | \mathbf{X}_U^1))(\mathbb{I}_{\{Y \le Q(\alpha | \mathbf{X}_U^1)\}} - \mathbb{I}_{\{Y \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}})| \right] \\
&\le C \mathbb{E} |\mathbb{I}_{\{Y \le Q(\alpha | \mathbf{X}_U^1)\}} - \mathbb{I}_{\{Y \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}}| = C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)) \\
&= C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1), |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| \le \varepsilon) \\
&\quad + C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1), |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| > \varepsilon) \\
&\le C \underbrace{\mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \varepsilon, Q(\alpha | \mathbf{X}_U^1) + \varepsilon))}_{T_1} + C \underbrace{\mathbb{P}(|Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| > \varepsilon)}_{T_2}
\end{align*}
$$

Under Assumption 4, $\lim_{\varepsilon \to 0} T_1 = \mathbb{P}(Y = Q(\alpha | \mathbf{X}_U)) = 0$. By Chebyshev's inequality, $T_2 \le \frac{1}{\varepsilon^2} \mathbb{E} \left[ |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|^2 \right]$. With Theorem 4.1, for any $\delta > 0$, we can find $n_0$ such that $\forall n > n_0$, $\mathbb{E} \left[ |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|^2 \right] < \delta \varepsilon^2$, which means $T_2 < \delta$. Finally $\mathbb{E}[|B|]$ goes to 0. Then B goes to 0 in probability by Markov's inequality. This achieves the proof of Theorem 4.2, that is the convergence in probability of $\hat{S}_U^\alpha$.

□

# 5 Numerical study and application

This section is devoted to numerical experiments. In order to evaluate the performance of our estimation procedure, we first consider models of the form $Y = m(\mathbf{X})$ for specific $m$ and input probability distribution $\mathbb{P}_{\mathbf{X}}$ for which the theoretical value of QOSA indices can be computed analytically. We first validate on simulated data our estimation procedure of first-order QOSA indices in Section 5.1, comparing to the state of the art. Then in Section 5.2 we focus on higher order QOSA index estimation for which, to our knowledge, no estimation procedure was proposed untill now. Finally in Section 5.3 we apply our estimation methodology to an environmental real dataset.

## 5.1 First-order QOSA index estimation

In this section, we aim at comparing our pinball estimation procedure for first-order QOSA indices with two alternative procedures from the state of the art:
The second line is due to equidistribution of $(Y^i, \mathbf{X}^i)$. Using Eq. (3.3) in Theorem 3.6 and Assumption 5, the upper bound in (4.1) tends to zero so that by Markov's inequality, A tends to zero in probability. We next consider the B term.

$$
\begin{align*}
\mathbb{E}[|B|] &= \mathbb{E} \left[ \frac{1}{n} \sum_{i=1}^{n} (Y^i - Q(\alpha | \mathbf{X}_U^i))(\mathbb{I}_{\{Y^i \le Q(\alpha | \mathbf{X}_U^i)\}} - \mathbb{I}_{\{Y^i \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^i)\}}) \right] \\
&\le \mathbb{E} \left[ |(Y - Q(\alpha | \mathbf{X}_U^1))(\mathbb{I}_{\{Y \le Q(\alpha | \mathbf{X}_U^1)\}} - \mathbb{I}_{\{Y \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}})| \right] \\
&\le C \mathbb{E} |\mathbb{I}_{\{Y \le Q(\alpha | \mathbf{X}_U^1)\}} - \mathbb{I}_{\{Y \le \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)\}}| = C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)) \\
&= C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1), |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| \le \varepsilon) \\
&\quad + C \mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1), |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| > \varepsilon) \\
&\le C \underbrace{\mathbb{P}(Y \in (Q(\alpha | \mathbf{X}_U^1) - \varepsilon, Q(\alpha | \mathbf{X}_U^1) + \varepsilon))}_{T_1} + C \underbrace{\mathbb{P}(|Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)| > \varepsilon)}_{T_2}
\end{align*}
$$

Under Assumption 4, $\lim_{\epsilon \to 0} T_1 = \mathbb{P}(Y = Q(\alpha | \mathbf{X}_U)) = 0$. By Chebyshev's inequality, $T_2 \le \frac{1}{\epsilon^2} \mathbb{E} [ |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|^2 ]$. With Theorem 4.1, for any $\delta > 0$, $\exists$ $n_0$ such that $\forall n > n_0$, $\mathbb{E} [ |Q(\alpha | \mathbf{X}_U^1) - \hat{Q}_n^{OOB}(\alpha | \mathbf{X}_U^1)|^2 ] < \delta \epsilon^2$, which means $T_2 < \delta$. Finally $\mathbb{E}[|B|]$ goes to 0. Then B goes to 0 in probability by Markov's inequality. This achieves the proof of Theorem 4.2, that is the convergence in probability of $\hat{S}_U^\alpha$. $\square$

**Comments** In this section, we have proposed a new estimator of the QOSA index $S_U^\alpha$. The main ingredients for the construction of the estimator are: (i) the generalization of the Bhat et al. (2015) splitting criterion based on the pinball loss to build an initial random forest from a learning sample; (ii) the application of the Projected Forest Algorithm introduced in Bénard et al. (2022a) to $\mathbf{X}_U$ and to the initial forest; (iii) the computation of the OOB estimator $\hat{S}_U^\alpha$. With our procedure, the estimation of the $p$ first-order QOSA indices $S_j^\alpha$ requires the learning of the initial forest in dimension $p$, with a complexity in $O(B \cdot mtry \cdot n \cdot \log(n))$, and the projection of the initial forest on each subspace spanned by $X_j$, with a total complexity in $O(p \cdot B \cdot n)$ for the $p$ projections (Bénard et al., 2022a, Appendix A). Whereas the estimation of a first-order QOSA index with the procedure introduced in Elie-Dit-Cosaque and Maume-Deschamps (2024) requires the training of a forest in dimension 1. Then the complexity to estimate all the $p$ first-order QOSA indices is in $O(p \cdot B \cdot n \cdot \log(n))$. Thus, by choosing mtry adequately (a common choice is mtry = $p/3$), we expect from the projection trick a reduction of the computational complexity. This projection trick is particularly relevant in the present work where the splitting criterion we use to learn the initial forest is based on the pinball loss, well-suited for quantile-oriented sensitivity analysis.
(i) the kernel-based method (Browne et al. (2017); Maume-Deschamps and Niang (2018)) available in the R package sensitivity, with the default setting, that is Gaussian kernel and the bandwidth selected following Wand et al. (1994);

(ii) the estimator introduced in (Elie-Dit-Cosaque and Maume-Deschamps, 2022b, Section 4.2.1) based on CART random forests constructed with only one variable, available in the Python package qosa, with the number of trees set to 100 and min_node_size optimized by cross-validation.

In the following we denote by pinball the new estimation procedure we introduced in this paper, and by Kernel, respectively CART, the alternative procedure described in (i), respectively in (ii). Concerning our approach based on projected pinball forest, for our experiments on simulated data, we fixed the number of trees to 100 and mtry parameter to p, the input space dimension. The γ parameter in Assumption (P2) was set to 0.2, and the π parameter in Assumption (P3) to 0.1. Finally, we set min_node_size = 10 for all experiments. One advantage of our approach is that results are good without any specific tuning of this last parameter. We also implemented the following additional rule to limit the number of splits: namely, we skip the split if sd$_l$/l < 0.03, where l, resp. sd$_l$, denotes the empirical mean, resp. standard deviation of (l$_i$)$_{i=1,...,n}$, with l$_i$ the loss function at each sample point. Concerning Kernel and CART, their implementation is based on a splitting of the original input-output sample D, using the first part (typically two third of the samples) to learn the conditional distribution and the second one (typically the remaining third of the samples) to estimate QOSA index by plugin. To be fair in the comparison, we implemented all the procedures with n samples (n varied from 500 to 2000). We present the comparison for the three examples. Note that in the second and third examples, Assumption 5 is violated. The results we obtain show that our procedure is robust to this assumption.

**Example 1:** As first example in this section, we consider $Y = X_1 + X_2, X_1, X_2 \sim U(0,1)$, independent. The analytical values of first-order QOSA indices are:

$$S_1^\alpha = S_2^\alpha = \begin{cases} 1 - \frac{\alpha/2 - \alpha^2/2}{\alpha - (\sqrt{2\alpha})^{3/3}} & \text{if } \alpha \ge 1/2 \\ 1 - \frac{\alpha/2 - \alpha^2/2}{\alpha - t_\alpha^2 + t_\alpha^3/3 + 1/3} & \text{if } \alpha < 1/2 \end{cases}$$

with $t_\alpha = 2 - \sqrt{2(1-\alpha)}$. In the experiments, we add a dummy variable $X_3$, that is Y is independent on $X_3$ which does not appear in the model, and test the ability of each method to detect it. We compute the root mean
# 5 Numerical study and application

This section is devoted to numerical experiments. In order to evaluate the performance of our estimation procedure, we first consider models of the form $Y = m(X)$ for specific $m$ and input probability distribution $P_X$ for which the theoretical value of QOSA indices can be computed analytically. We first validate on simulated data our estimation procedure of first-order QOSA indices in low dimension in Section 5.1, comparing to the state of the art. We also propose one rapid study of the impact of a noise on the estimation: i.e. the case where $Y = m(X) + \varepsilon$, with $\varepsilon$ a white noise. Section 5.2 is devoted to an example in higher dimension: $p = 10$ and $p = 15$. A non linear example with no analytic form of the theoretical QOSA is provided (Ishigami model) in Section 5.3. Then in Section 5.4 we focus on higher order QOSA index estimation for which, to our knowledge, no estimation procedure was proposed untill now. Finally in Section 5.5 we apply our estimation methodology to an environmental real dataset. Consistency results require the $\gamma$-regularity assumption (P2). In practice, if a node is split, the split is chosen so that each of the resulting sub-nodes contains at least a proportion $\gamma$ of the observations in the original node.

## 5.1 First-order QOSA index estimation in low dimension

In this section, we aim at comparing our pinball estimation procedure for first-order QOSA indices with two alternative procedures from the state of the art:

(i) the kernel-based method (Browne et al. (2017); Maume-Deschamps and Niang (2018)) available in the R package sensitivity, with the default setting, that is Gaussian kernel and the bandwidth selected following Wand et al. (1994);

(ii) the estimator introduced in (Elie-Dit-Cosaque and Maume-Deschamps, 2022b, Section 4.2.1) based on CART random forests constructed with only one variable, available in the Python package qosa, with the number of trees set to 100 and min_node_size optimized by cross-validation.

In the following we denote by *pinball* the new estimation procedure we introduced in this paper, and by *Kernel* and *CART* the alternative procedures described in (i) and (ii), respectively. It is important to note that the main advantage of *pinball* over *CART* is that only one forest has to be fitted and is then projected to estimate each Sobol' index of any order. We now detail the numerical settings for the different procedures.

For the *pinball* procedure, we fixed the number of trees to 100 and mtry parameter to *p*, the input space dimension. The $\gamma$ parameter in Assumption
squared error (RMSE) with respect to i-th feature as:

$$RMSE_i^\alpha = \sqrt{\frac{1}{M} \sum_{m=1}^{M} (\hat{S}_i^{\alpha,m} - S_i^\alpha)^2},$$

from M = 100 repetitions.

<figure><img src="image_6.png" /><figcaption>Figure 3: Estimation error for Example 1, at different quantile levels (0.2, 0.5, 0.8 and for different sample size (500, 1000, 2000).</figcaption></figure>

For this example with bounded inputs, the *Kernel* procedure has poor performance. In particular, it exhibits a large variance. It is well-known that kernel-based estimation procedures are prone to boundary issues. For this example, both *CART* and *pinball* show good performances. Note that the main advantage of *pinball* over *CART* is that only one forest has to be fitted and is then projected to estimate each first-order Sobol' index.
(P2) was set to 0.2, and the $\pi$ parameter in Assumption (P3) to 0.1. Finally, we set $min_node_size = 10$ for all experiments. One advantage of our approach is that results are good without any specific tuning of this last parameter. We also implemented the following additional rule to limit the number of splits: namely, we skip the split if $sd_l/\bar{l} < 0.03$, where $\bar{l}$, resp. $sd_l$, denotes the empirical mean, resp. standard deviation of $(l_i)_{i=1,...,n}$, with $l_i$ the loss function at each sample point.

The implementation of Kernel and CART procedures is based on a splitting of the original input-output sample $\mathcal{D}$, using the first part (typically two third of the samples) to learn the conditional distribution and the second one (typically the remaining third of the samples) to estimate QOSA index by plugin. To be fair in the comparison, we implemented all the procedures with $n$ samples ($n$ varied from 500 to 2000).

We present the comparison on three examples. Note that in the second and third examples, Assumption 5 is violated. The results we obtain show that our procedure is robust to this assumption.

### 5.1.1 Sum of uniform variables (Example 1)

As first example in this section, we consider $Y = X_1 + X_2, X_1, X_2 \sim \mathcal{U}(0, 1)$, independent. The analytical values of first-order QOSA indices are:

$$S_1^\alpha = S_2^\alpha = \begin{cases} 1 - \frac{\alpha/2 - \alpha^2/2}{\alpha - (\sqrt{2\alpha})^3/3} & \text{if } \alpha \ge 1/2 \\ 1 - \frac{\alpha/2 - \alpha^2/2}{\alpha - t_\alpha^2 + t_\alpha^3/3 + 1/3} & \text{if } \alpha < 1/2 \end{cases}$$

with $t_\alpha = 2 - \sqrt{2(1-\alpha)}$. In the experiments, we add a dummy variable $X_3$, that is $Y$ is independent on $X_3$ which does not appear in the model, and test the ability of each method to detect it. We compute the root mean squared error (RMSE) with respect to $i$-th feature as:

$$RMSE_i^\alpha = \sqrt{\frac{1}{M} \sum_{m=1}^{M} (\hat{S}_i^{\alpha,m} - S_i^\alpha)^2},$$

from $M = 100$ repetitions.
<table><caption>Table 1: RMSE for Example 1</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>truth</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td></tr><tr><td>500</td><td>0.038</td><td>0.038</td><td>0.027</td><td>0.234</td><td>0.234</td><td>0.115</td><td><strong>0.033</strong></td><td><strong>0.034</strong></td><td><strong>0.005</strong></td></tr><tr><td>1000</td><td><strong>0.024</strong></td><td><strong>0.028</strong></td><td>0.013</td><td>0.177</td><td>0.168</td><td>0.088</td><td>0.029</td><td>0.03</td><td><strong>0.005</strong></td></tr><tr><td rowspan="3">0.5</td><td>truth</td><td>0.25</td><td>0.25</td><td>0</td><td>0.25</td><td>0.25</td><td>0</td><td>0.25</td><td>0.25</td><td>0</td></tr><tr><td>500</td><td>0.045</td><td>0.048</td><td>0.019</td><td>0.469</td><td>0.382</td><td>0.178</td><td><strong>0.032</strong></td><td><strong>0.035</strong></td><td><strong>0.004</strong></td></tr><tr><td>1000</td><td><strong>0.028</strong></td><td>0.03</td><td>0.008</td><td>0.199</td><td>0.235</td><td>0.078</td><td><strong>0.028</strong></td><td><strong>0.027</strong></td><td><strong>0.004</strong></td></tr><tr><td rowspan="3">0.8</td><td>truth</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td></tr><tr><td>500</td><td>0.041</td><td>0.043</td><td>0.028</td><td>1.72</td><td>2.623</td><td>0.487</td><td><strong>0.033</strong></td><td><strong>0.033</strong></td><td><strong>0.005</strong></td></tr><tr><td>1000</td><td>0.027</td><td>0.027</td><td>0.014</td><td>0.451</td><td>0.411</td><td>0.212</td><td><strong>0.025</strong></td><td><strong>0.026</strong></td><td><strong>0.005</strong></td></tr><tr><td></td><td></td><td><strong>0.014</strong></td><td><strong>0.016</strong></td><td>0.007</td><td>0.286</td><td>0.219</td><td>0.11</td><td>0.024</td><td>0.025</td><td><strong>0.005</strong></td></tr></tbody></table>

**Example 2:** This is one of the examples in Fort et al. (2016). $Y = X_1 - X_2, X_i, i = 1, 2, 3$, are independent and $X_i \sim Exp(1)$. The analytical value of first-order QOSA indices for this model is:

$$S_1^\alpha = \begin{cases} \frac{(1-\alpha)(1-\log(2(1-\alpha)))+\alpha \log(\alpha)}{(1-\alpha)(1-\log(2(1-\alpha)))} & \text{if } \alpha \ge 1/2 \\ \frac{\alpha(1-\log(2\alpha))+\alpha \log(\alpha)}{\alpha(1-\log(2\alpha))} & \text{if } \alpha < 1/2 \end{cases},$$

$$S_2^\alpha = \begin{cases} \frac{(1-\alpha)(1-\log(2(1-\alpha)))+(1-\alpha)\log(1-\alpha)}{(1-\alpha)(1-\log(2(1-\alpha)))} & \text{if } \alpha \ge 1/2 \\ \frac{\alpha(1-\log(2\alpha))+(1-\alpha)\log(1-\alpha)}{\alpha(1-\log(2\alpha))} & \text{if } \alpha < 1/2 \end{cases},$$

and $X_3$ is a dummy variable.

<figure><img src="image_7.png" /><figcaption>Figure 4: Estimation error for Example 2 at different quantile levels and for different sample size. The first line is for $\alpha = 0.2$, the second line is for $\alpha = 0.5$ and the third line is for $\alpha = 0.8$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Boxplots for the stimulation error (RMSE) for Example 1, computed from $M = 100$ independent repetitions, for different quantile levels ($\alpha = 0.2, 0.5, 0.8$) and different sample sizes ($n = 500, 1000, 2000$).</figcaption></figure>

For this example with bounded inputs, the *Kernel* procedure has poor performance. In particular, it exhibits a large variance. It is well-known that kernel-based estimation procedures are prone to boundary issues. For this example, both *CART* and *pinball* show good performances. Note that the main advantage of *pinball* over *CART* is that only one forest has to be fitted and is then projected to estimate each first-order Sobol' index.

<table><caption>Table 1: RMSE for Example 1, computed from $M = 100$ independent repetitions, for different quantile levels ($\alpha = 0.2, 0.5, 0.8$) and different sample sizes ($n = 500, 1000, 2000$).</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>truth</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td></tr><tr><td>500</td><td>0.038</td><td>0.038</td><td>0.027</td><td>0.234</td><td>0.234</td><td>0.115</td><td><strong>0.033</strong></td><td><strong>0.034</strong></td><td><strong>0.005</strong></td></tr><tr><td>1000</td><td><strong>0.024</strong></td><td><strong>0.028</strong></td><td>0.013</td><td>0.177</td><td>0.168</td><td>0.088</td><td>0.029</td><td>0.03</td><td><strong>0.005</strong></td></tr><tr><td rowspan="3">0.5</td><td>truth</td><td>0.25</td><td>0.25</td><td>0</td><td>0.25</td><td>0.25</td><td>0</td><td>0.25</td><td>0.25</td><td>0</td></tr><tr><td>500</td><td>0.045</td><td>0.048</td><td>0.019</td><td>0.469</td><td>0.382</td><td>0.178</td><td><strong>0.032</strong></td><td><strong>0.035</strong></td><td><strong>0.004</strong></td></tr><tr><td>1000</td><td><strong>0.028</strong></td><td>0.03</td><td>0.008</td><td>0.199</td><td>0.235</td><td>0.078</td><td><strong>0.028</strong></td><td><strong>0.027</strong></td><td><strong>0.004</strong></td></tr><tr><td rowspan="3">0.8</td><td>truth</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td><td>0.308</td><td>0.308</td><td>0</td></tr><tr><td>500</td><td>0.041</td><td>0.043</td><td>0.028</td><td>1.72</td><td>2.623</td><td>0.487</td><td><strong>0.033</strong></td><td><strong>0.033</strong></td><td><strong>0.005</strong></td></tr><tr><td>1000</td><td>0.027</td><td>0.027</td><td>0.014</td><td>0.451</td><td>0.411</td><td>0.212</td><td><strong>0.025</strong></td><td><strong>0.026</strong></td><td><strong>0.005</strong></td></tr><tr><td></td><td>2000</td><td><strong>0.014</strong></td><td><strong>0.016</strong></td><td>0.007</td><td>0.286</td><td>0.219</td><td>0.11</td><td>0.024</td><td>0.025</td><td><strong>0.005</strong></td></tr></tbody></table>

We have also performed some more simulations for other quantile levels, using our pinball method. Below are the results for $Y = X_1 + X_2$, we do
From the boxplots presented in Fig.4, the random forest based methods (CART and pinball) have smaller variance than Kernel method. From Tab. 2, we conclude that for this example our method is the best in most cases, especially to detect the dummy variable $X_3$.

<table><caption>Table 2: RMSE for Example 2</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">0.2</td><td>truth</td><td>0.16</td><td>0.534</td><td>0</td><td>0.16</td><td>0.534</td><td>0</td><td>0.16</td><td>0.534</td><td>0</td></tr><tr><td>500</td><td>0.061</td><td><strong>0.036</strong></td><td>0.027</td><td>0.102</td><td>0.147</td><td>0.055</td><td><strong>0.029</strong></td><td>0.04</td><td><strong>0.004</strong></td></tr><tr><td>1000</td><td>0.042</td><td>0.031</td><td>0.013</td><td>0.078</td><td>0.109</td><td>0.038</td><td><strong>0.032</strong></td><td><strong>0.021</strong></td><td><strong>0.004</strong></td></tr><tr><td rowspan="3">0.5</td><td>truth</td><td>0.307</td><td>0.307</td><td>0</td><td>0.307</td><td>0.307</td><td>0</td><td>0.307</td><td>0.307</td><td>0</td></tr><tr><td>500</td><td>0.05</td><td>0.049</td><td>0.015</td><td>0.083</td><td>0.074</td><td>0.015</td><td><strong>0.032</strong></td><td><strong>0.03</strong></td><td><strong>0.002</strong></td></tr><tr><td>1000</td><td>0.033</td><td>0.034</td><td>0.006</td><td>0.063</td><td>0.065</td><td>0.008</td><td><strong>0.023</strong></td><td><strong>0.021</strong></td><td><strong>0.002</strong></td></tr><tr><td rowspan="3">0.8</td><td>truth</td><td>0.534</td><td>0.16</td><td>0</td><td>0.534</td><td>0.16</td><td>0</td><td>0.534</td><td>0.16</td><td>0</td></tr><tr><td>500</td><td><strong>0.036</strong></td><td>0.057</td><td>0.033</td><td>0.132</td><td>0.102</td><td>0.059</td><td>0.043</td><td><strong>0.03</strong></td><td><strong>0.006</strong></td></tr><tr><td>1000</td><td>0.027</td><td>0.039</td><td>0.013</td><td>0.088</td><td>0.084</td><td>0.043</td><td><strong>0.021</strong></td><td><strong>0.029</strong></td><td><strong>0.006</strong></td></tr><tr><td></td><td></td><td>2000</td><td>0.022</td><td><strong>0.024</strong></td><td>0.008</td><td>0.08</td><td>0.053</td><td>0.024</td><td><strong>0.015</strong></td><td>0.029</td><td><strong>0.006</strong></td></tr></tbody></table>

**Example 3:** This last example is one of the examples in Elie-Dit-Cosaque and Maume-Deschamps (2022a). $Y = X_1 + X_2$, with $(X_1, X_2)$ a gaussian vector, $X_1 \sim \mathcal{N}(0, 1)$, $X_2 \sim \mathcal{N}(0, 2)$ and correlation $\rho$. The analytical values of first-order QOSA indices are:

$$S_1^\alpha = 1 - \frac{\sigma_2\sqrt{1 - \rho^2}}{\sigma_Y}, \quad S_2^\alpha = 1 - \frac{\sigma_1\sqrt{1 - \rho^2}}{\sigma_Y}.$$

We also add a dummy variable $X_3$.
not include a dummy variable in these simulations; $\alpha = 0.1$ to 0.9 with a step of 0.1. We chose as a summary quality metric the square root of the mean of true value minus estimation squared, over all the variables, called mRMSE below, for 100 samples of size 2000.

<table><caption>Table 2: mRMSE for Example 1, computed from $M = 100$ independent repetitions of sample size $n = 2000$ with the pinball method, and without adding any dummy variable.</caption><thead><tr><th>α</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>mRMSE</th><td>0.026</td><td>0.024</td><td>0.026</td><td>0.027</td><td>0.026</td><td>0.025</td><td>0.025</td><td>0.025</td><td>0.023</td></tr></tbody></table>

Table 2 shows the stability of the method along different values of $\alpha$. A degradation of the estimation would be observed for more extreme quantile levels.

### 5.1.2 Difference of Exponential distributions (Example 2)

This is one of the examples in Fort et al. (2016). $Y = X_1 - X_2$, $X_i, i = 1, 2, 3$, are independent and $X_i \sim Exp(1)$. The analytical value of first-order QOSA indices for this model is:

$$S_1^\alpha = \begin{cases} \frac{(1-\alpha)(1-\log(2(1-\alpha)))+\alpha \log(\alpha)}{(1-\alpha)(1-\log(2(1-\alpha)))} & \text{if } \alpha \ge 1/2 \\ \frac{\alpha(1-\log(2\alpha))+\alpha \log(\alpha)}{\alpha(1-\log(2\alpha))} & \text{if } \alpha < 1/2 \end{cases},$$

$$S_2^\alpha = \begin{cases} \frac{(1-\alpha)(1-\log(2(1-\alpha)))+(1-\alpha)\log(1-\alpha)}{(1-\alpha)(1-\log(2(1-\alpha)))} & \text{if } \alpha \ge 1/2 \\ \frac{\alpha(1-\log(2\alpha))+(1-\alpha)\log(1-\alpha)}{\alpha(1-\log(2\alpha))} & \text{if } \alpha < 1/2 \end{cases},$$

and $X_3$ is a dummy variable.
<figure><img src="image_8.png" /><figcaption>Figure 5: Estimation error for Example 3 for different correlation coefficients, $\rho = -0.5$ (top line), $\rho = 0$ (middle line), $\rho = 0.5$ (bottom line).</figcaption></figure>

From the boxplots in Fig.5, the performance of the *CART* method decreases when dependence between variables increases. Its estimation accuracy for the dummy variable $X_3$ is poor, in comparison with both *Kernel* and *pinball* methods. The *pinball* approach we introduced reduces significantly the bias in most cases.

<table><caption>Table 3: RMSE for Example 3, scaled by $10^{-2}$</caption><thead><tr><th rowspan="2">ρ</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">-0.5</td><td>truth</td><td>0</td><td>0.5</td><td>0</td><td>0</td><td>0.5</td><td>0</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>500</td><td>13.255</td><td>8.075</td><td>12.881</td><td><strong>0.745</strong></td><td>10.939</td><td>0.833</td><td>1.662</td><td><strong>2.264</strong></td><td><strong>0.04</strong></td></tr><tr><td>1000</td><td>11.971</td><td>7.178</td><td>11.919</td><td><strong>0.434</strong></td><td>7.596</td><td>0.428</td><td>1.91</td><td><strong>1.897</strong></td><td><strong>0.043</strong></td></tr><tr><td rowspan="3">0</td><td>truth</td><td>0.106</td><td>0.553</td><td>0</td><td>0.106</td><td>0.553</td><td>0</td><td>0.106</td><td>0.553</td><td>0</td></tr><tr><td>500</td><td>4.324</td><td>3.61</td><td>1.722</td><td>3.91</td><td>10.558</td><td>0.814</td><td><strong>2.375</strong></td><td><strong>2.111</strong></td><td><strong>0.025</strong></td></tr><tr><td>1000</td><td>2.673</td><td>2.963</td><td>0.969</td><td>2.949</td><td>8.41</td><td>0.438</td><td><strong>2.244</strong></td><td><strong>1.417</strong></td><td><strong>0.028</strong></td></tr><tr><td rowspan="3">0.5</td><td>truth</td><td>0.345</td><td>0.673</td><td>0</td><td>0.345</td><td>0.673</td><td>0</td><td>0.345</td><td>0.673</td><td>0</td></tr><tr><td>500</td><td>22.271</td><td>10.304</td><td>13.071</td><td>8.669</td><td>11.101</td><td>0.956</td><td><strong>2.646</strong></td><td><strong>2.696</strong></td><td><strong>0.022</strong></td></tr><tr><td>1000</td><td>23.038</td><td>10.948</td><td>11.876</td><td>5.697</td><td>8.925</td><td>0.532</td><td><strong>2.049</strong></td><td><strong>1.16</strong></td><td><strong>0.025</strong></td></tr><tr><td></td><td>2000</td><td>23.346</td><td>11.026</td><td>11.523</td><td>5.177</td><td>6.141</td><td>0.328</td><td><strong>1.525</strong></td><td><strong>0.829</strong></td><td><strong>0.027</strong></td></tr></tbody></table>

## 5.2 Higher order QOSA index estimation

In this section, we aim to validate the estimation of Higher order QOSA indices defined in Section 2.2. Though *Kernel* and *CART* methods could be extended for this task, they are not implemented in respective packages. So we focus on results for our *pinball* projected forest method.
<figure><img src="image_5.png" /><figcaption>Figure 4: Boxplots for the estimation error (RMSE) for Example 2, computed from $M = 100$ independent repetitions, for different quantile levels, $\alpha = 0.2$ (first line), $\alpha = 0.5$ (second line) and $\alpha = 0.8$ (third line) and different sample sizes ($n = 500, 1000, 2000$).</figcaption></figure>

From the boxplots presented in Fig.4, the random forest based methods (CART and pinball) have smaller variance than Kernel method. From Tab. 3, we conclude that for this example our method is the best in most cases, especially to detect the dummy variable $X_3$.

<table><caption>Table 3: RMSE for Example 2, computed from $M = 100$ independent repetitions, for different sample sizes ($n = 500, 1000, 2000$) and different quantile levels ($\alpha = 0.2, 0.5, 0.8$).</caption><thead><tr><th rowspan="2">α</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="4">0.2</td><td>truth</td><td>0.16</td><td>0.534</td><td>0</td><td>0.16</td><td>0.534</td><td>0</td><td>0.16</td><td>0.534</td><td>0</td></tr><tr><td>500</td><td>0.061</td><td><strong>0.036</strong></td><td>0.027</td><td>0.102</td><td>0.147</td><td>0.055</td><td><strong>0.029</strong></td><td>0.04</td><td><strong>0.004</strong></td></tr><tr><td>1000</td><td>0.042</td><td>0.031</td><td>0.013</td><td>0.078</td><td>0.109</td><td>0.038</td><td><strong>0.032</strong></td><td><strong>0.021</strong></td><td><strong>0.004</strong></td></tr><tr><td>2000</td><td><strong>0.025</strong></td><td>0.02</td><td>0.007</td><td>0.051</td><td>0.082</td><td>0.026</td><td>0.029</td><td><strong>0.015</strong></td><td><strong>0.004</strong></td></tr><tr><td rowspan="4">0.5</td><td>truth</td><td>0.307</td><td>0.307</td><td>0</td><td>0.307</td><td>0.307</td><td>0</td><td>0.307</td><td>0.307</td><td>0</td></tr><tr><td>500</td><td>0.05</td><td>0.049</td><td>0.015</td><td>0.083</td><td>0.074</td><td>0.015</td><td><strong>0.032</strong></td><td><strong>0.03</strong></td><td><strong>0.002</strong></td></tr><tr><td>1000</td><td>0.033</td><td>0.034</td><td>0.006</td><td>0.063</td><td>0.065</td><td>0.008</td><td><strong>0.023</strong></td><td><strong>0.021</strong></td><td><strong>0.002</strong></td></tr><tr><td>2000</td><td>0.025</td><td>0.023</td><td>0.004</td><td>0.046</td><td>0.051</td><td>0.005</td><td><strong>0.017</strong></td><td><strong>0.02</strong></td><td><strong>0.001</strong></td></tr><tr><td rowspan="4">0.8</td><td>truth</td><td>0.534</td><td>0.16</td><td>0</td><td>0.534</td><td>0.16</td><td>0</td><td>0.534</td><td>0.16</td><td>0</td></tr><tr><td>500</td><td><strong>0.036</strong></td><td>0.057</td><td>0.033</td><td>0.132</td><td>0.102</td><td>0.059</td><td>0.043</td><td><strong>0.03</strong></td><td><strong>0.006</strong></td></tr><tr><td>1000</td><td>0.027</td><td>0.039</td><td>0.013</td><td>0.088</td><td>0.084</td><td>0.043</td><td><strong>0.021</strong></td><td><strong>0.029</strong></td><td><strong>0.006</strong></td></tr><tr><td>2000</td><td>0.022</td><td><strong>0.024</strong></td><td>0.008</td><td>0.08</td><td>0.053</td><td>0.024</td><td><strong>0.015</strong></td><td>0.029</td><td><strong>0.006</strong></td></tr></tbody></table>

We have also performed some more simulations for a noisy model, using our pinball method. More precisely, we consider $Y = X_1 - X_2 + \varepsilon$ with $X_i$ Exponential random variables, with rate parameter 1 and $\varepsilon$ a centered gaussian noise, with standard deviation $\sigma = \kappa \times \sqrt{2}$, $\sqrt{2}$ being the standard
**Example 4:** In this section, we consider a classical linear model $Y = \beta^\top X$, where $\beta = (1, 1, 1)^\top$ and $X \sim N(\mu, \Sigma)$ with $\mu = (0, 0, 0)^\top$,

$$ \Sigma = \begin{pmatrix} \sigma_1^2 & 0 & 0 \\ 0 & \sigma_2^2 & \rho\sigma_2\sigma_3 \\ 0 & \rho\sigma_2\sigma_3 & \sigma_3^2 \end{pmatrix}, -1 \le \rho \le 1, \sigma_1 = \sigma_2 = 1, \sigma_3 = 2. $$

The analytical value of order two QOSA indices for different sets $U$ are:

$$ S_{1,2}^{\alpha} = 1 - \frac{\sigma_3 \sqrt{1 - \rho^2}}{\sqrt{\beta^{\top} \Sigma \beta}}, S_{1,3}^{\alpha} = 1 - \frac{\sigma_2 \sqrt{1 - \rho^2}}{\sqrt{\beta^{\top} \Sigma \beta}}, S_{2,3}^{\alpha} = 1 - \frac{\sigma_1}{\sqrt{\beta^{\top} \Sigma \beta}}. $$

From Tab. 4, the error decreases as $n$ increases, as expected.

<table><caption>Table 4: RMSE for Example 4</caption><thead><tr><th>n</th><th>ρ</th><th>{X<sub>1</sub>, X<sub>2</sub>}</th><th>{X<sub>1</sub>, X<sub>3</sub>}</th><th>{X<sub>2</sub>, X<sub>3</sub>}</th></tr></thead><tbody><tr><td rowspan="3">500</td><td>-0.5</td><td>0.027</td><td>0.047</td><td>0.068</td></tr><tr><td>0</td><td>0.022</td><td>0.053</td><td>0.053</td></tr><tr><td>0.5</td><td>0.057</td><td>0.066</td><td>0.039</td></tr><tr><td rowspan="3">1000</td><td>-0.5</td><td>0.022</td><td>0.021</td><td>0.031</td></tr><tr><td>0</td><td>0.015</td><td>0.026</td><td>0.024</td></tr><tr><td>0.5</td><td>0.027</td><td>0.033</td><td>0.017</td></tr><tr><td rowspan="3">2000</td><td>-0.5</td><td>0.022</td><td>0.011</td><td>0.015</td></tr><tr><td>0</td><td>0.017</td><td>0.011</td><td>0.011</td></tr><tr><td>0.5</td><td>0.015</td><td>0.016</td><td>0.007</td></tr></tbody></table>

## 5.3 Application to a real dataset

In this section, we present the results obtained by applying our methodology to analyse MOCAGE data. This dataset was proposed and studied in Besse et al. (2007). The description is shown in Tab. 5. It contains 1041 observations with 10 variables (2 categorical and 8 continuous). The target is to predict O3obs (Observed ozone concentration) with the other 9 variables. The summary distribution of all continuous variables can be seen in Fig. 6.

For the implementation, we used 200 trees with $mtry = 4$, $\gamma = 0.2$ and $min_node_size = 10$. Also, there are 2 categorical variables in this dataset. The usual way to handle categorical variables is to transform them to continuous variables or transform them to dummy binary variables. We follow the dummy binary solution but with a modification. The STATION variable has 5 modalities, corresponding to 5 geographical sites, it is re-coded into 5 binary dummy variables. We bind them together is the sense that if one of these 5 variables is selected in *Mtry*, we add the left 4 other dummy variables into *Mtry*.
deviation of $X_1 - X_2$, $\kappa = 0, 0.1, 0.2$. We do not include a dummy variable and as before, we consider mRMSE as a quality metric, computed from $M = 100$ independent repetitions of size $n = 2000$.

<table><caption>Table 4: mRMSE for Example 2, computed from $M = 100$ independent repetitions, for different noise levels ($\kappa = 0, 0.1, 0.2$), different quantile levels ($\alpha = 0.1, 0.2, \dots, 0.9$) and sample size $n = 2000$, with the pinball method, and without adding any dummy variable.</caption><thead><tr><th>α</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td>κ = 0 (no noise)</td><td>0.031</td><td>0.019</td><td>0.017</td><td>0.017</td><td>0.016</td><td>0.016</td><td>0.018</td><td>0.019</td><td>0.029</td></tr><tr><td>κ = 0.1</td><td>0.043</td><td>0.027</td><td>0.022</td><td>0.018</td><td>0.016</td><td>0.019</td><td>0.021</td><td>0.029</td><td>0.042</td></tr><tr><td>κ = 0.2</td><td>0.070</td><td>0.051</td><td>0.038</td><td>0.028</td><td>0.022</td><td>0.028</td><td>0.036</td><td>0.052</td><td>0.069</td></tr></tbody></table>

As expected, Table 4 shows that the estimation error increases with the noise level.

### 5.1.3 A Gaussian example (Example 3)

This last example is one of the examples in Elie-Dit-Cosaque and Maume-Deschamps (2022a). We consider the input vector

$$X = (X_1, X_2, X_3)^\top \sim \mathcal{N} \left( (0, 0, 0)^\top, \begin{pmatrix} 1 & \sqrt{2}\rho & 0 \\ \sqrt{2}\rho & 2 & 0 \\ 0 & 0 & 1 \end{pmatrix} \right)$$

with $-1 < \rho < 1$ and the output $Y = X_1 + X_2$, so that $X_3$ is a dummy variable. Using Formulas (4.1) and (4.2) in (Elie-Dit-Cosaque and Maume-Deschamps, 2022a, Proposition 4.1), we get:

$$S_1^\alpha = 1 - \frac{\sqrt{2(1-\rho^2)}}{\sigma_Y} \quad S_2^\alpha = 1 - \frac{\sqrt{1-\rho^2}}{\sigma_Y} \quad S_3^\alpha = 0$$

$$ST_1^\alpha = \frac{\sqrt{1-\rho^2}}{\sigma_Y} \quad ST_2^\alpha = \frac{\sqrt{2(1-\rho^2)}}{\sigma_Y} \quad ST_3^\alpha = 0$$

with $\sigma_Y^2 = 3+2\sqrt{2}\rho$. The plot on Figure 5 shows the boxplots corresponding to the estimation of first-order QOSA indices in this framework. For positive integers $m$ and $n$, let $\mathbb{1}_{m,n}$ denote the matrix with entries $(\delta_{i,j})_{1\le i\le m, 1\le j\le n}$ and $\delta$ the Kronecker delta. We now we consider the input vector

$$X = (X_1, X_2, X_3, \dots, X_p)^\top \sim \mathcal{N} \left( \mathbf{0}_p, \begin{pmatrix} \left( \begin{pmatrix} 1 & \sqrt{2}\rho \\ \sqrt{2}\rho & 2 \end{pmatrix} \quad \mathbb{1}_{2,p-2} \\ \mathbb{1}_{p-2,2} & \mathbb{1}_{p-2,p-2} \end{pmatrix} \right) \right)$$

with $p > 2$ and the output $Y = X_1 + X_2$, so that $X_3, \dots, X_p$ are dummy variables. First-order and total QOSA indices of $X_1$ and $X_2$ remain unchanged, while first-order and total indices of the $p-2$ dummy variables are
The normalised QOSA index estimation at different quantile levels is shown in Fig. 7. In Besse et al. (2007); Broto et al. (2020), the variables selected as most influential are MOCAGE and TEMPE, then STATION and NO2. These two works focus on the impact of the different variables around the mean of O3obs (with GLM resp. Shapley values), while we are interested in the impact around quantiles at different levels. From Fig. 7, MOCAGE and TEMPE also are the most important variables, but the ranking is different depending on quantile levels. Also, RMH2O is the third important variable for high quantile levels (0.9, 0.8, 0.7) and STATION is more important at lowest quantile levels. This is consistent with the results in Elie-Dit-Cosque and Maume-Deschamps (2024). This example shows that QOSA indices give different informations than variance based sensitivity analysis.

<table><caption>Table 5: Summary of MOCAGE data</caption><thead><tr><th>Variable name</th><th>Type</th><th>Summary</th></tr></thead><tbody><tr><td>O3obs</td><td>Continuous</td><td>Observed ozone concentration (Response)</td></tr><tr><td>JOUR</td><td>Binary</td><td>holiday = 0, non-holiday = 1 (holiday: 724, non-holiday: 317)</td></tr><tr><td>MOCAGE</td><td>Continuous</td><td>Ozone concentration predicted by a fluid mechanics model<sup>a</sup></td></tr><tr><td>TEMPE</td><td>Continuous</td><td>Officially predicted temperatures</td></tr><tr><td>RMH2O</td><td>Continuous</td><td>Humidity ratio</td></tr><tr><td>NO2</td><td>Continuous</td><td>Nitrogen dioxide concentration</td></tr><tr><td>NO</td><td>Continuous</td><td>Nitric oxide concentration</td></tr><tr><td>STATION</td><td>Categorical</td><td>5 different sites (Aix:199, Als:222, Cad:202, Pla:208, Ram:210)<sup>b</sup></td></tr><tr><td>VentMOD</td><td>Continuous</td><td>Wind force</td></tr><tr><td>VentANG</td><td>Continuous</td><td>Wind direction</td></tr></tbody></table>

<footnote><sup>a</sup>Large Scale Atmospherical Chemestrial Model: MOCAGE (Modèle de Chimie Atmosphérique à Grande Echelle)</footnote>

<footnote><sup>b</sup>Aix=Aix-en-Provence, Ram=Rambouillet, Als=Munchhausen, Cad=Cadarache, Pla=Plan-de-Cuques</footnote>
all null. However, as the input space is of dimension $p$, it might be difficult to explore for large $p$, making estimation more difficult. Note also that for a dummy variable, only the total QOSA index is generically null. Indeed, the first-order QOSA may be different from zero if the dummy variable is correlated to other variables. For example, let us consider the input vector

$$X = (X_1, X_2, X_3)^\top \sim \mathcal{N} \left( (0, 0, 0)^\top, \begin{pmatrix} 1 & \sqrt{2}\rho & \theta \\ \sqrt{2}\rho & 2 & 0 \\ \theta & 0 & 1 \end{pmatrix} \right)$$

with $-1 < \rho < 1$ and $-\sqrt{1-\rho^2} < \theta < \sqrt{1-\rho^2}$ and the output $Y = X_1 + X_2$, so that $X_3$ is a dummy variable. Then, using once more Formulas (4.1) and (4.2) in (Elie-Dit-Cosaque and Maume-Deschamps, 2022a, Proposition 4.1), we get:

$$S_1^\alpha = 1 - \frac{\sqrt{2(1-\rho^2)}}{\sigma_Y} \qquad S_2^\alpha = 1 - \frac{\sqrt{1-\rho^2}}{\sigma_Y} \qquad S_3^\alpha = 1 - \frac{\sqrt{\sigma_Y^2 - \theta^2}}{\sigma_Y}$$

$$ST_1^\alpha = \frac{\sqrt{1-\rho^2-\theta^2}}{\sigma_Y} \qquad ST_2^\alpha = \frac{\sqrt{2(1-\frac{\rho^2}{1-\theta^2})}}{\sigma_Y} \qquad ST_3^\alpha = 0$$

with $\sigma_Y^2 = 3 + 2\sqrt{2\rho}$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Boxplots of the estimation error (RMSE) for Example 3, computed from $M = 100$ independent repetitions, for different sample sizes ($n = 500, 1000, 2000$) and different correlation coefficients ($\rho = -0.5$ in the first line, $\rho = 0$ in the second line, $\rho = 0.5$ in the third line). The indices do not depend on the quantile level for this example.</figcaption></figure>

From the boxplots in Fig.5, the performance of the CART method decreases when dependence between variables increases. Its estimation accu-
<figure><img src="image_9.png" /><figcaption>Figure 6: Summary of MOCAGE data</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: First-order QOSA indices at different quantile levels</figcaption></figure>

## References

Athey, S., Tibshirani, J., and Wager, S. (2019). Generalized random forests. *The Annals of Statistics*, 47(2):1148–1178.

Bénard, C., Biau, G., Da Veiga, S., and Scornet, E. (2022a). Shaff: Fast and consistent shapley effect estimates via random forests. In *International Conference on Artificial Intelligence and Statistics*, pages 5563–5582. PMLR.
racy for the dummy variable $X_3$ is poor, in comparison with both *Kernel* and *pinball* methods. The *pinball* approach we introduced reduces significantly the bias in most cases.

<table><caption>Table 5: RMSE for Example 3, scaled by $10^{-2}$, computed from $M = 100$ independent repetitions, for different sample sizes ($n = 500, 1000, 2000$) and different correlation coefficients ($ρ = -0.5, 0, 0.5$).</caption><thead><tr><th rowspan="2">ρ</th><th rowspan="2">n</th><th colspan="3">CART</th><th colspan="3">Kernel</th><th colspan="3">pinball</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">-0.5</td><td>truth</td><td>0</td><td>0.5</td><td>0</td><td>0</td><td>0.5</td><td>0</td><td>0</td><td>0.5</td><td>0</td></tr><tr><td>500</td><td>13.255</td><td>8.075</td><td>12.881</td><td><strong>0.745</strong></td><td>10.939</td><td>0.833</td><td>1.662</td><td><strong>2.264</strong></td><td><strong>0.04</strong></td></tr><tr><td>1000</td><td>11.971</td><td>7.178</td><td>11.919</td><td><strong>0.434</strong></td><td>7.596</td><td>0.428</td><td>1.91</td><td><strong>1.897</strong></td><td><strong>0.043</strong></td></tr><tr><td></td><td>2000</td><td>11.367</td><td>6.445</td><td>11.703</td><td><strong>0.259</strong></td><td>5.589</td><td>0.28</td><td>1.874</td><td><strong>1.351</strong></td><td><strong>0.032</strong></td></tr><tr><td rowspan="3">0</td><td>truth</td><td>0.106</td><td>0.553</td><td>0</td><td>0.106</td><td>0.553</td><td>0</td><td>0.106</td><td>0.553</td><td>0</td></tr><tr><td>500</td><td>4.324</td><td>3.61</td><td>1.722</td><td>3.91</td><td>10.558</td><td>0.814</td><td><strong>2.375</strong></td><td><strong>2.111</strong></td><td><strong>0.025</strong></td></tr><tr><td>1000</td><td>2.673</td><td>2.963</td><td>0.969</td><td>2.949</td><td>8.41</td><td>0.438</td><td><strong>2.244</strong></td><td><strong>1.417</strong></td><td><strong>0.028</strong></td></tr><tr><td></td><td>2000</td><td><strong>1.473</strong></td><td>1.744</td><td>0.458</td><td>2.166</td><td>6.757</td><td>0.281</td><td>1.822</td><td><strong>1.019</strong></td><td><strong>0.035</strong></td></tr><tr><td rowspan="3">0.5</td><td>truth</td><td>0.345</td><td>0.673</td><td>0</td><td>0.345</td><td>0.673</td><td>0</td><td>0.345</td><td>0.673</td><td>0</td></tr><tr><td>500</td><td>22.271</td><td>10.304</td><td>13.071</td><td>8.669</td><td>11.101</td><td>0.956</td><td><strong>2.646</strong></td><td><strong>2.696</strong></td><td><strong>0.022</strong></td></tr><tr><td>1000</td><td>23.038</td><td>10.948</td><td>11.876</td><td>5.697</td><td>8.925</td><td>0.532</td><td><strong>2.049</strong></td><td><strong>1.16</strong></td><td><strong>0.025</strong></td></tr><tr><td></td><td>2000</td><td>23.346</td><td>11.026</td><td>11.523</td><td>5.177</td><td>6.141</td><td>0.328</td><td><strong>1.525</strong></td><td><strong>0.829</strong></td><td><strong>0.027</strong></td></tr></tbody></table>

## 5.2 A higher dimension example

In order to explore the scalability of the method, we consider the generalized Erlang model, of the form $Y = \sum_{i=1}^{p} X_i$, with the $X_i$'s exponentially distributed with distinct rate parameters $\lambda_i$'s. The resulting output $Y$ is a generalized Erlang distribution also called Hypoexponential distribution. A semi-analytical formula for the QOSA indices is given in (Elie-Dit-Cosaque and Maume-Deschamps, 2024, Section 6.4). For $p = 10$ and $p = 15$, the $p$ values of the rate parameters $\lambda_i$ are selected regularly on the interval $[0.3, 1.25]$. Figure 6 shows the plot of theoretical QOSA values for $d = 10$.
Bénard, C., Da Veiga, S., and Scornet, E. (2022b). Mean decrease accuracy for random forests: inconsistency, and a practical solution via the sobolmda. *Biometrika*, 109(4):881-900.

Besse, P., Milhem, H., Mestre, O., Dufour, A., and Peuch, V.-H. (2007). Comparaison de techniques de « data mining » pour l'adaptation statistique des prévisions d'ozone du modèle de chimie-transport mocage. *Pollution atmosphérique*.

Bhat, H. S., Kumar, N., and Vaz, G. J. (2015). Towards scalable quantile regression trees. In *2015 IEEE International Conference on Big Data (Big Data)*, pages 53-60. IEEE.

Breiman, L. (2001). Random forests. *Machine learning*, 45:5-32.

Breiman, L., Friedman, J., Olshen, R.A., and Stone, C. (1984). *Classification and Regression Trees*. Chapman and Hall/CRC.

Broto, B., Bachoc, F., and Depecker, M. (2020). Variance reduction for estimation of shapley effects and adaptation to unknown input distribution. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):693-716.

Browne, T., Fort, J.-C., Iooss, B., and Le Gratiet, L. (2017). Estimate of quantile-oriented sensitivity indices.

Čevid, D., Michel, L., Näf, J., Meinshausen, N., and Bühlmann, P. (2022). Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. *Journal of Machine Learning Research*, 23(333):1-79.

Da Veiga, S., Gamboa, F., Iooss, B., and Prieur, C. (2021). *Basics and trends in sensitivity analysis: Theory and practice in R*. SIAM.

Efron, B. (1992). Bootstrap methods: another look at the jackknife. In *Breakthroughs in statistics: Methodology and distribution*, pages 569-593. Springer.

Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2022a). Goal-oriented shapley effects with special attention to the quantile-oriented case. *SIAM/ASA Journal on Uncertainty Quantification*, 10(3):1037-1069.

Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2022b). Random forest estimation of conditional distribution functions and conditional quantiles. *Electronic Journal of Statistics*, 16(2):6553-6583.

Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2024). Random forest based quantile-oriented sensitivity analysis indices estimation. *Computational Statistics*, pages 1-31.
<figure><img src="image_7.png" /><figcaption>Figure 6: Theoretical QOSA indices for a sum of 10 Exponential distributions, for different quantile levels $\alpha \in (0, 1)$.</figcaption></figure>

Estimations have been done for 100 samples of size 2000. As before, we have chosen the mRMSE as a quality metric.

<table><caption>Table 6: mRMSE for the generalized Erlang model in dimension $p = 10$ and $p = 15$, with the rate parameters equireparted on [0.3, 1.25]. mRMSE is computed from $M = 100$ independent repetitions for sample size $n = 2000$ and different quantile levels ($\alpha = 0.1, 0.3, 0.7, 0.9$).</caption><thead><tr><th>α</th><th>0.1</th><th>0.3</th><th>0.7</th><th>0.9</th></tr></thead><tbody><tr><th>p = 10</th><td>0.0095</td><td>0.0089</td><td>0.0091</td><td>0.0119</td></tr><tr><th>p = 15</th><td>0.0094</td><td>0.0079</td><td>0.0081</td><td>0.0099</td></tr></tbody></table>

Table 6 shows the scalability of the method to higher dimension.

## 5.3 A non linear model: the Ishigami model

All the models considered in the previous sections were linear models. The Ishigami model is widely used as a benchmark model in sensitivity analysis. In this model, input variables are $X_i$ are uniformly distributed over $[-\pi, \pi]$, $i = 1, 2, 3$, and the output is defined by

$$Y = \sin(X_1) + a \sin(X_2)^2 + b X_3^4 \sin(X_1).$$
Fort, J.-C., Klein, T., and Rachdi, N. (2016). New sensitivity analysis subordinated to a contrast. *Communications in Statistics-Theory and Methods*, 45(15):4349–4364.

Iooss, B. and Lemaître, P. (2015). A review on global sensitivity analysis methods. In *Uncertainty management in simulation-optimization of complex systems*, pages 101–122. Springer.

Kala, Z. (2019). Quantile-oriented global sensitivity analysis of design resistance. *Journal of Civil Engineering and Management*, 25(4):297–305.

Koenker, R. and Hallock, K. F. (2001). Quantile regression. *Journal of economic perspectives*, 15(4):143–156.

Lundberg, S. M. and Lee, S.-I. (2017). A unified approach to interpreting model predictions. *Advances in neural information processing systems*, 30.

Maume-Deschamps, V. and Niang, I. (2018). Estimation of quantile oriented sensitivity indices. *Statistics & Probability Letters*, 134:122–127.

Meinshausen, N. and Ridgeway, G. (2006). Quantile regression forests. *Journal of machine learning research*, 7(6).

Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., et al. (2004). *Sensitivity analysis in practice: a guide to assessing scientific models*, volume 1. Wiley Online Library.

Sobol', I. (1993). Sensitivity estimates for nonlinear mathematical models. *Math. Model. Comput. Exp.*, 1:407.

Wand, M. P., Jones, M. C., et al. (1994). Multivariate plug-in bandwidth selection. *Computational Statistics*, 9(2):97–116.

## A Technical results

This appendix contains some technical results that are needed for the proofs.

**Lemma A.1.** Assume the forest construction satisfies properties (P1-3) and that Assumption 1 is satisfied, then for any $U \subseteq \{1, \dots, p\}$, any $j \in U$, the diameter of leaves goes to 0 almost surely:

$$ \forall \mathbf{x} \in \mathcal{X}_U, \operatorname{diam}_j(A(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)) \xrightarrow[n \to \infty]{a.s.} 0, $$

$$ \text{where } \operatorname{diam}_j(A) = \sup_{\mathbf{z}, \mathbf{z}' \in A} |\mathbf{z}_j - \mathbf{z}'_j|. $$
A standard choice is $a = 7.0$, $b = 0.1$. With this standard choice, we have computed the QOSA indices for $\alpha = 0.1, \dots, 0.9$ by step of 0.1 with a Monte-Carlo method on a large sample (we call the results reference indices). More precisely, to compute the reference index $S_i^\alpha$, we used a $10^5$-sample of $(X_{i'})_{i' \neq i}$ to estimate the conditional expectations $\mathbb{E}(\psi_\alpha(Y, \theta)|X_i = X_i^{j'})$ for $X_i^j$, $1 \le j \le 2000$ a sample of $X_i$. The result is plotted in Figure 7.

<figure><img src="image_8.png" /><figcaption>Figure 7: Reference QOSA indices for the Ishigami model for quantile level $\alpha \in (0, 1)$.</figcaption></figure>

Then, we performed some simulations with this model, we computed the mRMSE (with respect to the reference indices) as a quality metric from $M = 100$ independent repetitions of size $n = 2000$.

<table><caption>Table 7: mRMSE for the Ishigami model, computed from $M = 100$ independent repetitions and sample size $n = 2000$, for different quantile levels ($\alpha = 0.1, 0.2, \dots, 0.9$).</caption><thead><tr><th>α</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td></td><td>0.0212</td><td>0.0151</td><td>0.0118</td><td>0.0110</td><td>0.0103</td><td>0.0110</td><td>0.0108</td><td>0.0142</td><td>0.0283</td></tr></tbody></table>

We see in Table 7 that the estimation errors for the Ishigami model have the same order of magnitude as the one for the linear models considered in the previous sections.
Proof. The proof follows partially the proof of Lemma 2 in Meinshausen and Ridgeway (2006) (see also Lemma 5 in Bénard et al. (2022b)). Let $j \in \{1, \dots, p\}$. We denote by $(C_k)_k$ the set of nodes from the root node to the leaf node containing $\boldsymbol{x}$, and by $S(\boldsymbol{x}, j)$ the number of such nodes for which the splitting is done along the $j$-th feature. We define the following event:

$$A = \{\exists \text{ infinitely many indices } k \text{ in } S(\boldsymbol{x}, j) \text{ for which, mtry } = 1 \text{ and } \Gamma = j\},$$

where $\Gamma$ (or $\Gamma(C_k)$ if we need to emphasize the dependency on $C_k$) is the random variable from (P3), uniform on $\{1, 2, \dots, p\}$ when mtry equals 1. It is easily seen that

$$A \subseteq \{S(\boldsymbol{x}, j) \to +\infty\}. \qquad (A.1)$$

Let mtry$(C_k)$ denote the value of mtry on the cell $C_k$. We have:

$$A = \{\forall N, \exists k \ge N \text{ with mtry}(C_k) = 1 \text{ and } \Gamma(C_k) = j\},$$

that is

$$A = \bigcap_{N} \bigcup_{k \ge N} \{\text{mtry}(C_k) = 1, \Gamma(C_k) = j\}.$$

Thus, by the Kolmogorov's zero-one law, event A has probability 0 or 1. For any fixed $N$ and $k_0 \ge N$, we have from (P3):

$$\mathbb{P}(\text{mtry}(C_{k_0}) = 1, \Gamma(C_{k_0}) = j) = \pi/p,$$

which implies $\mathbb{P}(A) \ge \pi/p$ and thus $\mathbb{P}(A) = 1$. Together with (A.1), it yields

$$S(\boldsymbol{x}, j) \xrightarrow{a.s.} \infty.$$

With the notation

$$A(\boldsymbol{x}; \mathcal{D}_n^*(\Theta), U)) = \prod_{j=1}^{p} A^{(j)}(\boldsymbol{x}; \mathcal{D}_n^*(\Theta), U))$$

and using (P1) and (P2), we get:

$$N^{(j)}(\boldsymbol{x}; \mathcal{D}_n^*(\Theta), U)) \le s_n(1 - \gamma)^{S(\boldsymbol{x}, j)}$$

with $N^{(m)}(\boldsymbol{x}; \mathcal{D}_n^*(\Theta), U))$ the number of observations whose $j$-th coordinate belongs to $A^{(j)}(\boldsymbol{x}; \mathcal{D}_n^*(\Theta), U))$. Then, from Assumption 1 and following the end of the proof of (Bénard et al., 2022b, Lemma 5), we get the result, that is $\text{diam}_j(A(\boldsymbol{x}; \mathcal{D}_n^*(\Theta))) \xrightarrow{n \to +\infty} 0$. $\square$

Then Corollary A.2 follows straightforwardly.
## 5.4 Higher order QOSA index estimation

In this section, we aim to validate the estimation of Higher order QOSA indices defined in Section 2.2. Though Kernel and CART methods could be extended for this task, they are not implemented in respective packages. So we focus on results for our pinball projected forest method.

**Example 4:** In this section, we consider a classical linear model $Y = \beta^T X$, where $\beta = (1, 1, 1)^T$ and $X \sim N(\mu, \Sigma)$ with $\mu = (0, 0, 0)^T$,

$$ \Sigma = \begin{pmatrix} \sigma_1^2 & 0 & 0 \\ 0 & \sigma_2^2 & \rho\sigma_2\sigma_3 \\ 0 & \rho\sigma_2\sigma_3 & \sigma_3^2 \end{pmatrix}, -1 \le \rho \le 1, \sigma_1 = \sigma_2 = 1, \sigma_3 = 2. $$

The analytical value of order two QOSA indices for different sets $U$ are:

$$ S_{1,2}^{\alpha} = 1 - \frac{\sigma_3 \sqrt{1 - \rho^2}}{\sqrt{\beta^T \Sigma \beta}}, S_{1,3}^{\alpha} = 1 - \frac{\sigma_2 \sqrt{1 - \rho^2}}{\sqrt{\beta^T \Sigma \beta}}, S_{2,3}^{\alpha} = 1 - \frac{\sigma_1}{\sqrt{\beta^T \Sigma \beta}}. $$

From Tab. 8, the error decreases as $n$ increases, as expected.

<table><caption>Table 8: RMSE for Example 4, computed from $M = 100$ independent repetitions, for different sample sizes ($n = 500, 1000, 2000$) and different correlation levels ($\rho = -0.5, 0, 0.5$). QOSA indices for this example do not depend on the quantile level.</caption><thead><tr><th>n</th><th>ρ</th><th>{X<sub>1</sub>, X<sub>2</sub>}</th><th>{X<sub>1</sub>, X<sub>3</sub>}</th><th>{X<sub>2</sub>, X<sub>3</sub>}</th></tr></thead><tbody><tr><td rowspan="3">500</td><td>-0.5</td><td>0.027</td><td>0.047</td><td>0.068</td></tr><tr><td>0</td><td>0.022</td><td>0.053</td><td>0.053</td></tr><tr><td>0.5</td><td>0.057</td><td>0.066</td><td>0.039</td></tr><tr><td rowspan="3">1000</td><td>-0.5</td><td>0.022</td><td>0.021</td><td>0.031</td></tr><tr><td>0</td><td>0.015</td><td>0.026</td><td>0.024</td></tr><tr><td>0.5</td><td>0.027</td><td>0.033</td><td>0.017</td></tr><tr><td rowspan="3">2000</td><td>-0.5</td><td>0.022</td><td>0.011</td><td>0.015</td></tr><tr><td>0</td><td>0.017</td><td>0.011</td><td>0.011</td></tr><tr><td>0.5</td><td>0.015</td><td>0.016</td><td>0.007</td></tr></tbody></table>

## 5.5 Application to a real dataset

In this section, we present the results obtained by applying our methodology to analyse MOCAGE data. This dataset was proposed and studied in Besse et al. (2007). The description is shown in Tab. 9. It contains 1041 observations with 10 variables (2 categorical and 8 continuous). The target is to predict O3obs (Observed ozone concentration) with the other 9 variables. The summary distribution of all continuous variables can be seen in Fig. 8.
**Corollary A.2.** Assume that Assumptions of Lemma A.1 are verified and that for any $y \in \mathbb{R}$, $F_Y(y|\mathbf{X}_U = \cdot)$ is continuous, then for each tree $b = 1, 2, ..., B$, the variation of the conditional cumulative distribution function in leaf goes to 0:

$$ \forall \mathbf{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, \quad \sup_{z \in A(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)} |F(y|\mathbf{X}_U = z) - F(y|\mathbf{X}_U = \mathbf{x})| \xrightarrow[n \to \infty]{a.s.} 0. $$

The number of elements of $\mathcal{D}_n^*(\Theta_b)$, resp. $\mathcal{D}_n^*(\Theta_b)$, in a leaf are compared in the following Lemma A.3 below.

**Lemma A.3.** For any $\varepsilon > 0$, we have

$$ \mathbb{P}(|N(\mathbf{x}; \mathcal{D}_n^*(\Theta), \mathcal{D}_{s_n}^*, U) - N(\mathbf{x}; \mathcal{D}_n^*(\Theta), \mathcal{D}_n^*(\Theta), U)| > \varepsilon) \le 24(s_n + 1)^{2p} e^{-\varepsilon^2/(188s_n)}. $$

*Proof.* The proof is similar to the one of Lemma 6.3 in Elie-Dit-Cosaque and Maume-Deschamps (2022b), where we use subsampling instead of bootstrap.

<figure><img src="image_11.png" /></figure>

We finish this appendix by sketching the proof of Lemma 3.4 which is used to obtain the consistency of the OOB estimator.

**Lemma A.4 (Lemma 4 in Bénard et al. (2022b)).** Assume the forest construction satisfies property (P1). Consider $\delta_{B,n}$ and $\gamma_{B,n}$ as follows:

$$ \delta_{B,n} = B^2 \mathbb{E} \left[ \frac{1}{|\Lambda_{n,i}|^2} \middle| 1, 2 \in \Lambda_{n,i} \right] \mathbb{P}(1, 2 \in \Lambda_{n,i}), $$

$$ \gamma_{B,n} = B^2 \mathbb{E} \left[ \frac{1}{|\Lambda_{n,i}|^2} \middle| 1 \in \Lambda_{n,i} \right] \mathbb{P}(1 \in \Lambda_{n,i}). $$

Then, for all $B \in \mathbb{N} \setminus \{0, 1\}$, we have

$$ \delta_{B,n} \le 1, \delta_{B,n} \le \gamma_{B,n} \le \frac{2}{1 - s_n/n}, $$

and for a fixed sample size $n$, $1 - \delta_{B,n} = O(B^{-1})$.

**Lemma 3.4.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $B \in \mathbb{N}^*, j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \\ \le \frac{2}{1-s_n/n} \sup_y \mathbb{E}[(\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] + O((s_n/n)^B). $$
For the implementation, we used 200 trees with mtry = 4, γ = 0.2 and min_node_size = 10. Also, there are 2 categorical variables in this dataset. The usual way to handle categorical variables is to transform them to continuous variables or transform them to dummy binary variables. We follow the dummy binary solution but with a modification. The STATION variable has 5 modalities, corresponding to 5 geographical sites, it is re-coded into 5 binary dummy variables. We bind them together is the sense that if one of these 5 variables is selected in mtry, we add the left 4 other dummy variables into mtry.

The normalised QOSA index estimation at different quantile levels is shown in Fig. 9. In Besse et al. (2007); Broto et al. (2020), the variables selected as most influential are MOCAGE and TEMPE, then STATION and NO2. These two works focus on the impact of the different variables around the mean of O3obs (with GLM resp. Shapley values), while we are interested in the impact around quantiles at different levels. From Fig. 9, MOCAGE and TEMPE also are the most important variables, but the ranking is different depending on quantile levels. Also, RMH2O is the third important variable for high quantile levels (0.9, 0.8, 0.7) and STATION is more important at lowest quantile levels. This is consistent with the results in Elie-Dit-Cosaque and Maume-Deschamps (2024). This example shows that QOSA indices give different informations than variance based sensitivity analysis.

<table><caption>Table 9: Summary of MOCAGE data</caption><thead><tr><th>Variable name</th><th>Type</th><th>Summary</th></tr></thead><tbody><tr><td>O3obs</td><td>Continuous</td><td>Observed ozone concentration (Response)</td></tr><tr><td>JOUR</td><td>Binary</td><td>holiday = 0, non-holiday = 1 (holiday: 724, non-holiday: 317)</td></tr><tr><td>MOCAGE</td><td>Continuous</td><td>Ozone concentration predicted by a fluid mechanics model<sup>a</sup></td></tr><tr><td>TEMPE</td><td>Continuous</td><td>Officially predicted temperatures</td></tr><tr><td>RMH2O</td><td>Continuous</td><td>Humidity ratio</td></tr><tr><td>NO2</td><td>Continuous</td><td>Nitrogen dioxide concentration</td></tr><tr><td>NO</td><td>Continuous</td><td>Nitric oxide concentration</td></tr><tr><td>STATION</td><td>Categorical</td><td>5 different sites (Aix:199, Als:222, Cad:202, Pla:208, Ram:210)<sup>b</sup></td></tr><tr><td>VentMOD</td><td>Continuous</td><td>Wind force</td></tr><tr><td>VentANG</td><td>Continuous</td><td>Wind direction</td></tr></tbody></table>

<footnote><sup>a</sup>Large Scale Atmospherical Chemestrial Model: MOCAGE (Modèle de Chimie Atmosphérique à Grande Echelle)</footnote>

<footnote><sup>b</sup>Aix=Aix-en-Provence, Ram=Rambouillet, Als=Munchhausen, Cad=Cadarache, Pla=Plan-de-Cuques</footnote>
Proof. We follow the proof of Lemma 2 in Bénard et al. (2022b). Consider $\mathbf{x} \in \mathcal{X}_U$. Recall the definition of standard estimator:

$$w_i(\mathbf{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}\{\mathbf{X}^i \in A_n(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)\}}{N(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)}, \quad (A.2)$$

$$\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \sum_{i=1}^{n} w_i(\mathbf{x}|U) \mathbb{I}_{(Y^i \le y)}.$$

Now we define $\hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x}; \Theta_l)$ from the following equation:

$$\hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \frac{1}{B} \sum_{l=1}^{B} \hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x}; \Theta_l).$$

Then, for fixed $y$,

$$\mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] = \mathbb{P}(|\Lambda_{n,j}| > 0)\mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j} > 0] \\ + \mathbb{P}(|\Lambda_{n,j}| = 0)\mathbb{E}[(F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j} = 0].$$

As $\mathbb{P}(|\Lambda_{n,j}| = 0) = (s_n/n)^B$, we have

$$\mathbb{P}(|\Lambda_{n,j}| = 0)\mathbb{E}[(F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j} = 0] \le \left(\frac{s_n}{n}\right)^B$$

Moreover, with the notation in Lemma A.4,

$$\mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j} > 0] \mathbb{P}(|\Lambda_{n,j}| > 0) \\ = \delta_{B,n} \mathbb{E} \left[ \frac{1}{B^2} \sum_{l,l'=1}^{B} \left( \hat{F}_{s_{n,n-1}}(y|\mathbf{X}_U^j; \Theta_l) - F(y|\mathbf{X}_U^j) \right) \left( \hat{F}_{s_{n,n-1}}(y|\mathbf{X}_U^j; \Theta_{l'}) - F(y|\mathbf{X}_U^j) \right) \right] \\ + (\gamma_{B,n} - \delta_{B,n}) \mathbb{E} \left[ \frac{1}{B^2} \sum_{l=1}^{B} \left( \hat{F}_{s_{n,n-1}}(y|\mathbf{X}_U^j; \Theta_l) - F(y|\mathbf{X}_U^j) \right)^2 \right].$$

The rest of the proof follows line by line the proof of Lemma 2 in Bénard et al. (2022b). □
<figure><img src="image_9.png" /><figcaption>Figure 8: Summary of MOCAGE data</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Normalized first-order QOSA indices for different quantile levels ($\alpha = 0.1, 0.2, \dots, 0.9$).</figcaption></figure>

Figure 9 shows that for $\alpha = 0.1, \dots, 0.5$, MOCAGE has the most important impact on quantiles, then TEMPE. The reverse holds true for $\alpha = 0.6, \dots, 0.9$. This result is coherent with what was observed in Elie-Dit-Cosaque and Maume-Deschamps (2024).
References

Athey, S., Tibshirani, J., and Wager, S. (2019). Generalized random forests. *The Annals of Statistics*, 47(2):1148–1178.

Bénard, C., Biau, G., Da Veiga, S., and Scornet, E. (2022a). Shaff: Fast and consistent shapley effect estimates via random forests. In *International Conference on Artificial Intelligence and Statistics*, pages 5563–5582. PMLR.

Bénard, C., Da Veiga, S., and Scornet, E. (2022b). Mean decrease accuracy for random forests: inconsistency, and a practical solution via the sobolmda. *Biometrika*, 109(4):881–900.

Besse, P., Milhem, H., Mestre, O., Dufour, A., and Peuch, V.-H. (2007). Comparaison de techniques de data mining pour l'adaptation statistique des prévisions d'ozone du modèle de chimie-transport mocage. *Pollution atmosphérique*.

Bhat, H. S., Kumar, N., and Vaz, G. J. (2015). Towards scalable quantile regression trees. In *2015 IEEE International Conference on Big Data (Big Data)*, pages 53–60. IEEE.

Breiman, L. (2001). Random forests. *Machine learning*, 45:5–32.

Breiman, L., Friedman, J., Olshen, R.A., and Stone, C. (1984). *Classification and Regression Trees*. Chapman and Hall/CRC.

Broto, B., Bachoc, F., and Depecker, M. (2020). Variance reduction for estimation of shapley effects and adaptation to unknown input distribution. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):693–716.

Browne, T., Fort, J.-C., Iooss, B., and Le Gratiet, L. (2017). Estimate of quantile-oriented sensitivity indices.

Ćevid, D., Michel, L., Näf, J., Meinshausen, N., and Bühlmann, P. (2022). Distributional random forests: Heterogeneity adjustment and multivariate distributional regression. *Journal of Machine Learning Research*, 23(333):1–79.

Da Veiga, S., Gamboa, F., Iooss, B., and Prieur, C. (2021). *Basics and trends in sensitivity analysis: Theory and practice in R*. SIAM.

Efron, B. (1992). Bootstrap methods: another look at the jackknife. In *Breakthroughs in statistics: Methodology and distribution*, pages 569–593. Springer.
Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2022a). Goal-oriented shapley effects with special attention to the quantile-oriented case. *SIAM/ASA Journal on Uncertainty Quantification*, 10(3):1037-1069.

Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2022b). Random forest estimation of conditional distribution functions and conditional quantiles. *Electronic Journal of Statistics*, 16(2):6553-6583.

Elie-Dit-Cosaque, K. and Maume-Deschamps, V. (2024). Random forest based quantile-oriented sensitivity analysis indices estimation. *Computational Statistics*, pages 1-31.

Fort, J.-C., Klein, T., and Rachdi, N. (2016). New sensitivity analysis subordinated to a contrast. *Communications in Statistics-Theory and Methods*, 45(15):4349-4364.

Iooss, B. and Lemaître, P. (2015). A review on global sensitivity analysis methods. In *Uncertainty management in simulation-optimization of complex systems*, pages 101-122. Springer.

Kala, Z. (2019). Quantile-oriented global sensitivity analysis of design resistance. *Journal of Civil Engineering and Management*, 25(4):297-305.

Koenker, R. and Hallock, K. F. (2001). Quantile regression. *Journal of economic perspectives*, 15(4):143-156.

Maume-Deschamps, V. and Niang, I. (2018). Estimation of quantile oriented sensitivity indices. *Statistics & Probability Letters*, 134:122-127.

Meinshausen, N. and Ridgeway, G. (2006). Quantile regression forests. *Journal of machine learning research*, 7(6).

Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., et al. (2004). *Sensitivity analysis in practice: a guide to assessing scientific models*, volume 1. Wiley Online Library.

Sobol', I. (1993). Sensitivity estimates for nonlinear mathematical models. *Math. Model. Comput. Exp.*, 1:407.

Wand, M. P., Jones, M. C., et al. (1994). Multivariate plug-in bandwidth selection. *Computational Statistics*, 9(2):97-116.

## A Technical results

This appendix contains some technical results that are needed for the proofs.
**Lemma A.1.** Assume the forest construction satisfies properties (P1-3) and that Assumption 1 is satisfied, then for any $U \subseteq \{1, \dots, p\}$, any $j \in U$, the diameter of leaves goes to 0 almost surely:

$$ \forall \mathbf{x} \in \mathcal{X}_U, \operatorname{diam}_j(A(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)) \xrightarrow[n \to \infty]{\text{a.s.}} 0, $$

$$ \text{where } \operatorname{diam}_j(A) = \sup_{\mathbf{z}, \mathbf{z}' \in A} |\mathbf{z}_j - \mathbf{z}'_j|. $$

*Proof.* The proof follows partially the proof of Lemma 2 in Meinshausen and Ridgeway (2006) (see also Lemma 5 in Bénard et al. (2022b)). Let $j \in \{1, \dots, p\}$. We denote by $(C_k)_k$ the set of nodes from the root node to the leaf node containing $\mathbf{x}$, and by $S(\mathbf{x}, j)$ the number of such nodes for which the splitting is done along the $j$-th feature. We define the following event:

$$ A = \{\exists \text{ infinitely many indices } k \text{ in } S(\mathbf{x}, j) \text{ for which mtry } = 1 \text{ and } \Gamma = j\}, $$

where $\Gamma$ (or $\Gamma(C_k)$ if we need to emphasize the dependency on $C_k$) is the random variable from (P3), uniform on $\{1, 2, \dots, p\}$ when mtry equals 1. It is easily seen that

$$ A \subseteq \{S(\mathbf{x}, j) \to +\infty\}. \qquad (A.1) $$

Let mtry$(C_k)$ denote the value of mtry on the cell $C_k$. We have:

$$ A = \{\forall N, \exists k \ge N \text{ with mtry}(C_k) = 1 \text{ and } \Gamma(C_k) = j\}, $$

that is

$$ A = \bigcap_{N} \bigcup_{k \ge N} \{\text{mtry}(C_k) = 1, \Gamma(C_k) = j\}. $$

Thus, by the Kolmogorov's zero-one law, event A has probability 0 or 1. For any fixed $N$ and $k_0 \ge N$, we have from (P3):

$$ \mathbb{P}(\text{mtry}(C_{k_0}) = 1, \Gamma(C_{k_0}) = j) = \pi/p, $$

which implies $\mathbb{P}(A) \ge \pi/p$ and thus $\mathbb{P}(A) = 1$. Together with (A.1), it yields

$$ S(\mathbf{x}, j) \xrightarrow[h \to \infty]{\text{a.s.}} \infty. $$

With the notation

$$ A(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)) = \prod_{j=1}^{p} A^{(j)}(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)) $$

and using (P1) and (P2), we get:

$$ N^{(j)}(\mathbf{x}; \mathcal{D}_n^*(\Theta), U)) \le s_n(1 - \gamma)^{S(\mathbf{x}, j)} $$
with $N^{(m)}(x; D_n^*(\Theta), U))$ the number of observations whose j-th coordinate belongs to $A^{(j)}(x; D_n^*(\Theta), U))$. Then, from Assumption 1 and following the end of the proof of (Bénard et al., 2022b, Lemma 5), we get the result, that is $\mathrm{diam}_j(A(x; D_n^*(\Theta))) \xrightarrow[n \to +\infty]{} 0$. $\square$

Then Corollary A.2 follows straightforwardly.

**Corollary A.2.** Assume that Assumptions of Lemma A.1 are verified and that for any $y \in \mathbb{R}$, $F_Y(y|\mathbf{X}_U = \cdot)$ is continuous, then for each tree $b = 1, 2, ..., B$, the variation of the conditional cumulative distribution function in leaf goes to 0:

$$ \forall \mathbf{x} \in \mathcal{X}_U, \forall y \in \mathbb{R}, \quad \sup_{z \in A(\mathbf{x}; D_n^*(\Theta), U)} |F(y|\mathbf{X}_U = z) - F(y|\mathbf{X}_U = \mathbf{x})| \xrightarrow{n \to \infty} 0. $$

The number of elements of $D_n^*(\Theta_b)$, resp. $D_n^*(\Theta_b)$, in a leaf are compared in the following Lemma A.3 below.

**Lemma A.3.** For any $\varepsilon > 0$, we have

$$ \mathbb{P}(|N(\mathbf{x}; D_n^*(\Theta), D_{s_n}^*, U) - N(\mathbf{x}; D_n^*(\Theta), D_n^*(\Theta), U)| > \varepsilon) \le 24(s_n + 1)^{2p}e^{-\varepsilon^2/(188s_n)}. $$

*Proof.* The proof is similar to the one of Lemma 6.3 in Elie-Dit-Cosaque and Maume-Deschamps (2022b), where we use subsampling instead of bootstrap. $\square$

<figure><img src="image_11.png" /></figure>

We finish this appendix by sketching the proof of Lemma 3.4 which is used to obtain the consistency of the OOB estimator.

**Lemma A.4 (Lemma 4 in Bénard et al. (2022b)).** Assume the forest construction satisfies property (P1). Consider $\delta_{B,n}$ and $\gamma_{B,n}$ as follows:

$$ \delta_{B,n} = B^2 \mathbb{E} \left[ \frac{1}{|\Lambda_{n,i}|^2} \middle| 1, 2 \in \Lambda_{n,i} \right] \mathbb{P}(1, 2 \in \Lambda_{n,i}), $$

$$ \gamma_{B,n} = B^2 \mathbb{E} \left[ \frac{1}{|\Lambda_{n,i}|^2} \middle| 1 \in \Lambda_{n,i} \right] \mathbb{P}(1 \in \Lambda_{n,i}). $$

Then, for all $B \in \mathbb{N} \setminus \{0, 1\}$, we have

$$ \delta_{B,n} \le 1, \delta_{B,n} \le \gamma_{B,n} \le \frac{2}{1 - s_n/n}, $$

and for a fixed sample size $n$, $1 - \delta_{B,n} = O(B^{-1})$.
**Lemma 3.4.** Assume that the forest construction properties (P1-3) and Assumptions 1 to 4 are verified, then for all $B \in \mathbb{N}^*, j \in \{1, 2, ..., n\}$ and $U \subseteq \{1, 2, ..., p\}$, we have

$$ \sup_y \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] \\ \leq \frac{2}{1-s_n/n} \sup_y \mathbb{E}[(\hat{F}_{B,s_n,n-1}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] + O((s_n/n)^B). $$

*Proof.* We follow the proof of Lemma 2 in Bénard et al. (2022b). Consider $\mathbf{x} \in \mathcal{X}_U$. Recall the definition of standard estimator:

$$ w_i(\mathbf{x}|U) = \frac{1}{B} \sum_{b=1}^{B} \frac{\mathbb{I}_{\{\mathbf{X}^i \in A_n(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)\}}}{N(\mathbf{x}; \mathcal{D}_n^*(b), \mathcal{D}_n^*(b), U)}, \quad (A.2) \\ \hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \sum_{i=1}^{n} w_i(\mathbf{x}|U) \mathbb{I}_{(Y^i \le y)}. $$

Now we define $\hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x}; \Theta_l)$ from the following equation:

$$ \hat{F}_{B,s_n,n}(y|\mathbf{X}_U = \mathbf{x}) = \frac{1}{B} \sum_{l=1}^{B} \hat{F}_{s_n,n}(y|\mathbf{X}_U = \mathbf{x}; \Theta_l). $$

Then, for fixed $y$,

$$ \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2] = \mathbb{P}(|\Lambda_{n,j}| > 0)\mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j}| > 0] \\ + \mathbb{P}(|\Lambda_{n,j}| = 0)\mathbb{E}[(F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j}| = 0]. $$

As $\mathbb{P}(|\Lambda_{n,j}| = 0) = (s_n/n)^B$, we have

$$ \mathbb{P}(|\Lambda_{n,j}| = 0)\mathbb{E}[(F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j}| = 0] \leq \left(\frac{s_n}{n}\right)^B $$

Moreover, with the notation in Lemma A.4,

$$ \mathbb{E}[(\hat{F}^{OOB}(y|\mathbf{X}_U^j) - F(y|\mathbf{X}_U^j))^2|\Lambda_{n,j}| > 0] \mathbb{P}(|\Lambda_{n,j}| > 0) \\ = \delta_{B,n} \mathbb{E} \left[ \frac{1}{B^2} \sum_{l,l'=1}^{B} \left( \hat{F}_{s_n,n-1}(y|\mathbf{X}_U^j; \Theta_l) - F(y|\mathbf{X}_U^j) \right) \left( \hat{F}_{s_n,n-1}(y|\mathbf{X}_U^j; \Theta_{l'}) - F(y|\mathbf{X}_U^j) \right) \right] \\ + (\gamma_{B,n} - \delta_{B,n}) \mathbb{E} \left[ \frac{1}{B^2} \sum_{l=1}^{B} \left( \hat{F}_{s_n,n-1}(y|\mathbf{X}_U^j; \Theta_l) - F(y|\mathbf{X}_U^j) \right)^2 \right]. $$

The rest of the proof follows line by line the proof of Lemma 2 in Bénard et al. (2022b). □