{'fraction_non_alphanumeric': 0.03170390730516788, 'fraction_numerical': 0.012202923962751463, 'mean_word_length': 4.3308071239105725, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Random Forest based Qantile Oriented Sensitivity Analysis indices estimation

We propose a random forest based estimation procedure for Quantile Oriented Sensitivity Analysis-QOSA. In order to be efficient, a cross validation step on the leaf size of trees is required. Our full estimation procedure is tested on both simulated data and a real dataset.


Introduction
Numerical models are ubiquitous in various fields, such as aerospace, economy, environment or insurance, they allow to approximate the behavior of physical phenomenon. Their main advantage is that they replace expensive, or even unachievable, real-life experiments and thus provide knowledge about the natural system. The extremely faithful representation of reality, made possible thanks to the increase in computing power, also explains this widespread use. However, this accuracy is often synonymous of complexity, ultimately leading to a difficult interpretation of models. Besides, model inputs are usually uncertain due to a lack of information or the random nature of factors, which means that the resulting output can be regarded as random. It is then important to assess the impact of this uncertainty on the model output. Global Sensitivity Analysis (GSA) methods solve these issues by studying how the uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model inputs (Saltelli et al., 2004). Hence, GSA allows to investigate input-ouput relationships by identifying the inputs that strongly influence the model response. Conversely, it may be of interest to see that although some inputs may not be very well established, they do not significantly contribute to output uncertainty.
Variance-based approaches are well-established and widely used for GSA. Among them, the sensitivity indices developed by Sobol (1993) are very popular. This last method stands on the assumption that the inputs are independent. Under this hypothesis, the overall variance of a scalar output can be split down into different partial variances using the so-called Hoeffding (1948) decomposition. Then, the first-order Sobol' index quantifies the individual contribution of an input to the output variance while the total Sobol' index (Jansen et al., 1994;Homma and Saltelli, 1996) measures the marginal and interaction effects. However, even if they are extremely popular and informative measures, variance-based approaches suffer from some limitation. Indeed, by definition, they study only the impact of the inputs on the expectation of the output since they consider the variance as distance measure.
A new class of sensitivity indices, generalizing the first-order Sobol' index to other quantities of interest than the expectation, has been introduced in Fort et al. (2016). These indices called Goal Oriented Sensitivity Analysis (GOSA) compare the minimum of a specific contrast function to its conditional counterpart when one of the inputs is fixed. The unconditional minimum being reached by the quantity of interest (for example a quantile).
In this paper, we focus on Quantile Oriented Sensitivity Analysis (QOSA) measuring the impact of the inputs on the α-quantile of the output distribution. Browne et al. (2017); Maume-Deschamps and Niang (2018) introduced a statistical estimator of the first-order QOSA index based on a kernel approach. Kala (2019) defined the second and higher order QOSA indices as well as a variancelike decomposition for quantiles in the case of independent inputs. Elie-Dit-Cosaque and Maume-Deschamps (2021) studied QOSA indices on various toy models in independent and dependent contexts. Despite these recents works, the question of the effective estimation of the first-order QOSA index remains open. Indeed, it turns out to be difficult to compute them in practice because it requires an accurate estimate of either the conditional quantile of the output given an input, or the minimum of a conditional expectation of the output given an input. Kala (2019) handles this feature with a brute force Monte-Carlo approach. As a matter of fact, for each value of an input, realizations of the other inputs are generated conditionally to the fixed value. Therefore, in this approach, the dependency structure of inputs has to be known, which is not always the case. Besides, the computational cost is too high to consider its use in an industrial context when dealing with costly models. Browne et al. (2017); Maume-Deschamps and Niang (2018) developed kernel-based estimators to avoid this double-loop issue. But, when using a small dataset, their performance is highly dependent of the bandwidth parameter. Browne et al. (2017) proposed a cumbersome algorithm for setting an efficient bandwidth that is not straighforward to implement in practice. As for the estimator of Maume-Deschamps and Niang (2018), a large dataset is needed in order to have a low estimation error, as no algorithm of bandwidth parameter selection is established.
To overcome these issues, we explore the random forest algorithm introduced by Breiman (2001) in order to estimate the conditional distribution of the output given an input. The main contribution of this paper is to provide different estimation strategies of the first-order QOSA index based on this method.
The paper is organized as follows. We recall in Section 2 the definition of the first-order QOSA index and initiate the estimation process. Section 3 presents the random forest algorithm and several estimators of the first-order QOSA index based on this method are described in Section 4. The entire process is summarized in Section 5. Then, the performance of the estimators is investigated in Section 6 on simulated data and the relevance of this index is highlighted on a real dataset in Section 7. Finally, a conclusion is given in Section 8.

Estimation of the QOSA index
Let us consider the input-output system where X = (X 1 , . . . , X d ) ∈ R d is a random vector of d independent inputs and Y = f (X) is the output random variable of a measurable deterministic function f : R d → R which can be a mathematical function or a computational code. Then, given a level α ∈ ]0, 1[, Fort et al. (2016) introduced the first-order Quantile Oriented Sensitivity Analysis (QOSA) index, related to the input X i , as with the contrast function ψ α : (y, θ) → (y − θ) α − 1 {y θ} . This function, also called pinball loss or check function in the literature is the cornerstone of the quantile regression (Koenker and Hallock, 2001). Quantile and conditional quantile are related to this loss function as follows Thus, the index S α i can be rewritten in the following way, Hence, as stated in Browne et al. (2017), the index S α i compares the mean distance between Y and its conditional quantile to the mean distance between Y and its quantile, where the pinball loss function ψ α is the considered distance. This index has some basic properties requested for a reasonable sensitivity index such as 0 S α It should be mentioned that Kucherenko et al. (2019) proposed new indices K α to assess the impact of inputs on the α-quantile of the output distribution. They directly quantify the mean distance between quantiles q α (Y ) and q α ( Y | X i ) rather than the mean distance between average contrast functions like in the first-order QOSA index. Different estimation strategies are investigated in their paper (brute force Monte Carlo and double-loop reordering approach). But a major limitation is that a large sample size is required to get an accurate computation of the index (samples of size 2 18 are used in their paper). Also, as mentioned in Elie-Dit-Cosaque and Maume-Deschamps (2021), the practical interpretation of the K α indices is questionable.
Let us now initiate the estimation procedure for the first-order QOSA index S α i , associated to a specific input X i and a level α. We consider an i.i.d n-sample D n = X j , Y j j=1,...,n such that Y j = f X j , j = 1, . . . , n. Then, a first natural estimator of the P term of the QOSA index based on the quantity E [ψ α (Y, q α (Y ))] is proposed with q α (Y ), the classical empirical estimator for q α (Y ) obtained from D n . The P term can be alternatively estimated as follows by using the quantity min where the minimum is reached for one of the elements of Y j j=1,...,n . As the function to minimize is decreasing then increasing, this estimator therefore requires to compute 1 . . . , n, until it increases, with Y (k) the order statistics of Y 1 , . . . , Y n . This process is much more time-consuming than the first estimator where we just need to compute the quantile and then plug it. Thus, in the sequel, we are going to use the P 1 estimator.
The O term of the QOSA index is trickier to estimate because a good approximation of the conditional distribution of Y given X i is necessary. Both existing estimators of the QOSA index currently provided in Browne et al. (2017); Maume-Deschamps and Niang (2018) handle this feature thanks to kernel-based methods. But in practice, with these methods, we are faced with determining the optimal bandwidth parameter or using large sample sizes in order to have a sufficiently low estimation error when employing a non optimal bandwidth. Thus, when dealing with costly computational models, a precise enough estimation of these indices can be difficult to achieve or even unfeasible.
We propose in this paper to address these issues by using the random forest method for estimating the conditional distribution. Therefore, several statistical estimators for the O term of the first-order QOSA index will be defined in Section 4. Let us first recall the random forest algorithm.

Random forests
Random forests are ensemble learning methods, first introduced by Breiman (2001), which can be used in classification or regression problems. We only focus on their use for regression task and assume to be given a training sample D n = X j , Y j j=1,...,n of i.i.d random variables distributed as the prototype pair (X, Y ).
Breiman's forest growns a collection of k regression trees based on the CART procedure described in Breiman et al. (1984). Building several different trees from a single dataset requires to randomize the tree building process. Randomness injected in each tree is denoted by Θ where (Θ ) =1,...,k are independent random variables distributed as Θ (independent of D n ). Θ = (Θ 1 , Θ 2 ) contains indices of observations selected to build the tree and indices of splitting candidate directions in each cell.
In more detail, the -th tree is built using a bootstrap sample D n (Θ ) from the original dataset. Only these observations are used to construct the tree and to make the tree prediction. Once the observations have been selected, the algorithm forms a recursive partitioning of the input space. In each cell, a number max f eatures of variables is selected uniformly at random among all inputs. Then, the best split is chosen as the one optimizing the CART splitting criterion only along the max f eatures preselected directions. This process is repeated in each cell. A stopping criterion, often implemented, is that a split point at any depth will only be considered if it leaves at least min samples leaf samples in each of the left and right child nodes. After tree partition has been completed, the prediction of the -th tree denoted by m b n (x; Θ , D n ) at a new point x is computed by averaging the N b n (x; Θ , D n ) observations falling into the cell A n (x; Θ , D n ) of the new point. Hence, the random forest prediction is the average of the k predicted values: By defining the random variable B j Θ 1 , D n as the number of times that the observation X j , Y j has been used from the original dataset for the -th tree construction, the conditional mean estimator in Equation (3.1) is rewritten as follows A variant of the Equation (3.2) provides another estimator of the conditional mean. Trees are still grown as in the standard random forest algorithm being based on the bootstrap samples but, for the tree prediction, the original dataset D n is used instead of the bootstrap sample D n (Θ ) associated to the -th tree and we get where the weights w o n,j (x; Θ 1 , . . . , Θ k , D n ) are defined by . (3.5) It has to be noted that contrary to Equation (3.3) where N b n (x; Θ , D n ) refers to the number of elements of D n (Θ ) falling into A n (x; Θ , D n ), in Equation (3.5), N o n (x; Θ , D n ) is the number of elements of D n that fall into A n (x; Θ , D n ).
Thus, both weighted approaches using, either the bootstrap samples (Equation (3.2)) or the original dataset (Equation (3.4)), allow to see the random forest method as a local averaging estimate (Lin and Jeon, 2006;Scornet, 2016) and will be at the heart of the strategies proposed for estimating the O term of the QOSA index. In the following, to lighten notation we will omit the dependence to Θ and D n in the weights.

Estimation of the O term of the QOSA index
By using the random forest method aforementioned, ten estimators of the O term may be defined.

The first four rely on the expression E
Since our aim is to estimate conditional expressions with respect to one input variable, say X i , we shall consider forests driven by X i , i.e. the random forest is built with the observations ..,n an independant copy of D i n .

Quantile-based O term estimators
In this section, the estimations of the O term of the QOSA index are based on the quantity ]. Using two training samples D i n and D i n , we define It is obtained thanks to two approaches based on the random forests, described in the sequel.

