key: cord-0287577-xi8mw5sh
authors: Piano, Samuele Lo; Ferretti, Federico; Puy, Arnald; Albrecht, Daniel; Saltelli, Andrea
title: Variance-based sensitivity analysis: The quest for better estimators and designs between explorativity and economy
date: 2022-03-01
journal: nan
DOI: 10.1016/j.ress.2020.107300
sha: 0f46899bbbc5b36075a6dd8b57b48b2732b828f7
doc_id: 287577
cord_uid: xi8mw5sh

Variance-based sensitivity indices have established themselves as a reference among practitioners of sensitivity analysis of model outputs. A variance-based sensitivity analysis typically produces the first-order sensitivity indices $S_j$ and the so-called total-effect sensitivity indices $T_j$ for the uncertain factors of the mathematical model under analysis. The cost of the analysis depends upon the number of model evaluations needed to obtain stable and accurate values of the estimates. While efficient estimation procedures are available for $S_j$, this availability is less the case for $T_j$. When estimating these indices, one can either use a sample-based approach whose computational cost depends on the number of factors or use approaches based on meta modelling/emulators. The present work focuses on sample-based estimation procedures for $T_j$ and tests different avenues to achieve an algorithmic improvement over the existing best practices. To improve the exploration of the space of the input factors (design) and the formula to compute the indices (estimator), we propose strategies based on the concepts of economy and explorativity. We then discuss how several existing estimators perform along these characteristics. We conclude that: a) sample-based approaches based on the use of multiple matrices to enhance the economy are outperformed by designs using fewer matrices but with better explorativity; b) among the latter, asymmetric designs perform the best and outperform symmetric designs having corrective terms for spurious correlations; c) improving on the existing best practices is fraught with difficulties; and d) ameliorating the results comes at the cost of introducing extra design parameters.

The present work focuses on sample-based estimation procedures for ! for independent inputs and tests different avenues to achieve an algorithmic improvement over the existing best practices. To improve the exploration of the space of the input factors (design) and the formula to compute the indices (estimator), we propose strategies based on the concepts of economy and explorativity. We then discuss how several existing estimators perform along these characteristics.

Numerical results are presented for a set of seven test functions corresponding to different settings (few important factors with low cross-factor interactions, all factors equally important with low cross-factor interactions, and all factors equally important with high cross-factor interactions). We conclude the following from these experiments: a) sample-based approaches based on the use of multiple matrices to enhance the economy are outperformed by designs using fewer matrices but with better explorativity; b) among the latter, asymmetric designs perform the best and outperform symmetric designs having corrective terms for spurious correlations; c) improving on the existing best practices is fraught with difficulties; and d) ameliorating the results comes at the cost of introducing extra design parameters. Economy of a given design, defined as the number of elementary effects useful to compute # (defined as ) )