Quantile estimation with a weighted approach
We consider the estimator of the Conditional Cumulative Distribution Function (C_CDF) introduced in Elie-Dit-Cosaque and Maume-Deschamps (2020) using D i n to construct the forest. The C_CDF estimator used to estimate the conditional quantile is As a result, the estimator of i . Another estimator of the C_CDF can be achieved by replacing the weights w b n,j (x i ) based on the bootstrap samples of the forest by those using the original dataset w o n,j (x i ) provided in Equation (3.5). That gives the following estimator which has been proposed in Meinshausen (2006), The conditional quantiles are then estimated by plugging

Quantile estimation within a leaf
Let us consider a set of k trees indexed by = 1, . . . , k constructed with the sample D i n . For the -th tree, the estimator The values from the k randomized trees are then agregated to obtain the following random forest estimate As for the conditional mean estimate defined in Section 3 or for the C_CDF approximation introduced in Subsection 4.1.1, we can provide a variant using the original sample. Thus, once the forest is constructed with the bootstrap samples, we may estimate the conditional quantiles in the leaves of the -th tree using the original sample as follows That gives us the following random forest estimate of the conditional quantile Thus, these two methods allow us to propose the following estimator using the bootstrap samples (resp. the original sample).

Minimum-based O term estimators
The estimators developped in Subsection 4.1, based on E [ψ α (Y, q α ( Y | X i ))], require to approximate the conditional quantile and then plug it to estimate the O term. As mentioned before, the model f could be time-consuming. Therefore, they may be inappropriate as two training samples are necessary. Hence, we propose in this part to develop estimators of the O term taking advantage for which we only need to find the minimum instead of plugging the quantile.

Minimum estimation with a weighted approach
First of all, a random forest is built with the observations D i n . Then, by considering an additional sample X j j=1,...,n independent of D n , the O term may be estimated as follows

Minimum estimation within a leaf
In this subsection, we are going to take advantage of the tree structure in order to propose a new estimator. To begin with, let us consider that a random forest is built with the observations D i n . Then, the key point is that an additional sample is no longer required in order to process the outer expectation of the O term. Indeed, for the -th tree, the observations falling into its m-th leaf node denoted by A n m; Θ , D i n approximate the conditional distribution of Y given a certain point X i = x i , which allows to estimate the minimum of the conditional expectation min Then, we make the average over all the leaves of the -th tree to deal with the outer expectation. Hence, let N b n (m; Θ , D i n ) be the number of observations of the bootstrap sample D i n (Θ ) falling into the m-th leaf node and N leaves be the number of leaves in the -th tree. We define the following tree estimator for the O term The approximations of the k randomized trees are then averaged to obtain the following random forest estimate As before, the entire procedure described above is still valid using the observations of the original sample falling into the leaves of the -th tree instead of the bootstrap ones. Thanks to this change, we get the following estimator of the O term It should be noted that looking for the minimum in the leaves directly implies that they are sufficiently sampled for the method to be valid.