The sensitivity analysis of mathematical models aims to 'apportion the output uncertainty to the uncertainty in the input factors' (Saltelli and Sobol', 1995) . Uses of sensitivity analysis are found in quality assurance, model calibration, model validation, uncertainty reduction, and model simplification, which are just a few among the possible applications.

Over the last three decades, sensitivity analysis (SA) has made steps to establish itself as a self-standing discipline with a community of practitioners gathering around the SAMO (Sensitivity Analysis of Modelling Output) international conferences since 1995. Special issues have been devoted to SA (Borgonovo and Tarantola, 2012; Ginsbourger et al., 2015; Helton et al., 2006; Saltelli, 2009; Tarantola and Saint-Geours, 2015; Tarantola and Saltelli, 2003; Turányi, 2008) , mostly in relation to the SAMO events. Available textbooks for sensitivity analysis include Borgonovo (2017) , Cacuci (2003) , de Rocquigny et al. (2008) , Fang et al. (2005) , and Saltelli et al. (2008 Saltelli et al. ( , 2004 Saltelli et al. ( , 2000 . SA is acknowledged as a useful practice in model development and applications. Its use in regulatory settings (e.g., in impact assessment studies) is prescribed in guidelines both in Europe and the United States (European Commission, 2015; Office of Management and Budget, 2006; US EPA, 2015) . SA is also an ingredient of sensitivity auditing (Saltelli et al., 2013; Saltelli and Funtowicz, 2014) , a procedure to investigate the relevance and plausibility of model-based inference as an input to policy (European Commission, 2015;  Science Advice for Policy by European Academies, 2019).

Tools such as sensitivity analysis and sensitivity auditing are particularly needed at this point in time when the accuracy, relevance and plausibility of the statistical and mathematical models used to support policy are often the subject of controversy (Jakeman et al., 2006; Padilla et al., 2018; Pilkey and Pilkey-Jarvis, 2007; Saltelli and Funtowicz, 2017; Saltelli and Giampietro, 2017) , including at the time of submitting the present article, the COVID-19 pandemic (Saltelli et al., 2020; Steinmann et al., 2020) .

As highlighted elsewhere (Lo Piano and Robinson, 2019; Saltelli et al., 2019; , part of the problem in the validation of models is that the quality of the accompanying SA is often wanting. Most SA applications still favour the use of a method known as OAT, where the sensitivity of factors is gauged by moving One-factor-At-a-Time (Ferretti et al., 2016; Saltelli et al., 2019) . When a sensitivity analysis is run in this fashion, it results in a perfunctory test of the robustness of the model predictions. While different methods exist for sensitivity analysis (see recent reviews in Becker and Saltelli (2015) , Borgonovo and Plischke (2016) , Iooss and Lemaître (2015) , Neumann (2012), Norton (2015), Pianosi et al. (2016) , Saltelli et al. (2012) , and Wei et al. (2015) ), the so-called 'variance-based' methods are considered to be a reference among practitioners. To make an example, when a new method for SA is introduced, its performance is investigated against variance-based measure (see, e.g., Mara et al. (2017) ). At present, the most widely used variance-based measures are Sobol' indices (Sobol', 1993) , particularly the Sobol' first-order sensitivity measures # and the socalled total sensitivity indices # (Homma and Saltelli, 1996) . In the following, we take the suggestion from Glen and Isaacs (2012) and for simplicity adopt the symbol # , rather than )# or # ) , for the total sensitivity indices, although these notations are also commonly found in the literature.

In the next section, we briefly describe how # and # are defined and computed for the case of independent input factors (Sections 2.1-2.2). Then, we present the set of estimators used (Section 2.3) and define the concepts of economy and explorativity in the estimation procedures for # ) (Section 2.4). The experimental set up, including the test functions, is outlined in Section 3.

Section 4 is dedicated to presenting and discussing our findings, while the general conclusions on the lessons learned are drawn in Section 5.

2 Variance-based sensitivity analysis

For a scalar model output = ( , , . , … , 1 ), where , to 1 are uncertain factors, the first-order sensitivity index # can be written as

where we assume, without any loss of generality and for the case of independent variables, that the factors are uniformly distributed over the k-dimensional unit hypercube W. (1) is taken over all-factors-but-# , (written as ~# ), while the outer variance is taken over factor # . ( ) is the unconditional variance of the output variable .

A short recap of this measure should mention the following without proof (the proof of which can be found in Glen and Isaacs (2012) ).

• An efficient way to estimate # is to obtain a curve corresponding to ~! P Q # R by smoothing or regressing the scatterplot of versus the sorted values of variable # and then compute the variance of this curve over # , as shown in Figure 1 .

• The Pearson's correlation ratio squared (Pearson, 1905 (Pearson, , 1903 , the fraction of the total variability of a response that can be explained by a given set of covariates, commonly indicated as η . , coincides with # .

• When the relationship between and # is linear, # reduces to the squared value of the standardised regression coefficient . , as shown in Figure 1 .

• # is a first-order term in the variance decomposition of (valid when the input factors are independent), which includes terms up to the order , i.e.,

• Terms higher than first-order indices are used sparingly in applications due to their multiplicity: a model with = 3 has just three second-order terms, but one with = 10 has as many as forty-five second-order terms. The total number of terms in (2) is 2 1 − 1.

• The meaning of # in plain language is 'the fractional reduction in the variance of which would be obtained on average if # could be fixed'. This is derived from another useful relationship:

The second term in (3) is the average of all partial variances obtained by fixing # to a given value over its uncertainty range. Thus, the first term in (3) is the average reduction. Note that while ~! P Q # R could be greater than ( ), * ! Y~!P Q # RZ is always smaller than ( ) because of (3). The total-order sensitivity indices # can be written as follows:

The following points are worth recalling for this measure.

• Unlike the case of # , the smoothing/interpolation of the inner variance is precluded by the impossibility to sort by ~# other than by using emulators or Fourier amplitude sensitivity testing, FAST (Saltelli et al., 1999) . However, this method requires parametric equations for the search-curve exploring the input space with factor-specific frequencies. Thus, it is more labourious to set up than purely Monte Carlo methods. FAST was the most efficient strategy to compute both # and # before the work of Saltelli (2002) .

• The meaning of # , in plain language (explicitly descending from Equation (4)), is 'the fraction of variance that would remain on average if one received perfect information on all other factors ~# '.

• Applying (3) again, one gets

Noting that the second term in (5) is the first-order effect on all-but-# , one derives that the first term in (5), i.e., the numerator in the Equation (4), is the total variance of all terms in decomposition

(2) that do include factor # . For example, for a model with just three factors, one can write 1 = , + . + ? + ,. + ,? + .? + ,.?

and , = , + ,. + ,? + ,.?

• Hence, a parsimonious description of the model = ( , , . , … , 1 ) can be obtained by computing all # s and all # s. This description tells us which factors behave additively P # = # R and which do not P # < # R. For an additive model, it holds that # = # for all ,

. A limitation of this parsimonious SA is that, in the case of non-negligible interactions, it does not provide information about which factors and which order of interactions are specifically involved.

Computing the couples # , # can become cumbersome when = ( , , . , … , 1 ) is computationally time consuming. This could be the case of a mathematical model involving large systems of differential equations, a labourious optimisation programme, natural system simulators involving spatially distributed grid points, and so on. This difficulty is especially relevant for # , and this is the reason why it is the focus of our contribution.

When estimating a sensitivity measure, two elements are generally demanded: the first is a design, i.e., a strategy to arrange the sample points into the multidimensional space of the input factors; and the second is an estimator, i.e., a formula to compute the selected sensitivity measures. Different authors have suggested different designs and estimators to compute the sensitivity measures (see the contributions reviewed below). In principle, different designs could be tested for a fixed estimator and vice versa, although this is not the most common approach in the literature. In the present work, we have strived to keep the inference and conclusions relative to the design (e.g., in terms of number of sampling matrices) distinct from those relative to the estimator. Ratto et al. (2007) , Saltelli (2002) , ), Sobol' (1993 , and Sobol' et al. (2007) . Some particularly efficient algorithms for the estimation of # belong to the class of spectral methods, which may be preferred in case the model has some regularity (Prieur and Tarantola, 2016) . These include random balance designs (Tarantola et al., 2006) , discrete cosine transformations (Plischke, 2012) and an "Effective Algorithm to compute Sensitivity Indices -EASI" (Plischke, 2010) .

All of these require a total number of model evaluations that does not depend on the number of factors.

In this paper, we focus solely on # and on those estimations based on the actual evaluation of the function at the sampled points without resorting to meta-modelling approaches. We move from the recipe given in , which is in turn derived from Saltelli (2002) . A quick recap of the main ingredients of this recipe is as follows.

• The computation is based on a quasi Monte Carlo (QMC) method and makes use of quasirandom (QR) points of Sobol' LPt sequences (Sobol, 1976; Sobol', 1967) . QR sequences possess desirable uniformity properties over the unit hypercube W. QR numbers are not random:

they are designed to optimally fill the unit hypercube in the sense of avoiding inhomogeneous (clustered) points. A useful concept in this respect is that of discrepancy (Kucherenko et al., 2015) .

• Given a set of points inside W, their discrepancy is the maximum deviation of the theoretical density ( times the volume of the parallelepiped) over all possible parallelepipeds drawn within the hypercube W against the actual density (the number of points in the parallelepiped).

Sobol' LPt sequences are designed to be 'low discrepancy' and perform well in existing QR method inter-comparisons (e.g., Sobol' et al. (2011) ).

• In this instance, we use two different sequence generators for Sobol' points: Algorithm 659 (Bratley and Fox, 1988) and SobolSeq16384, distributed by Broda Ltd. (2016) and based on Sobol' et al. (2011) . The different sequences are used in the implementations of our experiment with different programming languages: Python for the former, and Matlab® for the latter.

• A relevant characteristic of LPt sequences is that its uniformity properties deteriorate moving from left to right along a row of the sequence. This means that, for any given , one would expect that the left-most columns of the sample matrix have lower discrepancy than the rightmost (Kucherenko et al., 2015) .

• The estimation of # requires points in W that are separated by what is called a 'step in the # direction'. In other words, one needs two points in W that only differ in the value of factor # .

Note that # instead requires steps in the non-# direction, e.g., couples of points where all factors but # have differing values. As discussed elsewhere (Campolongo et al., 2011) , the estimation procedure for # resembles the method of Morris. However, # is preferred to Morris since the estimation of # requires less modelling assumptions and is easier to interpret.

• Given the two × matrices and , we build an additional set of matrices that we label as (#) , where column comes from matrix and all other − 1 columns come from matrix

the &' row of a matrix whose columns come from except for column j that comes from .

• The model = ( , , . , … , 1 ) is run for all ( % ) and _ % (#)` points at a cost of ( + 1) model runs, i.e., times for the ( % ) points and times for the _ % (#)` points.

• The numerator in Equation (4), which is needed to compute # , is obtained from the following estimator (Jansen, 1999) :

• The total variance, the denominator of Equation (4), has been estimated using independent runs, i.e., those corresponding to the rows of matrix .

• Each summand in Equation (6) constitutes an elementary effect fungible for the computation of the total sensitivity index # . Sobolʹ (2001) noted that this formula (6) was originally proposed by Šaltenis and Dzemyda (1982) (in Russian), and so in the following we shall call it the Šaltenis estimator.

In this contribution, we compare the Šaltenis estimator with Saltelli's design (Saltelli 2002) to those of Glen and Isaacs (2012) , Owen (Iooss et al., 2020) , and Lamboni (2018) 

The estimators used in Glen and Isaacs (2012) are symmetric (the two base matrices and are entrusted the same role and importance) and based on computing the Pearson correlation coefficients between vectors (not to be confused with the Pearson correlation ratio discussed in Section 2.1). This means that for each of the couples of vectors just described (see Table 2 ), one first computes the correlation coefficients. For example, for the first entry in Table 2 , instead of applying the Šaltenis estimator (6), one calculates the following:

where 〈 ( % )〉 is the mean of the ( % )s over the runs and P ( % )R is their variance. This is also the case for _ % (#)`. For simplicity, we have not indicated the dependence of r # upon the selected couples of function values ( % ) and _ % (#)`. The best performing estimator according to Glen and Isaacs (2012) -named D3 in their manuscript -has also been used in this study (Equation 8 ), where the

and the other terms are detailed in the Appendix (Table A1) .

Additionally, Glen and Isaacs (2012) note that supposedly uncorrelated vectors, such as ( % ) and

), may be affected by spurious correlations for finite values of . We say 'supposedly uncorrelated' since no columns are shared between ( % ) and ( % ), nor is this the case for _ % (#)` and _ % (#)`. These spurious correlations are explicitly computed in Glen and Isaacs (2012) and then used as correction terms in the computation of the sensitivity indices.

The main advantage of the symmetric design proposed over the asymmetric design of is that the coordinates of the base sample appear disproportionately with respect to the other coordinates in Saltelli (2002) while this is not the case with Glen and Isaacs (2012) .

To assess whether the use of multiple matrices can be beneficial, we tested the Owen estimator (Iooss et al., 2020) as Equation (9) 

). According to the author (Owen, 2013) , three-base matrices estimators offer better accuracy than two-base matrices ones.

The same case is made by Lamboni (2018) as regards the use of estimators having even more base matrices as per Equation (10) a~

where S,% is a generic row of a base matrix and S-P,% (#) the same row but for coordinate , as well as the asymmetric design of the Šaltenis estimator for a variable number of matrices.

As discussed in , there is a natural trade-off between the economy of a design (how many elementary effects we can obtain with a given number of runs, as in Equation 11) and how well the design fills or explores the input-factor space, i.e., in our case, the unit hypercube W. These aspects are taken here as measures of the quality of different estimators. In this contribution, economy ( ) and explorativity ( ) are defined in the context of the calculation of ! .

The trade-off between the economy ( ) and explorativity ( ) of the hypercube comes from the fact that one should strive to use any given point more than once to make a computation efficient. However, the more a given point is re-used, the less the k-dimensional space of the input is explored. In practice, when using matrices , , , etc. each of column length N, as well as hybrid matrices such as We know from the previous subsection that the estimator described in Equation (6) yields elementary effects per factor, i.e., differences ( % ) − _ % (#)` are used in Equation (6) at the cost of ( + 1) runs of the model. The economy e of the design is thus the following:

which is less than one.

From now on we shall call this design 'asymmetric' due to the different roles entrusted to matrices and : the coordinates of are used more than those of . 3 (1 + 2 ) runs has been generated. Each of the three matrices , and can be used to compute 2 effects, as shown in the first six rows of Table 2 below. 3 additional effects can be obtained by mixing the hybrid matrices (last three rows in Table 2 ). 

This gives a total of 9 effects for the case of = 3 matrices for an economy = 9 3 (1 + 2 ) 1 = 3 (1 + 2 ) 1 .

How can this be extended to a design with a generic number of matrices? Given matrices, there are _ 2` pairwise combinations; and for each of the 2 matrices that are produced, there is twice the number of factors (since for each couple of matrices, such as and , we shall have to consider both matrices (#) and (#) ). Since each matrix is composed of runs, ) will be Y + 2 _ 2`Z = P + ( − 1)R = P1 + ( − 1)R. With similar considerations, one derives that the number of

. ( − 1).

In summary, we have

and the resulting economy is defined in Eq. (14), whereby the value of e tends to 2 1 for a large or/and a large . 

Note that the same development made for # could be replicated for # , although first-order indices are not in the scope of this manuscript.

Different arrangements can be explored to calculate # to have the couples of points differing for the &' coordinate only. These settings can be compared against how many coordinates are used out of the maximum number of ) coordinates available (explorativity, ). The lower is, the less explorative the design, although the design's economy ( ) may be increased in these settings. Possible arrangements are detailed as follows. We recall from our legend ( Table 1) coordinates. Thus, the fraction of coordinates relative to the maximum ! is = " " ($%&)/)

which tends to ½ as k increases.

1T, points are initially generated in the hypercube, which is the core of the star, from which each of the available k dimensions is explored in turn. In this way, each star is made of k+1 points and needs k+k coordinates: k for the centre point of the star and one for each of its k rays. Thus, the fraction is now the follows:

which decreases as k increases.

In a winding-stair design, the hypercube is explored using a curve whereby each coordinate is increased in turn. This design needs k coordinates to generate the first point and ! − 1 additional coordinates for the remaining points 

which generally tends to 1/k as N>>k. This changes if one uses more than one trajectory. For example, if one uses trajectories of length k+1, the explorativity becomes identical to that of 'stars' above. (,) , (.) . .

(1) with column length N. This corresponds to a total of ! = ( + 1) points for a total of ( + 1)) coordinates. However, one only needs a total of 2Nk coordinates, Nk for each of the two matrices A and B. Hence,

1T,

which is the same as the 'stars' above. This design is one of the most widely used standards in the literature against which we will be benchmark the performance of the assessed estimators.

In this design, one makes use of both sets (,) , (.) . .

(1) and

(,) , (.) . .

(1) for a total of ! = 2 ( + 1) points, which in principle would correspond to a total of ! = 2 ( + 1) coordinates. However, only 2Nk coordinates are needed, Nk for each of the two matrices A and B. Hence,

Glen and Isaacs (2012) (estimator D3, symmetric). This is the same as Saltelli 2002, symmetric above.

). Only the matrices , , 

Generalisation of the symmetric case. The design now includes n base matrices A, B, …X, where X is the n th matrix, plus a total of two times _ 2` additional hybrid matrices of the type (,) , (.) . . 

which decreases as n increases and reduces to (15) for n=2. The overall cost of the sensitivity analysis in terms of ) and the number of factors are typically known to the modeller prior to the analysis. Hence, different couples of and can be chosen to meet the target ) value. To maximise ) , one would set as high as possible. In terms of discrepancy ( ) and , one would rather have lower values of to have less repeated coordinates. The lower the fraction of repeated coordinates (the higher ) is, the better the space-filling properties of the design, and hence is lower. Therefore, and have an inverse relation.

Note that the relation between the discrepancy and error is not simple. A given can be perfect for a smooth function and inadequate for a jigsaw-shaped one. With fixed and ) , there is an inverse quadratic relation between and , as shown in Equation (12), which describes the trade-off between and .

Let us examine the case where = 6 and the cost ) ≈ 500 runs (Table 3) . Since needs to be a power of two in quasi-random number sequences based on Sobol' LPt sequences, the value of ) may deviate from 500. has been estimated using the computational method provided by Jäckel (2002) and rounded to two digits ( Table 3 ).

The last row has the highest number of effects using the smallest number of random points -, i.e., just one row of a single QMC matrix in six dimensions. The opposite applies to the first row since it uses as many as 64 rows from two QMC matrices, but it gives the fewest effects. In other words, the first row is the most explorative while the last is the most economical in terms of the number of effects per run.

Using > 2 resulted in poorer convergence with respect to the case = 2. According to that paper, the Šaltenis estimator in conjunction with = 2 was the best available sample-based practice. Contrasting findings have been reported by Lamboni (2018) , according to whom the optimal number of matrices may be different from two depending on the function evaluated.

We also tested whether a variable explorativity across factors could improve the accuracy of the estimate. This experiment consisted of allocating more runs to the factors having the highest standard deviation of the elementary effects after an initial warm up. The number of residual runs is attributed according to the importance obtained at that given point in the simulation. The computational details of this experiment are described in section 4. and Glen and Isaacs (2012) base their analysis on a single function, namely, the widely used G function, which is defined as follows:

With this test function, one can modulate the importance of a factor # via the associated constant # , as shown in Table 4 . Although this function can be attributed to Davis and Rabinowitz (1984) , it was further developed by Saltelli and Sobol' (1995) and is known among practitioners as Sobol's function. It is reduced to the function used in Davis and Rabinowitz (1984) when all # = 0.

A six-dimensional version of the function with coefficients # = {0 0.5 3 9 99 99} is used here as in Glen and Isaacs (2012) , , and function A2 below (Equation 25). To test the effectiveness of the estimators with a wider typology of functions, we have used the taxonomy suggested by Kucherenko et al. (2011) in Equations (24 -30) .

The analytical values for the sensitivity indices are available from the GitHub repository: https://github.com/Confareneoclassico/New_estimator_algorithm. When the coefficients aj of function G are all equal and large, one is dealing with a B-type function (B3, Equation 28, for which all # =6.42). By contrast, the case of aj null coefficients (Davis and Rabinowitz, 1984) would correspond to a C-type function (function C1 above, Equation 29).

The Python code used for the computations reported in the present work is available from the following GitHub repository: https://github.com/Confareneoclassico/New_estimator_algorithm. A second Matlab® code was used in a separate set of computations limited to one test function (A2):

https://github.com/Confareneoclassico/Variance_SA_estimators_designs_explorativity_economy/tr ee/master/MatlabCode. The agreement of the independent computations coded and run by separate co-authors (SLP and FF) is offered as internal validation of the results presented in this paper.

Each test comparison is repeated 50 times to ensure reproducibility and reduce the stochastic variation in the results. Some of the experiments were also run 500 times to ensure stability. However, no major difference was observed between 50 and 500 repetitions:

https://github.com/Confareneoclassico/Variance_SA_estimators_designs_explorativity_economy/tree/ master/Extra_material. Each repetition uses an equal number of quasi-random rows from the Sobol' matrix with the input factors # uniformly distributed in (0,1). The total cost of the analysis is kept consistent across methods.

For each of the 50 repetitions, the randomisation procedure is based on the column permutations of the QMC matrix. The first thirty-six columns (that correspond to 6 since = 6 is the largest number of multiple-matrix design tests) of the Sobol' sequence are scrambled in each repetition. The k left-most ones are attributed to matrix A, the following k to matrix B and so forth depending on the number of matrices assessed in the estimator.

As customarily done in QMC computations using Sobol' sequences, the column dimension of each matrix is rounded to the nearest power of two: each power of two corresponds to a 'full scan' of the hyperspace for each block of the sequence. Different blocks of size 2 0 ( = 2, 3, … 14) have been tested for selected functions (classes A, B and C) against ) .

Following , the simulation results are presented in terms of the mean absolute error (MAE) versus the total cost of the analysis, where the MAE is defined as follows:

where # is the analytical value of the total sensitivity index, and a # is its estimated value. In other words, the total error over all factors is considered for the difference between the numerical estimate of the index (averaged over all available elementary effects) and its analytic value. The results are plotted on a decimal logarithmic scale.

The best-performing estimator suggested by Glen and Isaacs (2012) , named by them as estimator D3, is compared against Šaltenis with the Saltelli asymmetric design (from now on the Šaltenis asymmetric estimator for short) (Figure 3 ). These results have also been confirmed by the Matlab® implementation. dashed line). Functions: A1, A2, B1, B2, B3, C1, . Python implementation. Glen-Isaacs cannot achieve a very low cost due to its higher cost evaluating the same number of rows evaluation with respect to Šaltenis. As previously discussed, computing # requires couples of points where all factors but # have differing values. The logic of correcting these sets of points for spurious correlations is that we are considering vectors such as ( % ) and _ % (#)`, where = 1,2, … , when computing the correlation r # for the sensitivity index # . If any of these columns is spuriously correlated in the two matrices because of the finite value of , then this spurious correlation should be removed from r # as described in Glen and Isaacs (2012) .

These authors suggest that a similar correction is useful for computing # . However, the Šaltenis asymmetric estimator for # shows better convergence properties for type A and B functions and marginal better convergence properties for type C functions. This can be explained by considering that There are no chances of strong spurious correlations in this case. Looking back to Table 3 or Figure 2 tells us that one should not expect an improvement when changing from the asymmetric to the symmetric case for = 2 because we obtain the same number of effects at the cost of halving the explorativity of the design.

The Owen estimator (Iooss et al., 2020) has higher explorativity, as seen in section 2.4. This estimator makes different use of the set of base matrices , , and to improve the computation of the elementary effects. However, one can see in Figure 4 , that the Šaltenis asymmetric estimator systematically outperforms this estimator. Moving to the case of multiple-matrix-based designs, one would have hoped that moving to > 2 could improve because the decrease in is offset by an increased number of effects as per Table 3, but this does not appear to be the case ( Figure 5 ). The asymmetric design based on just two matrices is still the best option to compute the total sensitivity indices # , even when compared against the symmetric design for the Šaltenis asymmetric estimator and those based on multiple matrices ( Figure 5 ). It would thus appear that is more important than : the increased number of effects is outperformed by the The results of the comparison with the Lamboni estimator are examined in Figure 6 . Note that the twomatrix symmetric design of the Šaltenis estimator corresponds to Lamboni's for this number of matrices. The lower distance of the Lamboni estimator from the frontrunner (the Šaltenis asymmetric estimator) confirms that this design outperforms the multiple-matrix-based design for the Šaltenis estimator (comparing Figure 5 and Figure 6 ). However, it is still beaten by the Šaltenis asymmetric estimator. The result is also confirmed when the difference with the analytic error is measured as the root-mean squared error (not shown). As shown above, the trade-off between and demonstrates that > 2 is not a convenient design choice. Another way to look at these results is to assess them in terms of stars, which are computationally equivalent to the Saltelli asymmetric design , as seen in section 2.4. The basic design is the one where each star is made of k+1 points using 2k coordinates. To increase , one must increase the number of points in the stars, although this results in decreasing since one uses more coordinates of the core. The cases where the number of matrices is greater than 2 fall into this class. This approach led to worse results. In other words, increasing does not seem to pay off.

We have also tried to compute # using matrix alone, i.e., instead of computing # from ( % ) and

_ % (#)`, we used ( % ) and _ ,(%T,) (#)`. In this approach, when the last N-th row is reached, one uses ( @ ) and _ , (#)` for # , i.e., the system closes on itself. However, this approach did not lead to improvements.

Another sampling procedure we have tested to improve the Šaltenis asymmetric estimator consisted of variably investing the computational budget by improving # 's estimation for the subset of the most important factors while devoting less computational resources to the least important ones for each subsequent model execution. In this adaptive sampling strategy, the choices of the factor to estimate are made using increasing blocks of power of two to fully take advantage of the properties of the lowdiscrepancy Sobol' sequence.

The number of design parameters of the adaptive sampling strategy is reduced to a modicum. Let one assume that one has ) = ( + 1)2 0 model runs available.

• The algorithm is run as per the Šaltenis asymmetric estimator to 'warm up' to a sample size =2 0T.-1 at a cost of ( + 1)2 0T.-1 .

• The factors are then ordered in decreasing order of the standard deviation of the elementary effect bb # .

• At every following block of rows =2 cT0T.-1 (with in the range 1, − 1), it is decided whether the computation of the elementary effect can be stopped at the ( − ) &' factor per order of decreasing importance, thus saving runs.

• The condition upon which this decision is made is the ratio between bb # -s. It is assumed that bb # would be reduced by a factor of √2 by doubling the sample size.

• Hence, the main assumption of the computation is that if bb 1-c-,

• > bb 1-c , this latter factor (and all those having lower bb # s) can be removed from the calculation in the following block.

The computational details are available from the dedicated Jupyter notebook. The results are here presented for test functions A1 and A2, for which # differs across parameters. Another type of function G has also been tested in this experiment, where A3 ( # = {1 2 4 8 16 32}) corresponds to various degrees of importance across parameters.

In Figure 7 , one can appreciate how our method outperforms the Šaltenis asymmetric estimator by up to a factor of two for functions A2 and A3. In this context, the importance of input variables on the output uncertainty can be easily disentangled due to the difference in magnitude across factors' # . This is the setting of a typical real-world model, where the importance of the input factors on the output uncertainty obeys the Pareto principle (Pareto, 1906) with few factors responsible for most of the output variance. However, the case where the sensitivity indices of the input factors are closer in magnitude is more challenging. This adaptive sampling strategy does not outperform the Šaltenis asymmetric estimator in the case of function A1 (Figure 7 ). 

Taking the works of Glen and Isaacs (2012) , Lamboni (2018) , and as our points of departure, we have explored different estimators to improve the computation of the total-effect index # for independent factors using a taxonomy of test functions proposed by . We have seen that the estimator of Glen and Isaacs (2012) is outperformed by Šaltenis and Dzemyda (1982) and the Saltelli asymmetric design . Furthermore, we did not observe improvements in the computational results by extending the symmetric matrix arrangement to values > 2. The larger number of effects obtained with > 2 does not compensate for the loss of explorativity, as is also evidenced by our discrepancy calculation. The increase in economy by using more matrices is offset by the loss of explorativity due to the higher share of repeated coordinates.

To increase the explorativity, one would need to rely on a 'stars' design with centres having less than k rays, which decreases the economy. The latter approach has led to an improvement in the setting of factors receiving a number of estimates proportional to their importance. However, this comes at the cost of introducing an extra design parameter. 

Arnald Puy has worked on this paper on a Marie Sklodowska-Curie Global Fellowship

Handbook of Design and Analysis of Experiments

Sensitivity Analysis: an Introduction for the Management Scientist

Sensitivity analysis: A review of recent advances

Advances in sensitivity analysis

Algorithm 659: Implementing Sobol's Quasirandom Sequence Generator

Sobol' sequences in 16384 dimension

Theory: Theory Vol

From screening to quantitative sensitivity analysis. A unified approach

Methods of Numerical Integration. Courier Corporation

Uncertainty in Industrial Practice: A Guide to Quantitative Uncertainty Management

Design and Modeling for Computer Experiments

Trends in sensitivity analysis practice in the last decade

Estimating Sobol sensitivity indices using correlations

Sensitivity Analysis of Model Output: SAMO

Importance measures in global sensitivity analysis of nonlinear models

Uncertainty Management in Simulation-Optimization of Complex Systems: Algorithms and Applications

2020. sensitivity: Global Sensitivity Analysis of Model Outputs

Monte Carlo Methods in Finance

Ten iterative steps in development and evaluation of environmental models

Asymptotic normality and efficiency of two Sobol index estimators

Analysis of variance designs for model output

Exploring multi-dimensional spaces: a Comparison of Latin Hypercube and Quasi Monte Carlo Sampling Techniques

The identification of model effective dimensions using global sensitivity analysis

Global sensitivity analysis: a generalized, unbiased and optimal estimator of total-effect variance

Sensitivity analysis of spatial models

Nutrition and public health economic evaluations under the lenses of post normal science

Addressing factors fixing setting from given data: A comparison of different methods

Comparison of some efficient methods to evaluate the main effect of computer model factors

Evaluating prediction uncertainty (No. NUREG/CR--6311). Nuclear Regulatory Commission

Comparison of sensitivity analysis methods for pollutant degradation modelling: A case study from drinking water treatment

An introduction to sensitivity assessment of simulation models

Proposed Risk Assessment Bulletin

Better estimation of small sobol' sensitivity indices

Observations on the practice and profession of modeling and simulation: A survey approach

Manual of Political Economy: A Critical and Variorum Edition

Mathematical contributions to the theory of evolution.-On homotyposis in homologous but differentiated organs

Sensitivity analysis of environmental models: A systematic review with practical workflow

Useless Arithmetic: Why Environmental Scientists Can't Predict the Future

How to compute variance-based sensitivity indicators with your spreadsheet software

An effective algorithm for computing global sensitivity indices (EASI)

Global sensitivity measures from given data

Variance-Based Sensitivity Analysis: Theory and Estimation Algorithms

State Dependent Parameter metamodelling and sensitivity analysis

Making best use of model evaluations to compute sensitivity indices

Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices

How to avoid a perfunctory sensitivity analysis

Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index

Five ways to ensure that models serve society: a manifesto

Sensitivity Analysis: Gauging the Worth of Scientific Models

What is science's crisis really about? Futures, Post-Normal science in practice 91

When All Models Are Wrong

What is wrong with evidence based policy, and how can it be improved? Futures, Post-Normal science in practice 91

What do I make of your latinorum? Sensitivity auditing of mathematical modelling

Global Sensitivity Analysis: The Primer

Update 1 of: Sensitivity Analysis for Chemical Models

About the use of rank transformation in sensitivity analysis of model output

Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models

A Quantitative Model-Independent Method for Global Sensitivity Analysis of Model Output

Structure analysis of extremal problems using an approximation of characteristics

Making sense of science for policy under conditions of complexity and uncertainty

Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates

Sensitivity analysis for non-linear mathematical models

Uniformly distributed sequences with an additional uniform property

On the distribution of points in a cube and the approximate evaluation of integrals

Construction and Comparison of High-Dimensional Sobol' Generators. Wilmott

Estimating the approximation error when fixing unessential factors in global sensitivity analysis

Don't try to predict COVID-19. If you must, use Deep Uncertainty methods. Review of Artificial Societies and Social Simulation

Random balance designs for the estimation of first order global sensitivity indices

SAMO 2001: methodological advances and innovative applications of sensitivity analysis

Sensitivity analysis in chemical kinetics

Guidance Document on the Development, Evaluation, and Application of Environmental Models

How to Conduct a Proper Sensitivity Analysis in Life Cycle Assessment: Taking into Account Correlations within LCI Data and Interactions within the LCA Calculation Model

 Table A1 -Terms composing the Glen & Isaacs (2012) 

Explicit formula