Minimum estimation with a weighted approach and complete trees
In Subsections 4.2.1 and 4.2.2, the conditional distribution of Y given X i is obtained from trees grown with D i n . Instead of using this approach, we propose in this part to build a forest with complete trees, i.e. grown with all the model's inputs and then adjust the weights to recover the conditional expectation Thus, as noticed, a full random forest is constructed with the whole dataset D n . Then, by using an additional sample X j j=1,...,n independent of D n , the conditional expectation (4.1) The notation X −i indicates the set of all variables except X i and we note that the conditional expectation given X i = x i is recovered by averaging over the components X −i . Thus, having independent inputs is very convenient. Otherwise, it would be necessary to know the dependency structure in order to generate the observations X l −i l=1,...,n for each new point X i = x i , which would make this estimator very cumbersome.
In addition to being used to recover the conditional expectation given X i = x i , the sample X j j=1,...,n is also used to estimate the outer expectation and we finally obtain the following estimator for the O term By using the weights w o n,j (x) instead of w b n,j (x), we may define the estimator Q 3,o i .

Overall estimation procedure
After defining the respective estimators for each term of the first-order QOSA index in Sections 2 and 4, the overall estimators are set in the following. In order to improve their accuracy, different strategies are also presented to tune hyperparameters of the random forest.

Issues with the leaf size
When using a random forest method for a regression task, a prediction is generally obtained by using the default values proposed in the packages for the max f eatures and min samples leaf hyperparameters. There are some empirical studies on the impact of these hyperparameters such as Díaz-Uriarte and De Andres (2006); Scornet (2017); Duroux and Scornet (2018) but no theoretical guarantee to support the default values.
Concerning the estimation methods of the O term of the QOSA index proposed in Section 4, except for Q 3,b i and Q 3,o i , it turns out that the values of the hyperparameters must be chosen carefully. First of all, as a forest explaining Y by X i is built for each model's input, the max f eatures hyperparameter has no impact in our procedures because it equals 1. Regarding the min samples leaf hyperparameter, its impact on the quality of the estimators is investigated through the following toy example To illustrate the influence of this hyperparameter, we present in Figure 1 the boxplot of R 1,o 1 made with 100 values for different leaf sizes. For each value of min samples leaf , an estimation R 1,o 1 is computed using two samples of size n = 10 4 and a forest grown with n trees = 500. Then, the boxplots are compared with the analytical value given below and represented with the dotted orange line on each graph in Figure 1: Based on the results obtained in Figure 1, we see that for each level α, the performance of R 1,o 1 depends highly on the choice of the min samples leaf hyperparameter. Indeed, with the grid proposed for the values of min samples leaf , the optimum value seems to be 258 for α = 0.1, 83 for α = 0.3, 47 for α = 0.7 and 27 for α = 0.9. This issue about the leaf size is only highlighted for R 1,o i but is also encountered for both methods, stated in Subsection 4.1, computing the conditional quantile with either the bootstrap samples or the original sample.
By using the same setting as in Figure 1, the distribution of Q 1,o 1 is presented in Figure 2 in order to assess the impact of the min samples leaf hyperparameter for a method where the minimum is estimated instead of plugging the quantile. The quality of Q 1,o 1 also seems to depend on the leaf size and the optimum value, allowing to well estimate E min θ∈R E [ ψ α (Y, θ)| X 1 ] for each level α, is the same as in Figure 1. As before, this concern about the leaf size was only emphasized for Q 1,o i but is also encountered for both methods, detailed in Subsections 4.2.1 and 4.2.2, approximating the minimum with either the bootstrap samples or the original sample.
For the methods Q 3,b i and Q 3,o i , based on complete trees, it seems that the tuning of the leaf size is less important as observed in Figure 3. Indeed, whatever the α level, the best results are observed for almost fully developed trees.
Thus, for all other estimators of the O term proposed in Section 4, a method giving us the optimal value of the leaf size for each level α is required to properly estimate the first-order QOSA index.

Tuning the leaf size
In order to tune the leaf size of our estimators, two methods are presented in this part. They lead to significatively improve the efficiency of the estimation. The first one rests on a classical cross-validation procedure and the second one uses the Out-Of-Bag samples.

Cross-validation procedure
The estimators of the O term developed in Subsection 4.1 are part of the conditional quantile estimation problem. Indeed, in a regression scheme, the conditional mean minimizes the expected squared error loss, while the conditional quantile q α ( Y | X i = x i ) minimizes the following expected Thus, estimators of E [ψ α (Y, q α ( Y | X i ))] established in Subsection 4.1 allow to assess the quality of the approximation of the true conditional quantile function. The smaller they are, the better the estimate of the conditional quantile function is. That is verified in Figure 1 and explains why we have this convex shape depending on the leaf size. As a matter of fact, when the value of the min samples leaf hyperparameter is incorrectly chosen, the approximation of the true conditional quantile function is wrong and so, this of Hence, in order to estimate well the conditional quantile function q α ( Y | X i ) and therefore, E [ψ α (Y, q α ( Y | X i ))] (which is our goal), the optimum value of the leaf size will be chosen within a predefined grid containing potential values as being the one minimizing the empirical generalization error computed with a K-fold cross-validation procedure. A detailed description of this process is given in Algorithm 1 with R 1,o i for instance. The principle is the same for all estimators defined in Subsection 4.1.
Algorithm 1: K-fold cross-validation procedure explained with R 1,o i Input: • Datasets: It has to be noted that the number of folds K should be chosen carefully. Indeed, a lower value of K results in a more biased estimation of the generalization error, and hence undesirable. In contrast, a larger value of K is less biased, but can suffer from large variability. The choice of K is usually 5 or 10, but there is no formal rule.

Out-Of-Bag quantile error
The estimators detailed in Subsection 4.1 deserve special attention. Indeed, another less cumbersome approach than cross-validation can be used to tune the leaf size. It is based on an adaptation to our context of the widespread "Out-Of-Bag" (OOB) error (Breiman, 1996) in regression and classification to estimate the generalization error.
• We first adapt the calculation of the OOB error for the conditional quantiles estimated with local averaging estimate of the C_CDF proposed in Subsection 4.1.1. For this purpose, we start by defining the OOB quantile error for R 1,b i . Let us fix an observation (X m i , Y m ) from D i n and consider I m as the set of trees built with the bootstrap samples not containing this observation, i.e. for which this one is "Out-Of-Bag". The conditional quantile given that where the weights are tailored to our context as follows , j = 1, . . . , n .
After this operation is carried out for all data in D i n , we calculate the error related to the approximation of the true conditional quantile function, i.e. the empirical generalization error We may use the original sample (rather than the bootstrap one) in the definition of the weights: This leads to define F o k,n as follows The estimations of the conditional quantile q o,α oob and of the OOB quantile error OOB o i follow.
• Secondly, we adapt the calculation of the OOB error for conditional quantiles estimated directly in tree leaves as introduced in Subsection 4.1.2. In that sense, define OOB quantile error for R 2,b i . Let us fix an observation (X m i , Y m ) from D i n and consider the set of trees built with the bootstrap samples not containing this observation. We then aggregate only the predictions of these trees to make our prediction After this operation carried out for all the data in D i n , we calculate the error related to the approximation of the true conditional quantile function, i.e. the empirical generalization error Again, using the original sample instead of the bootstrap one lead to define OOB o i . The advantage of these methods, compared to cross-validation techniques, is that they do not require cutting out the training sample D i n and take place during the forest construction process. Thus, given the dataset D i n and a grid containing potential values of the min samples leaf hyperparameter, a random forest is built for each one and the OOB quantile error associated is computed. Then, the optimal hyperparameter is chosen as the one with the smallest OOB error.

Full estimation procedure
Now, we have all the components in order to set the estimators of the first-order QOSA index S α i . These are separated in two classes according to the estimation method adopted for the O term. First of all, with the methods plugging the quantile, we define The whole procedure integrating the cross-validation process for these methods is detailed in Algorithm 2 (see Appendix A.1).
On the other hand, regarding the methods based on the minimum to compute the O term, we set The estimation process based on the minimum is formalized in Algorithms 3, 4 and 5. For the sake of clarity, they are all gathered in Appendix A.1 . Algorithm 3 (resp. 5) estimating the QOSA index with Q 1,b , needs a full training sample D n as well as a partial one X j j=1,...,n . While estimating the QOSA index with Q 2,b i or Q 2,o i only requires one training sample D n . This is a major advantage over methods plugging the quantile that need two full training samples.
So far, no consistency result has been proved for S α i . These various estimators are reviewed in the next section in order to establish their efficiency in practice. Moreover, all these algorithms are implemented within a python package named qosa-indices available at Elie-Dit-Cosaque (2020), it can be also freely downloaded on the PyPI website.

Numerical illustrations
Let us now carry out some simulations in order to investigate the influence of the hyperparameter optimization algorithm, the impact of the number of trees on our estimators and compare the decrease of the estimation error of each one in function of the train sample-size. From these results, the performance of the two best estimators as well as those based on kernel methods defined in Browne et al. (2017); Maume-Deschamps and Niang (2018) is assessed. Then, their scalability is tested on a toy example.

Comparison of hyperparameter optimization algorithms
We start by studying the influence of the hyperparameter optimization algorithm on the performance of our estimators plugging the conditional quantile (i.e. using This survey is carried out with the model introduced in Equation (5.1) and the following setting.
The estimators of the QOSA index are computed with samples of size n = 10 4 . The leaf size is tuned for each estimator over a grid with 20 numbers evenly spaced ranging from 5 to 300 by using either the strategy based on the OOB quantile error developed in Subsection 5.2.2 or a 3-fold crossvalidation procedure. Then, to assess the efficiency of each method (CV vs OOB), the experiment is repeated s = 100 times and the following metrics are computed with S α i , the analytical values that were provided in Fort et al. (2016). In Figure 4, for three levels α, we present the evolution of the different metrics related to the variable X 1 of our toy example in function of the number of trees ranging from 1 to 200 (in log scale). More precisely, sub-figures at the top of Figure 4 show the Root Mean Square Error (RMSE), in the middle, the bias and the variance at the bottom.
We observe that regardless of the level α and the number of trees, our estimators plugging the quantile have globally the same performance when calculated with either the OOB strategy or the cross-validation procedure. But, the run time is faster when using the OOB strategy rather than the cross-validation procedure.

Convergence with the number of trees and the train sample-size
We analyze in this part the impact of the number of trees on the performance of all our estimators except for those using Q 3,b i and Q 3,o i because of the computational cost. This survey is also carried out with the model introduced in Equation (5.1) and the following setting.  Figure 4: Evolution of RMSE, bias and variance of the estimators associated with X 1 , calculated with either the OOB strategy or the Cross-Validation procedure, in function of the number of trees for three levels α.
The estimators of the QOSA index are computed with samples of size n = 10 4 . The leaf size is tuned over a grid with 20 numbers evenly spaced ranging from 5 to 300 by using a 3-fold crossvalidation procedure for R 1,b i and R 1,o i while the strategy based on the OOB samples, developed in Subsection 5.2.2, is used for R 2,b i and R 2,o i . Regarding the minimum based estimators, the optimal leaf size is obtained via R 1,o i during the 3-fold cross-validation process. Then, the efficiency of our estimators is assessed with the metrics introduced in Equation (6.1) by repeating the experiment s = 200.
In Figure 5, for three levels α, we present the evolution of the different metrics related to the variable X 1 of our toy example in function of the number of trees ranging from 1 to 200 (in log scale). More precisely, sub-figures at the top of Figure 5 show the Root Mean Square Error (RMSE), in the middle, the bias and the variance at the bottom.
We observe that regardless of the level α, RMSE of our estimators is small. The number of trees seems to have no impact for those using Q 1,o i and Q 2,o i as the RMSE value is almost always the same. RMSE of the others decreases in function of the number of trees until it reaches a threshold starting at about 50 trees. Indeed, it is well known that from a certain number, increasing the number of trees becomes useless but results in higher calculation costs. However, we did not expect to have a stable estimation error with so few trees. Besides, still from the RMSE curves, it first appears that the estimators using the original sample (plain lines) have a lower error compared to those using the bootstrap samples (dotted lines). On  Figure 5: Evolution of RMSE, bias and variance of the estimators associated with X 1 in function of the number of trees for three levels α.
the other hand, the performance of the minimum based estimators (green and red lines) seems better than those based on the quantile (blue and orange lines). That might be explained by the additional error due to the estimation of the conditional quantile. Variance of all estimators is close to 0 and the bias curves have the same behavior as RMSE curves. This means that bias is the main/only source of error in the RMSE. This bias could be reduced by taking a larger grid where looking for the optimal leaf size during the cross-validation or using another more efficient method to find the optimum.
Let us now compare the decrease of the estimation error in function of the train sample-size. As observed in Figure 5, take a very large number of trees is not required in order to have a stable estimation error. Thus, we take n trees = 100 and the same setting as before for other parameters in the next study and observe the evolution of the metrics introduced in Equation (6.1) in function of the sample size. Figure 6 presents RMSE, bias and variance of our estimators for different sample sizes. We observe that all the metrics associated with the various estimators converge to 0 at different rates. Indeed, the convergence rates of the metrics of the quantile-based estimators are slower than those based on the minimum.
Hence, from our experiments, it turns out that the minimum-based estimators give the best results. This is an interesting feature because they need less data than those plugging the quantile. Furthermore, few trees are necessary in order to reduce the estimation error. It therefore allows to  Figure 6: Evolution of RMSE, bias and variance of the estimators associated with X 1 in function of the train sample-size for three levels α.
get a good estimation of the indices with a reasonable computational cost.

Comparison with kernel methods
In this subsection, we compare on the toy example introduced in Equation (5.1): • the kernel-based estimators proposed in Browne et al. (2017); Maume-Deschamps and Niang (2018) denoted byŠ α i and S α i , • the minimum-based QOSA index estimators building one forest for each input and using the original sample, • and the minimum-based QOSA index estimators using a forest grown with trees fully developed.
The estimators of the QOSA indices are computed with samples of size n = 10 4 . Forest methods are grown with n trees = 100. The optimal leaf size for the minimum-based estimators building one forest for each input is obtained with R 1,o i during the 3-fold cross-validation process over a grid containing 20 numbers evenly spaced ranging from 5 to 300. Regarding the minimum-based estimators using a forest grown with trees fully developed, the min samples leaf hyperparemeter equals 2. In order to have comparable methods, a cross-validation procedure is also implemented for the kernel-based estimators to choose the optimal bandwidth parameter. It is selected within over a grid containing 20 potential values ranging from 0.001 to 1. Then, we assess the performance of the different estimators by computing their empirical root mean squared error with 100 experiments.   Table 1 contains the empirical root mean square error of the different estimators associated to each input as well as the overall run time requested to obtain them. About their performance, it seems that the random forest-based estimators are better than the kernel methods. Nevertheless, as regards the methods using Q 3,b i and Q 3,o i , while they have a low error and do not need to tune the leaf size, their run time with the current implementation is too long to be used in practice. Accordingly, we recommend to compute the indices with Q 1,o i and Q 2,o i in order to get good estimations of the first-order QOSA indices in a reasonable time.

Scalability of the methods
The influence of the model's dimension d over the performance of the estimators using Q 1,o i and Q 2,o i is investigated in this subsection with the following additive exponential framework Independent inputs X i , i = 1, . . . , d, follow an Exponential distribution E(λ i ), with distinct λ i . The resulting output Y is a generalized Erlang distribution also called Hypoexponential distribution. By taking advantage of the other expression of the first-order QOSA index given in Maume-Deschamps and Niang (2018), we obtain the following semi closed-form analytical formula with Xs (−i) = j =i X j that also follows a Hypoexponential distribution. Knowing the cumulative distribution function of the Hypoexponential distribution, quantiles q α (Y ) and q α Xs (−i) are computed by numeric inversion and the analytical expression of the truncated expectations is derived from Marceau (2013).
For a specific dimension d, d values evenly spaced are selected from the interval [0.3, 1.25] and then each one represents the λ i parameter of an input X i , i = 1, . . . , d. QOSA index estimations are then computed with samples of size n = 10 4 , a forest grown with n trees = 100 and the setting defined hereafter. The leaf size is tuned with R 1,o i over a grid with 20 numbers evenly spaced ranging from 5 to 300 by using a 3-fold cross-validation. Each experiment is done 100 times in order to compute the RMSE defined in Equation (6.1) for each input, and then we take the weighted mean by the analytical values of the QOSA indices over all dimensions in order to get a global measure. Figure 7 presents the weighted RMSE as a function of the increasing dimension of our model for several levels α. For each one, we observe that the error increases slowly at the beginning until the dimension 6 for both methods then decreases. This phenomenon is due to the chosen parametrization. Indeed, when increasing the dimension of the model, the respective impact of each input is reduced. Thus, from a certain dimension, all the analytical values of the first-order QOSA indices become small and even close to 0 for some inputs. Our estimators properly capture this trend as they decrease by increasing the dimension. However, the estimator using Q 1,o i seems better than this one employing Q 1,o 2 as its error is lower.

Practical case study
We propose to apply our methodology to a practical case study. It concerns bias between the predictions from MOCAGE (Modèle de Chimie Atmosphérique à Grande Echelle) 1 and the observed ozone concentration.
This dataset has been proposed and studied in Besse et al. (2007), the full data description may be found there. It contains 10 variables with 1041 observations. O3obs: observed ozone concentration will be explained by the 9 other variables described below.
JOUR: type of day (0 for holiday vs 1 for non holiday) STATION: site of observations (5 different sites) RMH2O: humidity ratio NO2: nitrogen dioxide concentration VentMOD: wind force VentANG: wind direction NO: nitric oxide concentration TEMPE: officially predicted temperatures MOCAGE: ozone concentration predicted by a fluid mechanics model  In Besse et al. (2007); Broto et al. (2020) the central effects have been studied and lead to consider MOCAGE and TEMPE as the most influencial variables, then STATION and NO2. Our study considers various quantile level effects and shows that for quantile levels greater than 0.6, wind and RMH2O are also important, more than STATION.

Conclusion
In this paper, we introduced several estimators for the first-order QOSA index by using the random forest method. Some of them use the original sample while the others use the bootstrap samples generated during the forest construction. Both classes of estimator seem to be efficient even if we observe in our experiments that the methods using the original sample have a lower estimation error than those based on the bootstrap ones. Thus, supplementary studies should be conducted to inquire into this difference. Furthermore, the performance of these methods is highly dependent on the leaf size. This parameter could be compared to the bandwidth parameter of kernel estimators as it controls the bias of the method. But, it turns out to be easier to calibrate and we propose two methods to do this: K-fold cross-validation and Out of Bag samples based selection method.
It is also well known for random forest methods that the number of trees k should be chosen large enough to reach the desired statistical precision and small enough to make the calculations feasible as the computational cost increases linearly with k as mentioned in Scornet (2017). But, we have seen on our "toy example" that estimators proposed herein require few trees in order to have a low estimation error. This makes possible to estimate the indices correctly while maintaining a reasonable computation time.
Besides, we obtain in our application better results for our estimators when comparing with the kernel methods. A major advantage is that we have developed an estimator that requires only one training sample, whereas kernel methods require two training samples or a full one plus a partial. This feature is interesting when dealing with costly models. Another significant asset of our estimators is that their efficiency seems maintained when increasing the model dimension.
Despite these benefits, the proof for the estimators' consistency as well as the asymptotic analysis to establish the convergence rates and confidence intervals remains a major wish for the future. It is also important to remember that these indices do not have an analogue to the variance decomposition offered by Sobol indices through the theorem of Hoeffding (1948). Thus, using the values of Shapley (1953) could be interesting to get condensed and easy-to-interpret indices with a good apportionment of the interaction and dependences contributions between the inputs involved.