# A comprehensive comparison of total-order estimators for global sensitivity
analysis

Arnald Puy Corresponding author Department of Ecology and Evolutionary
Biology, M31 Guyot Hall, Princeton University, New Jersey 08544, USA. E-Mail:
<EMAIL_ADDRESS>Centre for the Study of the Sciences and the Humanities
(SVT), University of Bergen, Parkveien 9, PB 7805, 5020 Bergen, Norway.
William Becker  European Commission, Joint Research Centre, Via Enrico Fermi,
2749, 21027 Ispra VA, Italy Samuele Lo Piano  School of the Built Environment,
JJ Thompson Building, University of Reading, Whiteknights Campus, Reading, RG6
6AF, United Kingdom Andrea Saltelli  Open Evidence Research, Universitat
Oberta de Catalunya (UOC), Barcelona, Spain.

###### Abstract

Sensitivity analysis helps identify which model inputs convey the most
uncertainty to the model output. One of the most authoritative measures in
global sensitivity analysis is the Sobol’ total-order index, which can be
computed with several different estimators. Although previous comparisons
exist, it is hard to know which estimator performs best since the results are
contingent on the benchmark setting defined by the analyst (the sampling
method, the distribution of the model inputs, the number of model runs, the
test function or model and its dimensionality, the weight of higher order
effects or the performance measure selected). Here we compare several total-
order estimators in an eight-dimension hypercube where these benchmark
parameters are treated as random parameters. This arrangement significantly
relaxes the dependency of the results on the benchmark design. We observe that
the most accurate estimators are Razavi and Gupta’s, Jansen’s or Janon/Monod’s
for factor prioritization, and Jansen’s, Janon/Monod’s or Azzini and Rosati’s
for approaching the “true” total-order indices. The rest lag considerably
behind. Our work helps analysts navigate the myriad of total-order formulae by
reducing the uncertainty in the selection of the most appropriate estimator.

Keywords: Uncertainty analysis; sensitivity analysis; modeling; Sobol’
indices; variance decomposition, benchmarking analysis

##  1 Introduction

Sensitivity analysis, i.e. the assessment of how much uncertainty in a given
model output is conveyed by each model input, is a fundamental step to judge
the quality of model-based inferences [1, 2, 3]. Among the many sensitivity
indices available, variance-based indices are widely regarded as the gold
standard because they are model-free (no assumptions are made about the
model), global (they account for interactions between the model inputs) and
easy to interpret [4, 5, 6]. Given a model of the form $y=f(\bm{x})$,
$\bm{x}=(x_{1},x_{2},...,x_{i},...,x_{k})\in\mathbb{R}^{k}$, where $y$ is a
scalar output and $x_{1},...,x_{k}$ are the $k$ independent model inputs, the
variance of $y$ is decomposed into conditional terms as

$V(y)=\sum_{i=1}^{k}V_{i}+\sum_{i}\sum_{i<j}V_{ij}+...+V_{1,2,...,k}\\\ ,$ (1)

where

$\displaystyle V_{i}=V_{x_{i}}\big{[}E_{\bm{x}_{\sim i}}(y|x_{i})\big{]}\hskip
11.38109ptV_{ij}$ $\displaystyle=V_{x_{i},x_{j}}\big{[}E_{\bm{x}_{\sim
i,j}}(y|x_{i},x_{j})\big{]}\hskip 11.38109pt$ (2) $\displaystyle-
V_{x_{i}}\big{[}E_{\bm{x}_{\sim i}}(y|x_{i})\big{]}$ $\displaystyle-
V_{x_{j}}\big{[}E_{\bm{x}_{\sim j}}(y|x_{j})\big{]}$

and so on up to the $k$-th order. The notation $\bm{x}_{\sim i}$ means _all-
but- $x_{i}$_. By dividing each term in Equation 1 by the unconditional model
output variance $V(y)$, we obtain the first-order indices for single inputs
($S_{i}$), pairs of inputs ($S_{ij}$), and for all higher-order terms. First-
order indices thus provide the proportion of $V(y)$ caused by each term and
are widely used to rank model inputs according to their contribution to the
model output uncertainty, a setting known as _factor prioritization_ [1].

[7] also proposed the calculation of the total-order index $T_{i}$, which
measures the first-order effect of a model input jointly with its interactions
up to the $k$-th order:

$T_{i}=1-\frac{V_{\bm{x}_{\sim i}}\big{[}E_{x_{i}}(y|\bm{x}_{\sim
i})\big{]}}{V(y)}=\frac{E_{\bm{x}_{\sim i}}\big{[}V_{x_{i}}(y|\bm{x}_{\sim
i})\big{]}}{V(y)}\,.$ (3)

When $T_{i}\approx 0$, it can be concluded that $x_{i}$ has a negligible
contribution to $V(y)$. For this reason, total-order indices have been applied
to distinguish influential from non-influential model inputs and reduce the
dimensionality of the uncertain space, a setting known as _factor-fixing_ [1].

The most direct computation of $T_{i}$ is via Monte Carlo (MC) estimation
because it does not impose any assumption on the functional form of the
response function, unlike metamodeling approaches [8, 9]. The Fourier
Amplitude Sensitivity Test (FAST) may also be used to calculate $T_{i}$, which
involves transforming input variables into periodic functions of a single
frequency variable, sampling the model and analysing the sensitivity of input
variables using Fourier analysis in the frequency domain [10, 11]. While an
innovative approach, FAST is sensitive to the characteristic frequencies
assigned to input variables, and is not a very intuitive method - for these
reasons it has mostly been superseded by Monte Carlo approaches, or by
metamodels when computational expense is a serious issue. In this work we
focus on the former.

MC methods require generating a $(N,2k)$ base sample matrix with either random
or quasi-random numbers (e.g. Latin Hypercube Sampling, Sobol’ quasi-random
numbers [12, 13]), where each row is a sampling point and each column a model
input. The first $k$ columns are allocated to an $\bm{A}$ matrix and the
remaining $k$ columns to a $\bm{B}$ matrix, which are known as the “base
sample matrices”. Any point in either $\bm{A}$ or $\bm{B}$ can be indicated as
$x_{vi}$, where $v$ and $i$ respectively index the row (from 1 to $N$) and the
column (from 1 to $k$). Then, $k$ additional $\bm{A}_{B}^{(i)}$
($\bm{B}_{A}^{(i)}$) matrices are created, where all columns come from
$\bm{A}$ ($\bm{B}$) except the $i$-th column, which comes from $\bm{B}$
($\bm{A}$). The numerator in Equation 3 is finally estimated using the model
evaluations obtained from the $\bm{A}$ ($\bm{B}$) and $\bm{A}_{B}^{(i)}$
($\bm{B}_{A}^{(i)}$) matrices. Some estimators may also use a third or
$\bm{X}$ base sample matrices (i.e. $\bm{A},\bm{B},\bm{C},\ldots,\bm{X}$),
although the use of more than three matrices has been recently proven
inefficient by [14].

###  1.1 Total-order estimators and uncertainties in the benchmark settings

The search for efficient and robust total-order estimators is an active field
of research [7, 15, 1, 16, 17, 18, 19, 20]. Although some works have compared
their asymptotic properties (i.e. [16]), most studies have promoted empirical
comparisons where different estimators are benchmarked against known test
functions and specific sample sizes. However valuable these empirical studies
may be, [21] observed that their results are very much conditional on the
choice of model, its dimensionality and the selected number of model runs. It
is hard to say from previous studies whether an estimator outperforming
another truly reflects its higher accuracy or simply its better performance
under the narrow statistical design of the study. Below we extend the list of
factors which [21] regards as influential in a given benchmarking exercise and
discuss how they affect the relative performance of sensitive estimators.

  * •

The sampling method: The creation of the base sample matrices can be done
using Monte-Carlo (MC) or quasi Monte-Carlo (QMC) methods [12, 13]. Compared
to MC, QMC allows to more effectively map the input space as it leaves smaller
unexplored volumes (Fig. S1). However, [22] observed that MC methods might
help obtain more accurate sensitivity indices when the model under examination
has important high-order terms. Both MC and QMC have been used when
benchmarking sensitivity indices [15, 23].

  * •

The form of the test function: some of the most commonly used functions in SA
are the [24]’s, the Sobol’ G and its variants [25, 23], the [26]’s or the set
of functions presented in [22] [16, 18, 23, 14]. Despite being analytically
tractable, these functions capture only one possible interval of model
behaviour, and the effects of nonlinearities and nonadditivities is typically
unknown in real models. This _black-box_ nature of models has become more of a
concern in the last decades due to the increase in computational power and
code complexity (which prevents the analyst from intuitively grasping the
model’s behaviour [27]), and to the higher demand for model transparency [3,
28, 29]. This renders the functional form of the model similar to a random
variable [21], something not accounted for by previous works [16, 18, 23, 14].

  * •

The function dimensionality: many studies focus on low-dimensional problems,
either by using test functions that only require a few model inputs (e.g. the
Ishigami function, where $k=3$), or by using test functions with a flexible
dimensionality, but setting $k$ at a small value of e.g. $k\leq 8$ ([25]’s G
or [26] functions). This approach trades computational manageability for
comprehensiveness: by neglecting higher dimensions, it is difficult to tell
which estimator might work best in models with tens or hundreds of parameters.
Examples of such models can be readily found in the Earth and Environmental
Sciences domain [30], including the Soil and Water Assessment Tool (SWAT)
model, where $k=50$ [31], or the Modélisation Environmentale-Surface et
Hydrologie (MESH) model, where $k=111$ [32].

  * •

The distribution of the model inputs: the large majority of benchmarking
exercises assume uniformly-distributed inputs $p(\bm{x})\in U(0,1)^{k}$ [16,
33, 23, 14]. However, there is evidence that the accuracy of $T_{i}$
estimators might be sensitive to the underlying model input distributions, to
the point of overturning the model input ranks [34, 35]. Furthermore, in
uncertainty analysis – e.g. in decision theory, the analysts may use
distributions with peaks for the most likely values derived, for instance,
from an experts elicitation stage.

  * •

The number of model runs: sensitivity test functions are generally not
computationally expensive and can be run without much concern for
computational time. This is frequently not the case for real models, whose
high dimensionality and complexity might set a constraint on the total number
of model runs available. Under such restrictions, the performance of the
estimators of the total-order index depends on their efficiency (how accurate
they are given the budget of runs that can be allocated to each model input).
There are no specific guidelines as to which total-order estimator might work
best under these circumstances [21].

  * •

The performance measure selected: typically, a sensitivity estimator has been
considered to outperform the rest if, on average, it displays a smaller mean
absolute error (MAE), computed as

$\text{MAE}=\frac{1}{p}\sum_{v=1}^{p}\left(\frac{\sum_{i=1}^{k}|T_{i}-\hat{T}_{i}|}{k}\right)\,,$
(4)

where $p$ is the number of replicas of the sample matrix, and $T_{i}$ and
$\hat{T}_{i}$ the analytical and the estimated total-order index of the $i$-th
input. The MAE is appropriate when the aim is to assess which estimator better
approaches the true total-order indices, because it averages the error for
both influential and non-influential indices. However, the analyst might be
more interested in using the estimated indices
$\bm{\hat{T}}=\\{\hat{T}_{1},\hat{T}_{2},...,\hat{T}_{i},...,\hat{T}_{k}\\}$
to accurately rank parameters or screen influential from non-influential model
inputs [1]. In such context, the MAE may be best substituted or complemented
with a measure of rank concordance between the vectors $\bm{r}$ and
$\bm{\hat{r}}$, which reflect the ranks in $\bm{T}$ and $\bm{\hat{T}}$
respectively, such as the Spearman’s $\rho$ or the Kendall’s $W$ coefficient
[36, 37, 21]. It can also be the case that disagreements on the exact ranking
of low-ranked parameters may have no practical importance because the interest
lies in the correct identification of top ranks only [30]. [38] scores or
other measures that emphasize this top-down correlation are then a more
suitable choice.

Here we benchmark the performance of eight different MC-based formulae
available to estimate $T_{i}$ (Table 1). While the list is not exhaustive,
they reflect the research conducted on $T_{i}$ over the last 20 years: from
the classic estimators of [7, 15, 1] up to the new contributions by [16],
[17], [33] and [20, 39]. In order to reduce the influence of the benchmarking
design in the assessment of the estimators’ accuracy, we treat the sampling
method $\tau$, the underlying model input distribution $\phi$, the number of
model runs $N_{t}$, the test function $\varepsilon$, its dimensionality and
degree of non-additivity ($k,k_{2},k_{3}$) and the performance measure
$\delta$ as random parameters. This better reflects the diversity of models
and sensitivity settings available to the analyst. By relaxing the dependency
of the results on these benchmark parameters111We refer to the set of
benchmarking assumptions as _benchmarking parameters_ or _parameters_. This is
intended to distinguish them from the inputs of each test function generated
by the metafunction, which we refer to as inputs., we define an
unprecedentedly large setting where all formulae can prove their accuracy. We
therefore extend [21]’s approach by testing a wider set of Monte Carlo
estimators, by exploring a wider range of benchmarking assumptions and by
performing a formal SA on these assumptions. The aim is therefore to provide a
much more global comparison of available MC estimators than is available in
the existing literature, and investigate how the benchmarking parameters may
affect the relative performance of estimators. Such information can help point
to estimators that are not only efficient on a particular case study, but
efficient and robust to a wide range of practical situations.

Table 1: Formulae to compute $T_{i}$. $f_{0}$ and $V(y)$ are estimated according to the original papers. For estimators 2 and 5, $f_{0}=\frac{1}{N}\sum_{v=1}^{N}f(\bm{A})_{v}$. For estimators 1, 2 and 5, $V(y)=\frac{1}{N}\sum_{v=1}^{N}\left[f(\bm{A})_{v}-f_{0}\right]^{2}$ [Eq. 4.16]Saltelli2008[Eqs. 15, 20]Homma1996a. For estimator 3, $f_{0}=\frac{1}{N}\sum_{v=1}^{N}\frac{f(\bm{A})_{v}+f(\bm{A}_{B}^{(i)})_{v}}{2}$ and $V(y)=\frac{1}{N}\sum_{v=1}^{N}\frac{f(\bm{A})_{v}^{2}+f(\bm{A}_{B}^{(i)})_{v}^{2}}{2}-f_{0}^{2}$ [16, Eq. 15]. In estimator 4, $\langle f(\bm{A})_{v}\rangle$ is the mean of $f(\bm{A})_{v}$. We use a simplified version of the Glen and Isaacs estimator because spurious correlations are zero by design. As for estimator 7, we refer to it as pseudo-Owen given its use of a $\bm{C}$ matrix and its identification with [40] in [41], where we retrieve the formula from. $V(y)$ in Estimator 7 is computed as in Estimator 3 following [41], whereas $V(y)$ in Estimator 8 is computed as in Estimator 1. Nº | Estimator |  Author
---|---|---
1 | $\frac{\frac{1}{2N}\sum_{v=1}^{N}\left[f(\bm{A})_{v}-f(\bm{A}_{B}^{(i)})_{v}\right]^{2}}{V(y)}$ |  [15]
2 | $\frac{V(y)-\frac{1}{N}\sum_{v=1}^{N}f(\bm{A})_{v}f(\bm{A}_{B}^{(i)})_{v}+f_{0}^{2}}{V(y)}$ |  [7]
3 | $1-\frac{\frac{1}{N}\sum_{v=1}^{N}f(\bm{A})_{v}f(\bm{A}_{B}^{(i)})_{v}-f_{0}^{2}}{V(y)}$ |  [16]
[19]
4 | $1-\left[\frac{1}{N-1}\sum_{v=1}^{N}\frac{\left[f(\bm{A})_{v}-\left\langle f(\bm{A})_{v}\right\rangle\right]\left[f(\bm{A}_{B}^{(i)})_{v}-\left\langle f(\bm{A}_{B}^{(i)})_{v}\right\rangle\right]}{\sqrt{V\left[f(\bm{A})_{v}\right]V\left[f(\bm{A}_{B}^{(i)})_{v}\right]}}\right]$ |  [17]
5 | $1-\frac{\frac{1}{N}\sum_{v=1}^{N}f(\bm{B})_{v}f(\bm{B}^{(i)}_{A})_{v}-f_{0}^{2}}{V(y)}$ |  [1]
6 | $\frac{\sum_{v=1}^{N}[f(\bm{B})_{v}-f(\bm{B}^{(i)}_{A})_{v}]^{2}+[f(\bm{A})_{v}-f(\bm{A}^{(i)}_{B})_{v}]^{2}}{\sum_{v=1}^{N}[f(\bm{A})_{v}-f(\bm{B})_{v}]^{2}+[f(\bm{B}^{(i)}_{A})_{v}-f(\bm{A}^{(i)}_{B})_{v}]^{2}}$ |  [33, 18]
7 | $\frac{V(y)-\left[\frac{1}{N}\sum_{v=1}^{N}\left\\{\left[f(\bm{B})_{v}-f(\bm{C}_{B}^{(i)})_{v}\right]\left[f(\bm{B}_{A}^{(i)})_{v}-f(\bm{A})_{v}\right]\right\\}\right]}{V(y)}$ |  pseudo-Owen
8 | $\frac{E_{x^{*}_{\sim_{i}}}\left[\gamma_{x^{*}{{}_{\sim}i}}(h_{i})\right]+E_{x^{*}{{}_{\sim}i}}\left[C_{x^{*}{{}_{\sim}i}}(h_{i})\right]}{V(y)}$ |  [20, 39] (see SM).

##  2 Assessment of the uncertainties in the benchmarking parameters

In this section we formulate the benchmarking parameters as random variables
and assess how the performance of estimators is dependent on them by
performing a sensitivity analysis. In essence this is a _sensitivity analysis
of sensitivity analyses_ [42], and a natural extension of a similar
uncertainty analysis in a recent work by [21]. The use of global sensitivity
analysis tools to better understand the properties of estimators can give
insights into how estimators behave in different scenarios that are not
available through analytical approaches.

###  2.1 The setting

The variability in the benchmark settings
($\tau,N_{t},k,k_{2},k_{3},\phi,\epsilon,\delta$) is described by probability
distributions (Table 2). We assign uniform distributions (discrete or
continuous) to each parameter. In particular, we choose
$\tau\sim\mathcal{DU}(1,2)$ to check how the performance of $T_{i}$ estimators
is conditioned by the use of Monte-Carlo ($\tau=1$) or Quasi Monte-Carlo
($\tau=2$) methods in the creation of the base sample matrices. For $\tau=2$
we use the Sobol’ sequence scrambled according to [43] to avoid repeated
coordinates at the beginning of the sequence. The total number of model runs
and inputs is respectively described as $N_{t}\sim\mathcal{DU}(10,1000)$ and
$k\sim\mathcal{DU}(3,100)$ to explore the performance of the estimators in a
wide range of $N_{t},k$ combinations. Given the sampling constraints set by
the estimators’ reliance on either a $\bm{B}$, $\bm{B}_{A}^{(i)}$,
$\bm{A}_{B}^{(i)}$ or $\bm{C}_{B}^{(i)}$ matrices (Table 1), we modify the
space defined by ($N_{t},k$) to a non-rectangular domain (we provide more
information on this adjustment in Section 2.2).

Table 2: Summary of the parameters and their distributions. $\mathcal{DU}$ stands for discrete uniform. Parameter | Description | Distribution
---|---|---
$\tau$ | Sampling method | $\mathcal{DU}(1,2)$
$N_{t}$ | Total number of model runs | $\mathcal{DU}(10,1000)$
$k$ | Number of model inputs | $\mathcal{DU}(3,100)$
$\phi$ | Probability distribution of the model inputs | $\mathcal{DU}(1,8)$
$\varepsilon$ | Randomness in the test function | $\mathcal{DU}(1,200)$
$k_{2}$ | Fraction of pairwise interactions | $\mathcal{U}(0.3,0.5)$
$k_{3}$ | Fraction of three-wise interactions | $\mathcal{U}(0.1,0.3)$
$\delta$ | Selection of the performance measure | $\mathcal{DU}(1,2)$

For $\phi$ we set $\phi\sim\mathcal{DU}(1,8)$ to ensure an adequate
representation of the most common shapes in the $(0,1)^{k}$ domain. Besides
the normal distribution truncated at $(0,1)$ and the uniform distribution, we
also take into account four beta distributions parametrized with distinct
$\alpha$ and $\beta$ values and a logitnormal distribution (Fig. 1a). The aim
is to check the response of the estimators under a wide range of probability
distributions, including U-shaped distributions and distributions with
different degrees of skewness.

Figure 1: The metafunction approach. a) Probability distributions incuded in
$\phi$. $N_{T}$ stands for truncated normal distribution. b) Univariate
functions included in the metafunction ($f_{1}(x)=$ cubic, $f_{2}(x)=$
discontinuous, $f_{3}(x)=$ exponential, $f_{4}(x)=$ inverse, $f_{5}(x)=$
linear, $f_{6}(x)=$ no effect, $f_{7}(x)=$ non-monotonic, $f_{8}(x)=$
periodic, $f_{9}(x)=$ quadratic, $f_{10}(x)=$ trigonometric).

We link each distribution in Fig. 1a to an integer value from 1 to 7. For
instance, if $\phi=1$, the joint probability distribution of the model inputs
is described as $p(x_{1},\ldots,x_{k})=\mathcal{U}(0,1)^{k}$. If $\phi=8$, we
create a vector $\bm{\phi}=\\{\phi_{1},\phi_{2},...,\phi_{i},...,\phi_{k}\\}$
by randomly sampling the seven distributions in Fig. 1a, and use the $i$-th
distribution in the vector to describe the uncertainty of the $i$-th input.
This last case examines the behavior of the estimators when several
distributions are used to characterize the uncertainty in the model input
space.

####  2.1.1 The test function

The parameter $\varepsilon$ operationalizes the randomness in the form and
execution of the test function. Our test function is an extended version of
[21]’s metafunction, which randomly combines $p$ univariate functions in a
multivariate function of dimension $k$. Here we consider the 10 univariate
functions listed in Fig. 1b, which represent common responses observed in
physical systems and in classic SA test functions (see [21] for a discussion
on this point). We note that an alternative approach would be to construct
orthogonal basis functions which could allow analytical evaluation of true
sensitivity indices for each generated function; however, this extension is
left for future work.

We construct the test function as follows:

  1. 1.

Let us consider a sample matrix such as

$\bm{M}=\begin{bmatrix}x_{11}&x_{12}&\cdots&x_{1i}&\cdots&x_{1k}\\\
x_{21}&x_{22}&\cdots&x_{2i}&\cdots&x_{2k}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
x_{v1}&x_{v2}&\cdots&x_{vi}&\cdots&x_{vk}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
x_{N1}&x_{N2}&\cdots&x_{Ni}&\cdots&x_{Nk}\\\ \end{bmatrix}$ (5)

where every point $\bm{x}_{v}=x_{v1},x_{v2},\ldots,x_{vk}$ represents a given
combination of values for the $k$ inputs and $x_{i}$ is a model input whose
distribution is defined by $\phi$.

  2. 2.

Let $\bm{u}=\\{u_{1},u_{2},...,u_{k}\\}$ be a $k$-length vector formed by
randomly sampling with replacement the ten functions in Fig. 1b. The $i$-th
function in $\bm{u}$ is then applied to the $i$-th model input: for instance,
if $k=4$ and $\bm{u}=\\{u_{3},u_{4},u_{8},u_{1}\\}$, then
$f_{3}(x_{1})=\frac{e^{x_{1}}-1}{e-1}$,
$f_{4}(x_{2})=(10-\frac{1}{1.1})^{-1}(x_{2}+0.1)^{-1}$,
$f_{8}(x_{3})=\frac{\sin(2\pi x_{3})}{2}$, and $f_{1}(x_{4})=x_{4}^{3}$. The
elements in $\bm{u}$ thus represent the first-order effects of each model
input.

  3. 3.

Let $\bm{V}$ be a $(n,2)$ matrix, for $n=\frac{k!}{2!(k-2)!}$, the number of
pairwise combinations between the $k$ inputs of the model. Each row in
$\bm{V}$ thus specifies an interaction between two columns in $\bm{M}$. In the
case of $k=4$ and the same elements in $\bm{u}$ as defined in the previous
example,

$\bm{V}=\begin{bmatrix}1&2\\\ 1&3\\\ 1&4\\\ 2&3\\\ 2&4\\\ 3&4\\\
\end{bmatrix}$ (6)

e.g., the first row promotes $f_{3}(x_{1})\cdot f_{4}(x_{2})$, the second row
$f_{3}(x_{1})\cdot f_{8}(x_{3})$, and so on until the $n$-th row. In order to
follow the _sparsity of effects principle_ (most variations in a given model
output should be explained by low-order interactions [44]), the metafunction
activates only a fraction of these effects: it randomly samples $\llceil
k_{2}n\rrceil$ rows from $\bm{V}$, and computes the corresponding interactions
in $\bm{M}$. $\llceil k_{2}n\rrceil$ is thus the number of pairwise
interactions present in the function. We make $k_{2}$ an uncertain parameter
described as $k_{2}\sim\mathcal{U}(0.3,0.5)$ in order to randomly activate
only between 30% and 50% of the available second-order effects in $\bm{M}$.

  4. 4.

Same as before, but for third-order effects: let $\bm{W}$ be a ($m,3$) matrix,
for $m=\frac{k!}{3!(k-3)!}$, the number of three-wise combinations between the
$k$ inputs in $\bm{M}$. For $k=4$ and $\bm{u}$ as before,

$\bm{W}=\begin{bmatrix}1&2&3\\\ 1&2&4\\\ 1&3&4\\\ 2&3&4\\\ \end{bmatrix}$ (7)

e.g. the first row leads to $f_{3}(x_{1})\cdot f_{4}(x_{2})\cdot
f_{8}(x_{3})$, and so on until the $m$-th row. The metafunction then randomly
samples $\llceil k_{3}m\rrceil$ rows from $\bm{W}$ and computes the
corresponding interactions in $\bm{M}$. $\llceil k_{3}m\rrceil$ is therefore
the number of three-wise interaction terms in the function. We also make
$k_{3}$ an uncertain parameter described as $k_{3}\sim\mathcal{U}(0.1,0.3)$ to
activate only between 10% and 30% of all third-order effects in $\bm{M}$. Note
that $k_{2}>k_{3}$ because third-order effects tend to be less dominant than
two-order effects (Table 2).

  5. 5.

Three vectors of coefficients ($\bm{\alpha},\bm{\beta},\bm{\gamma}$) of length
$k$, $n$ and $m$ are defined to represent the weights of the first, second and
third-order effects respectively. These coefficients are generated by sampling
from a mixture of two normal distributions
$\Psi=0.3\mathcal{N}(0,5)+0.7\mathcal{N}(0,0.5)$. This coerces the
metafunction into replicating the [45] principle (around 80% of the effects
are due to 20% of the parameters), found to widely apply in SA [46, 1].

  6. 6.

The metafunction can thus be formalized as

$\displaystyle y=$
$\displaystyle\sum_{i=1}^{k}\alpha_{i}f^{u_{i}}\phi_{i}(x_{i})$ (8)
$\displaystyle+\sum_{i=1}^{\llceil
k_{2}n\rrceil}\beta_{i}f^{u_{V_{i,1}}}\phi_{i}(x_{V_{i,1}})f^{u_{V_{i,2}}}\phi_{i}(x_{V_{i,2}})$
$\displaystyle+\sum_{i=1}^{\llceil
k_{3}m\rrceil}\gamma_{i}f^{u_{W_{i,1}}}\phi_{i}(x_{W_{i,1}})f^{u_{W_{i,2}}}\phi_{i}(x_{W_{i,2}})f^{u_{W_{i,3}}}\phi_{i}(x_{W_{i,3}})\,.$

Note that there is randomness in the sampling of $\bm{\phi}$, the univariate
functions in $\bm{u}$ and the coefficients in
$(\bm{\alpha},\bm{\beta},\bm{\gamma})$. The parameter $\varepsilon$ assesses
the influence of this randomness by fixing the starting point of the pseudo-
random number sequence used for sampling the parameters just mentioned. We use
$\varepsilon\sim\mathcal{U}(1,200)$ to ensure that the same seed does not
overlap with the same value of $N_{t}$, $k$ or any other parameter, an issue
that might introduce determinism in a process that should be stochastic. In
Figs. S2–S3we show the type of $T_{i}$ indices generated by this metafunction.

Finally, we describe the parameter $\delta$ as $\delta\sim\mathcal{DU}(1,2)$.
If $\delta=1$, we compute the Kendall $\tau$-b correlation coefficient between
$\bm{\hat{r}}$ and $\bm{r}$, the estimated and the “true” ranks calculated
from $\bm{\hat{T}}$ and $\bm{T}$ respectively. This aims at evaluating how
well the estimators in Table 1 rank all model inputs. If $\delta=2$, we
compute the Pearson correlation between $\bm{r}$ and $\bm{\hat{r}}$ after
transforming the ranks to Savage scores [38]. This setting examines the
performance of the estimators when the analyst is interested in ranking only
the most important model inputs. Savage scores are given as

$Sa_{i}=\sum_{j=i}^{k}\frac{1}{j}\,,$ (9)

where $j$ is the rank assigned to the $j$th element of a vector of length $k$.
If $x_{1}>x_{2}>x_{3}$, the Savage scores would then be
$Sa_{1}=1+\frac{1}{2}+\frac{1}{3}$, $Sa_{2}=\frac{1}{2}+\frac{1}{3}$, and
$Sa_{3}=\frac{1}{3}$. The parameter $\delta$ thus assesses the accuracy of the
estimators in properly ranking the model inputs; in other words, when they are
used in a factor prioritization setting [1].

In order to examine also how accurate the estimators are in approaching the
“true” indices, we run an extra round of simulations with the MAE as the only
performance measure, which we compute as

$\text{MAE}=\frac{\sum_{i=1}^{k}|T_{i}-\hat{T}_{i}|}{k}\,.$ (10)

Note that, unlike Equation 4, Equation 10 does not make use of replicas. This
is because the effect of the sampling is averaged out in our design by
simultaneously varying all parameters in many different simulations.

###  2.2 The execution of the algorithm

We examine how sensitive the performance of total-order estimators is to the
uncertainty in the benchmark parameters
$\tau,N_{t},k,k_{2},k_{3},\phi,\epsilon,\delta$ by means of a global SA. We
create an $\bm{A}$, $\bm{B}$ and $k-1$ $\bm{A}^{(i)}_{B}$ matrices, each of
dimension $(2^{11},k)$, using Sobol’ quasi-random numbers. In these matrices
each column is a benchmark parameter described with the probability
distributions of Table 2 and each row is a simulation with a specific
combination of $\tau,N_{t},k,\ldots$ values. Note that we use $k-1$
$\bm{A}^{(i)}_{B}$ matrices because we group $N_{t}$ and $k$ and treat them
like a single benchmark parameter given their correlation (see below).

Our algorithm runs rowwise over the $\bm{A}$, $\bm{B}$ and $k-1$
$\bm{A}^{(i)}_{B}$ matrices, for $v=1,2,\ldots,18,432$ rows. In the $v$-th row
it does the following:

  1. 1.

It creates five $(N_{t_{v}},k_{v})$ matrices using the sampling method defined
by $\tau_{v}$. The need for these five sub-matrices responds to the five
specific sampling designs requested by the estimators of our study (Table 1).
We use these matrices to compute the vector of estimated indices
$\bm{\hat{T}}_{i}$ for each estimator:

    1. (a)

An $\bm{A}$ matrix and $k_{v}$ $\bm{A}_{B}^{(i)}$ matrices, each of size
$(N_{v},k_{v})$, $N_{v}=\llceil\frac{N_{t{{}_{v}}}}{k_{v}+1}\rrceil$
(Estimators 1–4 in Table 1).

    2. (b)

An $\bm{A}$, $\bm{B}$ and $k_{v}$ $\bm{A}_{B}^{(i)}$ matrices, each of size
$(N_{v},k_{v})$, $N_{v}=\llceil\frac{N_{t{{}_{v}}}}{k_{v}+2}\rrceil$
(Estimator 5 in Table 1).

    3. (c)

An $\bm{A}$, $\bm{B}$ and $k_{v}$ $\bm{A}_{B}^{(i)}$ and $\bm{B}_{A}^{(i)}$
matrices, each of size $(N_{v},k_{v})$,
$N_{v}=\llceil\frac{N_{t_{v}}}{2k_{v}+2}\rrceil$ (Estimator 6 in Table 1).

    4. (d)

An $\bm{A}$, $\bm{B}$ and $k_{v}$ $\bm{B}_{A}^{(i)}$ and $\bm{C}_{B}^{(i)}$
matrices, each of size $(N_{v},k_{v})$,
$N_{v}=\llceil\frac{N_{t{{}_{v}}}}{2k_{v}+2}\rrceil$ (Estimator 7 in Table 1).

    5. (e)

A matrix formed by $N_{v}$ stars, each of size $k_{v}(\frac{1}{\Delta
h}-1)+1$. Given that we set $\Delta h$ at 0.2 (see Supplementary Materials),
$N_{v}=\llfloor\frac{N_{t_{v}}}{4k+1}\rrfloor$ (Estimator 8 in Table 1).

The different sampling designs and the value for $k_{v}$ constrains the total
number of runs $N_{t_{v}}$ that can be allocated to each estimator.
Furthermore, given the probability distributions selected for $N_{t}$ and $k$
(Table 2), specific combinations of ($N_{t_{v}},k_{v}$) lead to $N_{v}\leq 1$,
which is computationally unfeasible. To minimize these issues we force the
comparison between estimators to approximate the same $N_{t_{v}}$ value. Since
the sampling design structure of Razavi and Gupta is the most constraining, we
use $N_{v}=\frac{2(4k+1)}{k+1}$ (for estimators 1–4),
$N_{v}=\frac{2(4k+1)}{k+2}$ (for estimator 5) and $N_{v}=\frac{2(4k+1)}{2k+2}$
(for estimators 6–7) when $N_{v}\leq 1$ in the case of Razavi and Gupta. This
compels all estimators to explore a very similar portion of the ($N_{t},k$)
space, but $N_{t}$ and $k$ become correlated, which contradicts the
requirement of independent inputs characterizing variance-based sensitivity
indices [1]. This is why we treat ($N_{t},k$) as a single benchmark parameter
in the SA.

  2. 2.

It creates a sixth matrix, formed by an $\bm{A}$ and $k_{v}$
$\bm{A}_{B}^{(i)}$ matrices, each of size $(2^{11},k_{v})$. We use this sub-
matrix to compute the vector of “true” indices $\bm{T}$, which could not be
calculated analytically due to the wide range of possible functional forms
created by the metafunction. Following [21], we assume that a fairly accurate
approximation to $\bm{T}$ could be achieved with a large Monte Carlo
estimation.

  3. 3.

The distribution of the model inputs in these six sample matrices is defined
by $\phi_{v}$.

  4. 4.

The metafunction runs over these six matrices simultaneously, with its
functional form, degree of active second and third-order effects as set by
$\varepsilon_{v}$, $k_{2_{v}}$ and $k_{3_{v}}$ respectively.

  5. 5.

It computes the estimated sensitivity indices $\bm{\hat{T}}_{v}$ for each
estimator and the “true” sensitivity indices $\bm{T}_{v}$ using the [15]
estimator, which is currently best practice in SA.

  6. 6.

It checks the performance of the estimators. This is done in two ways:

    1. (a)

If $\delta=1$, we compute the correlation between $\bm{\hat{r}}_{v}$ and
$\bm{r}_{v}$ (obtained respectively from $\bm{\hat{T}}_{v}$ and $\bm{T}_{v}$)
with Kendall tau, and if $\delta=2$ we compute the correlation between
$\bm{\hat{r}}_{v}$ and $\bm{r}_{v}$ on Savage scores. The model output in both
cases is the correlation coefficient $r$, with higher $r$ values indicating a
better performance in properly ranking the model inputs.

    2. (b)

We compute the MAE between $\bm{\hat{T}}_{v}$ and $\bm{T}_{v}$. In this case
the model output is the MAE, with lower values indicating a better performance
in approaching the “true” total-order indices.

##  3 Results

###  3.1 Uncertainty analysis

Under a factor prioritization setting (e.g. when the aim is to rank the model
inputs in terms of their contribution to the model output variance), the most
accurate estimators are Jansen, Razavi and Gupta, Janon/Monod and Azzini and
Rosati. The distribution of $r$ values (the correlation between estimated and
”true” ranks) when these estimators are used is highly negatively skewed, with
median values of $\approx 0.9$. Glen and Isaacs, Homma and Saltelli, Saltelli
and pseudo-Owen lag behind and display median $r$ values of $\approx 0.35$,
with pseudo-Owen ranking last ($r\approx 0.2$). The range of values obtained
with these formulae is much more spread out and include a significant number
of negative $r$ values, suggesting that they overturned the true ranks in
several simulations (Figs. 2a, S4).

Figure 2: Boxplots summarizing the results of the simulations. a) Correlation
coefficient between $\bm{\hat{r}}$ and $\bm{r}$, the vector of estimated and
“true” ranks. b) Mean Absolute Error (MAE).

When the goal is to approximate the “true” indices, Janon/Monod, Jansen and
Azzini and Rosati also offer the best performance. The median MAE obtained
with these estimators is generally smaller than Glen and Isaacs’ and pseudo-
Owen’s, and the distribution of MAE values is much more narrower than that
obtained with Homma and Saltelli, Saltelli or Razavi and Gupta. These three
estimators are the least accurate and produce several MAE values larger than
$10^{2}$ in several simulations (Fig. 2b). The volatility of Razavi and Gupta
under the MAE is reflected in the numerous outliers produced and sharply
contrasts with its very good performance in a factor prioritization setting
(Fig. 2a).

To obtain a finer insight into the structure of these results, we plot the
total number of model runs $N_{t}$ against the function dimensionality $k$
(Fig. 3). This maps the performance of the estimators in the input space
formed by all possible combinations of $N_{t}$ and $k$ given the specific
design constraints of each formulae. Under a factor prioritization setting,
almost all estimators perform reasonably well at a very small dimensionality
($k\leq 10,r>0.7$), regardless of the total number of model runs available.
However, some differences unfold at higher dimensions: Saltelli, Homma and
Saltelli, Glen and Isaacs and especially pseudo-Owen swiftly become inaccurate
for $k>10$, even with large values for $N_{t}$. Azzini and Rosati display a
very good performance overall except in the upper $N_{t},k$ boundary, where
most of the orange dots concentrate. The estimators of Jansen, Janon/Monod and
Razavi and Gupta rank the model inputs almost flawlessly regardless of the
region explored in the $N_{t},k$ domain (Fig. 3a).

With regards to the MAE, Janon/Monod, Jansen and Azzini and Rosati maintain
their high performance regardless of the $N_{t},k$ region explored. The
accuracy of Razavi and Gupta, however, drops at the upper-leftmost part of the
$N_{t},k$ boundary, where most of the largest MAE scores are located
($\mbox{MAE}>10$). In the case of Saltelli and Homma and Saltelli, the largest
MAE values concentrate in the region of small $k$ regardless of the total
number of model runs, a domain in which they achieved a high performance when
the focus was on properly ranking the model inputs.

Figure 3: Number of runs $N_{t}$ against the function dimensionality $k$. Each
dot is a simulation with a specific combination of the benchmark parameters in
Table 2. The greener (blacker) the color, the better (worse) the performance
of the estimator. a) Accuracy of the estimators when the goal is to properly
rank the model inputs, e.g. a factor prioritization setting. b) Accuracy of
the estimators when the goal is to approach the “true” total-order indices.

The presence of a non-negligible proportion of model runs with $r<0$ suggests
that some estimators significantly overturned the true ranks (Figs 3a, S4). To
better examine this phenomenon, we re-plot Fig 3b with just the simulations
yielding $r<0$ (Fig. S5). We observe that $r<0$ values not only appear in the
region of small $N_{t}$, a foreseeable miscalculation derived from allocating
an insufficient number of model runs to each model input: they also emerge at
a relatively large $N_{t}$ and low $k$ in the case of pseudo-Owen, Saltelli
and Homma and Saltelli. The Saltelli estimator actually concentrates in the
$k<10$ zone most of the simulations with the lowest negative $r$ values (Fig.
S5). This suggests that rank reversing is not an artifact of our study design
as much as a by-product of the volatility of these estimators when stressed by
the sources of computational uncertainty listed in Table 2. Such strain may
lead these estimators to produce a significant fraction of negative indices or
indices beyond 1, thus effectively promoting $r<0$.

We calculate the proportion of $T_{i}<0$ and $T_{i}>1$ in each simulation that
yielded $r<0$. In the case of Glen and Isaacs and Homma and Saltelli, $r<0$
values are caused by the production of a large proportion of $T_{i}<0$
(25%–75%, the $x$ axis in Fig. 4). Pseudo-Owen and Saltelli suffer this bias
too and in several simulations they also generate a large proportion of
$T_{i}>1$ (up to 100% of the model inputs, the $y$ axis in Fig. 4). The
production of $T_{i}<0$ and $T_{i}>1$ is caused by numerical errors and
fostered by the values generated at the numerator of Equation 3: $T_{i}<0$ may
either derive from $E_{\bm{x}_{\sim i}}\big{[}V_{x_{i}}(y|\bm{x}_{\sim
i})\big{]}<0$ (e.g. Homma and Saltelli and pseudo-Owen) or $V_{\bm{x}_{\sim
i}}\big{[}E_{x_{i}}(y|\bm{x}_{\sim i})\big{]}>V(y)$ (e.g. Saltelli), whereas
$T_{i}>1$ from $E_{\bm{x}_{\sim i}}\big{[}V_{x_{i}}(y|\bm{x}_{\sim
i})\big{]}>V(y)$ (e.g. Homma and Saltelli and pseudo-Owen) or $V_{\bm{x}_{\sim
i}}\big{[}E_{x_{i}}(y|\bm{x}_{\sim i})\big{]}<0$ (e.g. Saltelli).

Figure 4: Scatterplot of the proportion of $T_{i}<0$ against the proportion of
$T_{i}>1$ mapped against the model output $r$. Each dot is a simulation. Only
simulations with $r<0$ are displayed.

To better examine the efficiency of the estimators, we summarized their
performance as a function of the number of runs available per model input
$N_{t}/k$ [21] (Fig. 5, S6). This information is especially relevant to take
an educated decision on which estimator to use in a context of a high-
dimension, computationally expensive model. Even when the budget of runs per
input is low $\left[(N_{t}/k)\in[2,20]\right]$, Razavi and Gupta, Jansen and
Janon/Monod are very good at properly ranking model inputs ($r\approx 0.9$),
and are followed very close by Azzini and Rosati ($r\approx 0.8$). Saltelli,
Homma and Saltelli and Glen and Isaacs come after ($r\approx 0.3$), with
pseudo-Owen scoring last ($r\approx 0.2$). When the $N_{t}/k$ ratio is
increased, all estimators improve their ranking accuracy and some quickly
reach the asymptote: this is the case of Razavi and Gupta, Janon/Monod and
Jansen, whose performance becomes almost flawless from $(N_{t}/k)\in[40,60]$
onwards, and of Azzini and Rosati, which reaches its optimum at
$(N_{t}/k)\in[60,80]$. The accuracy of the other estimators does not seem to
fully stabilize within the range of ratios examined. In the case of Homma and
Saltelli and Saltelli, their performance oscillates before plummeting at
$(N_{t}/k)\in[200,210]$, $(N_{t}/k)\in[240,260]$ and $(N_{t}/k)\in[260,280]$
due to several simulations yielding large $r<0$ values (Fig. 5a).

Janon/Monod and Jansen are also the most efficient estimators when the MAE is
the measure of choice, followed closely by Azzini and Rosati, Razavi and Gupta
and Glen and Isaacs. Saltelli and Homma and Saltelli gain accuracy at higher
$N_{t}/k$ ratios yet their precision diminishes all the same from
$(N_{t}/k)\in[200,210]$ onwards (Fig. 5b).

Figure 5: Scatterplot of the model output $r$ against the number of model runs
allocated per model input $(Nt/k)$. See Fig. S6 for a visual display of all
simulations and Fig. S7 for an assessment of the number of model runs that
each estimator has in each $N_{t}/k$ compartment.

###  3.2 Sensitivity analysis

When the aim is to rank the model inputs, the selection of the performance
measure ($\delta$) has the highest first-order effect in the accuracy of the
estimators (Fig. 6a). The parameter $\delta$ is responsible for between 20%
(Azzini and Rosati) and 30% (Glen and Isaacs) of the variance in the final $r$
value. On average, all estimators perform better when the rank is conducted on
Savage scores ($\delta=2$), i.e. when the focus is on ranking the most
important model inputs only (Figs. S8–S15). As for the distribution of the
model inputs ($\phi$), it has a first-order effect in the accuracy of Azzini
and Rosati ($\approx 10$%), Jansen and Janon / Monod ($\approx 15$%) and
Razavi and Gupta ($\approx 20$%) regardless of whether the aim is a factor
prioritization ($r$) or approaching the “true” indices (MAE). The performance
of these estimators drops perceptibly when the model inputs are distributed as
$Beta(8,2)$ or $Beta(2,8)$ ($\phi=3$ and $\phi=4$, Figs. S8-S23), suggesting
that they may be especially stressed by skewed distributions. The selection of
random or quasi-random numbers during the construction of the sample matrix
($\tau$) also directly conditions the accuracy of several estimators. If the
aim is to approach the “true” indices (MAE), $\tau$ conveys from 17% (Azzini
and Rosati) to $\approx 30$% (Glen and Isaacs) of the model output variance,
with all estimators except Razavi and Gupta performing better on quasi-random
numbers ($\tau=2$, Figs. S16–S23). In a factor prioritization setting, $\tau$
is mostly influential through interactions. Interestingly, the proportion of
active second and third-order interactions ($k_{2},k_{3}$) does not alter the
performance of any estimator in any of the settings examined.

Figure 6: Sobol’ indices. a) Individual parameters. b) Clusters of parameters.
The cluster $f(x)$ includes all parameters that describe the uncertainty in
the functional form of the model ($\epsilon,k_{2},k_{3},\phi$). $N_{t}$ and
$k$ are assessed simultaneously due to their correlation. Note that the MAE
facet does not include the group ($\delta\tau$) because $\delta$ (the
performance measure used) is no longer an uncertain parameter in this setting.

To better understand the structure of the sensitivities, we compute Sobol’
indices after grouping individual parameters in three clusters, which we
define based on their commonalities: the first group includes $(\delta,\tau)$
and reflects the influence of those parameters that can be defined by the
sensitivity analyst during the setting of the benchmark exercise. The second
combines ($\varepsilon,k_{2},k_{3},\phi$) and examines the overall impact of
the model functional form, referred to as $f(x)$, which is often beyond the
analyst’s grasp. Finally, the third group includes $(N_{t},k)$ only and
assesses the influence of the sampling design in the accuracy of the
estimators (we assume that the total number of model runs, besides being
conditioned by the computing resources at hand, is also partially determined
by the joint effect of the model dimensionality and the use of either a
$\bm{B}$, $\bm{A}_{B}^{(i)})$, $\bm{B}_{A}^{(i)}$ or $\bm{C}_{B}^{(i)}$
matrices) (Fig 6b).

The uncertainty in the functional form of the model [$f(x)$] is responsible
for approximately 20% of the variance in the performance of Azzini and Rosati,
Janon/Monod or Jansen in a factor prioritization setting. For Glen and Isaacs,
Homma and Saltelli, pseudo-Owen or Saltelli, $f(x)$ is influential only
through interactions with the other clusters. When the MAE is the performance
measure of interest, $f(x)$ has a much stronger influence in the accuracy of
the estimators than the couple $(N_{t},k)$, especially in the case of Glen and
Isaacs ($\approx 40$%). In any case, the accuracy of the estimators is
significantly conditioned by interactions between the benchmark parameters.
The sum of all individual $S_{i}$ indices plus the $S_{i}$ index of the
$(N_{t},k)$ cluster only explains from $\approx 45$% (Saltelli) to $\approx
70$% (Glen and Isaacs) of the estimators’ variance in ranking the model
inputs, and from $\approx 24$% (pseudo-Owen) to $\approx 60$% (Razavi and
Gupta) of the variance in approaching the “true” indices.

##  4 Discussion and conclusions

Here we design an eight-dimension background for variance-based total-order
estimators to confront and prove their value in an unparalleled range of SA
scenarios. By randomizing the parameters that condition their performance, we
obtain a comprehensive picture of the advantages and disadvantages of each
estimator and identify which particular benchmark factors make them more prone
to error. Our work thus provides a thorough empirical assessment of state-of-
the-art total-order estimators and contributes to define best practices in
variance-based SA. The study also aligns with previous works focused on
testing the robustness of the tools available to sensitivity analysts, a line
of inquiry that can be described as a _sensitivity analysis of a sensitivity
analysis_ (SA of SA) [42].

Our results provide support to the assumption that the scope of previous
benchmark studies is limited by the plethora of non-unique choices taken
during the setting of the analysis [21]. We have observed that almost all
decisions have a non-negligible effect: from the selection of the sampling
method to the choice of the performance measure, the design prioritized by the
analyst can influence the performance of the estimator in a non-obvious way,
namely through interactions. The importance of non-additivities in
conditioning performance suggests that the benchmark of sensitivity estimators
should no longer rely on statistical designs that change one parameter at a
time (usually the number of model runs and, more rarely, the test function
[16, 33, 18, 23, 14, 20, 39, 40, 42]). Such setting reduces the uncertain
space to a minimum and misses the effects that the interactions between the
benchmark parameters have in the final accuracy of the estimator. If global SA
is the recommended practice to fully explore the uncertainty space of models,
sensitivity estimators, being algorithms themselves, should be likewise
validated [42].

Our approach also compensates the lack of studies on the theoretical
properties of estimators in the sensitivity analysis literature (see for
instance [47, 15]), and allows a more detailed examination of their
performance than theoretical comparisons. Empirical studies like ours mirror
the numerical character of sensitivity analysis when the indices can not be
analytically calculated, which is most of the time in “real-world”
mathematical modeling.

Two recommendations emerge from our work: the estimators by Razavi and Gupta,
Jansen, Janon / Monod or Azzini and Rosati should be preferred when the aim is
to rank the model inputs. Jansen, Janon/Monod or Azzini and Rosati should also
be prioritized if the goal is to estimate the “true” total-order indices. The
drop in performance of Razavi and Gupta in the second setting may be explained
by a bias at a lower sample sizes, i.e. a consistent over-estimation of all
total-order indices. This is because their estimator relies on a constant mean
assumption whose validity degrades with larger values of $\Delta h$ [20, 39].
In order to remove this bias, $\Delta h$ should take very small values (e.g.,
$\Delta h=0.01$), which may not be computationally feasible. Since the
direction of this bias is the same for all parameters it only affects the
calculation of the “true” total-order indices, not the capacity of the
estimator to properly rank the model inputs.

It is also worth stating that Razavi and Gupta is the only estimator studied
here that require the analyst to define a tuning parameter, $\Delta h$. In
this paper we have set $\Delta h=0.2$ after some preliminary trials with the
estimator; other works have used different values (e.g. $\Delta h=0.002$,
$\Delta h=0.1$, $\Delta h=0.3$; [21, 20, 39]). Selecting the most appropriate
value for a given tuning parameter is not an obvious choice and this
uncertainty can make an estimator volatile, as shown by [42] in the case of
the PAWN index.

The fact that Glen and Isaacs, Homma and Saltelli, Saltelli and pseudo-Owen do
not perform as well in properly ranking the model inputs and approaching the
“true” total-order indices may be partially explained by their less efficient
computation of elementary effects: by allowing the production of negative
terms in the numerator these estimators also permit the production of negative
total-order indices, thus leading to biased rankings or sensitivity indices.
In the case of Saltelli, the use of a $\bm{B}$ matrix at the numerator and an
$\bm{A}$ matrix at the denominator exacerbates its volatility (Table 1, Nº 5).
Such inconsistency was corrected in [23].

The consistent robustness of Jansen, Janon/Monod and Azzini and Rosati makes
their sensitivity to the uncertain parameters studied here almost negligible.
They are already highly optimized estimators with not much room for
improvement. Most of their performance is conditioned by the first and total-
order effects of the model form jointly with the underlying probability
distributions ($f(x)$ in Fig. 6b), as well as by their sampling design
($N_{t},k$), which are in any case beyond the analyst’s control. As for the
rest, their accuracy might be enhanced by allocating a larger number of model
runs per input (if computationally affordable), and especially in the case of
Homma and Saltelli, Saltelli and Glen and Isaacs, by restricting their use to
low-dimensional models ($k<10$) and sensitivity settings that only require
ranking the most important parameters (a _restricted_ factor prioritisation
setting; [1]). Nevertheless, their substantial volatility is considerably
driven by non-additivities, a combination that makes them hard to tame and
should raise caution about their use in any modeling exercise.

Our results slightly differ from [21]’s, who observed that Jansen outperformed
Janon/Monod under a factor prioritization setting. We did not find any
significant difference between these estimators. Although our metafunction
approach is based on [21]’s, our study tests the accuracy of estimators in a
larger uncertain space as we also account for the stress introduced by changes
in the sampling method $\tau$, the underlying probability distributions $\phi$
or the performance measure selected $\delta$. These differences may account
for the slightly different results obtained between the two papers.

Our analysis can be extended to other sensitivity estimators (i.e. moment-
independent like entropy-based [48]; the $\delta$-measure [49]; or the PAWN
index, [50, 51]). Moreover, it holds potential to be used overall as a
standard crash test every time a new sensitivity estimator is introduced to
the modeling community. One of its advantages is its flexibility: [21]’s
metafunction can be easily extended with new univariate functions or
probability distributions, and the settings modified to check performance
under different degrees of non-additivities or in a larger $(N_{t},k)$ space.
With some slight modifications it should also allow to produce functions with
dominant low-order or high-order terms, labeled as Type B and C by [22]. This
should prompt developers of sensitivity indices to severely stress their
estimators so the modeling community and decision-makers fully appraise how
they deal with uncertainties.

##  5 Code availability

The R code to replicate our results is available in [52] and in GitHub
(https://github.com/arnaldpuy/battle_estimators). The uncertainty and
sensitivity analysis have been carried out with the R package sensobol [53],
which also includes the test function used in this study.

##  6 Acknowledgements

We thank Saman Razavi for his insights on the behavior of the Razavi and Gupta
estimator. This work has been funded by the European Commission (Marie
Skłodowska-Curie Global Fellowship, grant number 792178 to A.P.).

## References

  * [1] Andrea Saltelli, Marco Ratto, Terry Andres, Francesca Campolongo, Jessica Cariboni, Debora Gatelli, Michaela Saisana and Stefano Tarantola “Global Sensitivity Analysis. The Primer” Chichester, UK: John Wiley & Sons, Ltd, 2008 DOI: 10.1002/9780470725184
  * [2] A.J. Jakeman, R.A. Letcher and J.P. Norton “Ten iterative steps in development and evaluation of environmental models” In _Environmental Modelling & Software_ 21.5, 2006, pp. 602–614 DOI: 10.1016/j.envsoft.2006.01.004
  * [3] Sibel Eker, Elena Rovenskaya, Michael Obersteiner and Simon Langan “Practice and perspectives in the validation of resource management models” In _Nature Communications_ 9.1 Springer US, 2018, pp. 1–10 DOI: 10.1038/s41467-018-07811-9
  * [4] A. Saltelli “Sensitivity analysis for importance assessment” In _Risk Analysis_ 22.3, 2002, pp. 579–590 DOI: 10.1111/0272-4332.00040
  * [5] B. Iooss and P. Lemaître “A review on global sensitivity analysis methods” In _Uncertainty Management in Simulation-Optimization of Complex Systems. Operations Research/Computer Science Interfaces Series, vol 59_ Boston: Springer, 2015, pp. 101–122 DOI: 10.1007/978-1-4899-7547-8˙5
  * [6] William Becker and Andrea Saltelli “Design for sensitivity analysis” In _Handbook of Design and Analysis of Experiments_ Boca Ratón: CRC Press, Taylor & Francis, 2015, pp. 627–674 DOI: 10.1201/b18619
  * [7] Toshimitsu Homma and Andrea Saltelli “Importance measures in global sensitivity analysis of nonlinear models” In _Reliability Engineering & System Safety_ 52.1, 1996, pp. 1–17 DOI: 10.1016/0951-8320(96)00002-6
  * [8] Loïc Le Gratiet, Stefano Marelli and Bruno Sudret “Metamodel-Based Sensitivity Analysis: Polynomial Chaos Expansions and Gaussian Processes” In _Handbook of Uncertainty Quantification_ Cham: Springer International Publishing, 2017, pp. 1289–1325 DOI: 10.1007/978-3-319-12385-1˙38
  * [9] Andrea Saltelli, S. Tarantola and K..-S. Chan “A quantitative model-independent method for global sensitivity analysis of model output” In _Technometrics_ 41.1, 1999, pp. 39 DOI: 10.2307/1270993
  * [10] R I Cukier, C M Fortuin, Kurt E Shuler, A G Petschek and J H Schaibly “Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I Theory” In _The Journal of Chemical Physics_ 59.8 American Institute of Physics, 1973, pp. 3873–3878
  * [11] R I Cukier, H B Levine and K E Shuler “Nonlinear sensitivity analysis of multiparameter model systems” In _Journal of Computational Physics_ 26.1 Elsevier, 1978, pp. 1–42
  * [12] Ilya.. Sobol’ “On the distribution of points in a cube and the approximate evaluation of integrals” In _USSR Computational Mathematics and Mathematical Physics_ 7.4, 1967, pp. 86–112 DOI: 10.1016/0041-5553(67)90144-9
  * [13] Ilya.. Sobol’ “Uniformly distributed sequences with an additional uniform property” In _USSR Computational Mathematics and Mathematical Physics_ 16.5, 1976, pp. 236–242 DOI: 10.1016/0041-5553(76)90154-3
  * [14] Samuele Lo Piano, Federico Ferretti, Arnald Puy, Daniel Albrecht and Andrea Saltelli “Variance-based sensitivity analysis: The quest for better estimators and designs between explorativity and economy” In _Reliability Engineering & System Safety_ 206 Elsevier Ltd, 2021, pp. 107300 DOI: 10.1016/j.ress.2020.107300
  * [15] M. Jansen “Analysis of variance designs for model output” In _Computer Physics Communications_ 117.1-2, 1999, pp. 35–43 DOI: 10.1016/S0010-4655(98)00154-4
  * [16] A. Janon, T. Klein, A. Lagnoux, M. Nodet and C. Prieur “Asymptotic normality and efficiency of two Sobol index estimators” In _ESAIM: Probability and Statistics_ 18.3, 2014, pp. 342–364 DOI: 10.1051/ps/2013040
  * [17] Graham Glen and Kristin Isaacs “Estimating Sobol sensitivity indices using correlations” In _Environmental Modelling and Software_ 37 Elsevier Ltd, 2012, pp. 157–166 DOI: 10.1016/j.envsoft.2012.03.014
  * [18] Ivano Azzini, Thierry Mara and Rossana Rosati “Monte Carlo estimators of first-and total-orders Sobol’ indices”, 2020 arXiv: http://arxiv.org/abs/2006.08232
  * [19] Hervé Monod, Cédric Naud and David Makowski “Uncertainty and Sensitivity Analysis for Crop Models” In _Working with Dynamic Crop Models_ Elsevier, 2006, pp. 35–100 URL: http://reseau-mexico.fr/sites/reseau-mexico.fr/files/06%7B%5C_%7Dzebook%7B%5C_%7DCH-03.pdf
  * [20] Saman Razavi and Hoshin V. Gupta “A new framework for comprehensive, robust, and efficient global sensitivity analysis: 2. Application” In _Water Resources Research_ 52.1, 2016, pp. 440–455 DOI: 10.1002/2015WR017558
  * [21] William Becker “Metafunctions for benchmarking in sensitivity analysis” In _Reliability Engineering and System Safety_ 204 Elsevier Ltd, 2020, pp. 107189 DOI: 10.1016/j.ress.2020.107189
  * [22] S. Kucherenko, B. Feil, N. Shah and W. Mauntz “The identification of model effective dimensions using global sensitivity analysis” In _Reliability Engineering & System Safety_ 96.4 Elsevier, 2011, pp. 440–449 DOI: 10.1016/j.ress.2010.11.003
  * [23] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto and S. Tarantola “Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index” In _Computer Physics Communications_ 181.2 Elsevier B.V., 2010, pp. 259–270 DOI: 10.1016/j.cpc.2009.09.018
  * [24] T. Ishigami and T. Homma “An importance quantification technique in uncertainty analysis for computer models” In _Proceedings. First International Symposium on Uncertainty Modeling and Analysis_ 12, 1990, pp. 398–403
  * [25] Ilya.. Sobol’ “On quasi-Monte Carlo integrations” In _Mathematics and Computers in Simulation_ 47.2-5, 1998, pp. 103–112 DOI: 10.1016/S0378-4754(98)00096-2
  * [26] P. Bratley and B.. Fox “ALGORITHM 659: implementing Sobol’s quasirandom sequence generator” In _ACM Transactions on Mathematical Software (TOMS)_ 14.1, 1988, pp. 88–100
  * [27] Emanuele Borgonovo and Elmar Plischke “Sensitivity analysis: A review of recent advances” In _European Journal of Operational Research_ 248.3 Elsevier Ltd., 2016, pp. 869–887 DOI: 10.1016/j.ejor.2015.06.032
  * [28] Andrea Saltelli “A short comment on statistical versus mathematical modelling” In _Nature Communications_ 10.1 Springer US, 2019, pp. 8–10 DOI: 10.1038/s41467-019-11865-8
  * [29] Andrea Saltelli, Gabriele Bammer, Isabelle Bruno, Erica Charters, Monica Di Fiore, Emmanuel Didier, Wendy Nelson Espeland, John Kay, Samuele Lo Piano, Deborah Mayo, Roger Pielke Jr, Tommaso Portaluri, Theodore M Porter, Arnald Puy, Ismael Rafols, Jerome R Ravetz, Erik Reinert, Daniel Sarewitz, Philip B Stark, Andrew Stirling, Jeroen Sluijs and Paolo Vineis “Five ways to ensure that models serve society: a manifesto” In _Nature_ 582.7813, 2020, pp. 482–484 DOI: 10.1038/d41586-020-01812-9
  * [30] R. Sheikholeslami, S. Razavi, H.. Gupta, W. Becker and A. Haghnegahdar “Global sensitivity analysis for high-dimensional problems: How to objectively group factors and measure robustness and convergence while reducing computational cost” In _Environmental Modelling and Software_ 111 Elsevier, 2019, pp. 282–299 DOI: 10.1016/j.envsoft.2018.09.002
  * [31] F. Sarrazin, F. Pianosi and T. Wagener “Global Sensitivity Analysis of environmental models: Convergence and validation” In _Environmental Modelling and Software_ 79 Elsevier Ltd, 2016, pp. 135–152 DOI: 10.1016/j.envsoft.2016.02.005
  * [32] Amin Haghnegahdar, Saman Razavi, Fuad Yassin and Howard Wheater “Multicriteria sensitivity analysis as a diagnostic tool for understanding model behaviour and characterizing model uncertainty” In _Hydrological Processes_ 31.25, 2017, pp. 4462–4476 DOI: 10.1002/hyp.11358
  * [33] Ivano Azzini and Rossana Rosati “The IA-Estimator for Sobol’ sensitivity indices” In _Ninth International Conference on Sensitivity Analysis of Model Output_ , 2019
  * [34] Mun Ju Shin, Joseph H A Guillaume, Barry F W Croke and Anthony J. Jakeman “Addressing ten questions about conceptual rainfall-runoff models with global sensitivity analyses in R” In _Journal of Hydrology_ 503 Elsevier B.V., 2013, pp. 135–152 DOI: 10.1016/j.jhydrol.2013.08.047
  * [35] L. Paleari and R. Confalonieri “Sensitivity analysis of a sensitivity analysis: We are likely overlooking the impact of distributional assumptions” In _Ecological Modelling_ 340 Elsevier B.V., 2016, pp. 57–63 DOI: 10.1016/j.ecolmodel.2016.09.008
  * [36] C. Spearman “The proof and measurement of association between two things” In _The American Journal of Psychology_ 15.1, 1904, pp. 72 DOI: 10.2307/1412159
  * [37] M.. Kendall and B. Smith “The problem of m rankings” In _The Annals of Mathematical Statistics_ 10.3, 1939, pp. 275–287 DOI: 10.1214/aoms/1177732186
  * [38] I. Savage “Contributions to the theory of rank order statistics - the two sample case” In _Annals of Mathematical Statistics_ 27, 1956, pp. 590–615
  * [39] Saman Razavi and Hoshin V. Gupta “A new framework for comprehensive, robust, and efficient global sensitivity analysis: 1. Theory” In _Water Resources Research_ 52.1, 2016, pp. 423–439 DOI: 10.1002/2015WR017559
  * [40] A.. Owen “Better estimation of small sobol’ sensitivity indices” In _ACM Transactions on Modeling and Computer Simulation_ 23.2, 2013, pp. 1–17 DOI: 10.1145/2457459.2457460
  * [41] Bertrand Iooss, Alexandre Janon, Gilles Pujol, with contributions from Baptiste Broto, Khalid Boumhaout, Sebastien Da Veiga, Thibault Delage, Reda El Amri, Jana Fruth, Laurent Gilquin, Joseph Guillaume, Loic Le Gratiet, Paul Lemaitre, Amandine Marrel, Anouar Meynaoui, Barry L Nelson, Filippo Monari, Roelof Oomen, Oldrich Rakovec, Bernardo Ramos, Olivier Roustant, Eunhye Song, Jeremy Staum, Roman Sueur, Taieb Touati and Frank Weber “sensitivity: Global Sensitivity Analysis of Model Outputs” R package version 1.22.1, 2020 URL: https://cran.r-project.org/package=sensitivity
  * [42] Arnald Puy, Samuele Lo Piano and Andrea Saltelli “A sensitivity analysis of the PAWN sensitivity index” In _Environmental Modelling and Software_ 127 Elsevier Ltd, 2020, pp. 104679 DOI: 10.1016/j.envsoft.2020.104679
  * [43] Art B. Owen “Randomly permuted (t, m, s)-nets and (t, s)-sequences” In _Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing. Lecture Notes in Statistics, vol. 106_ , 1995, pp. 299–317
  * [44] George E.. Box, J. Hunter and William G. Hunter “Statistics for Experimenters: Design, Innovation, and Discovery” Wiley, 2005
  * [45] Vilfredo Pareto “Manuale di Economia Politica. Vol. 13” Societa Editrice, 1906
  * [46] George E.. Box and R. Meyer “An analysis for unreplicated fractional factorials” In _Technometrics_ 28.1, 1986, pp. 11–18 DOI: 10.1080/00401706.1986.10488093
  * [47] I. Azzini, G. Listorti, T.. Mara and R. Rosati “Uncertainty and Sensitivity Analysis for policy decision making. An introductory guide” Luxembourg: Joint Research Centre, European Commission, 2020, pp. 1–66 DOI: 10.2760/922129
  * [48] Huibin Liu, Wei Chen and Agus Sudjianto “Relative entropy based method for probabilistic sensitivity analysis in engineering design” In _Journal of Mechanical Design, Transactions of the ASME_ 128.2, 2006, pp. 326–336 DOI: 10.1115/1.2159025
  * [49] E. Borgonovo “A new uncertainty importance measure” In _Reliability Engineering and System Safety_ 92.6, 2007, pp. 771–784 DOI: 10.1016/j.ress.2006.04.015
  * [50] F. Pianosi and T. Wagener “A simple and efficient method for global sensitivity analysis based on cumulative distribution functions” In _Environmental Modelling and Software_ 67 Elsevier LtdX_i=, 2015, pp. 1–11 DOI: 10.1016/j.envsoft.2015.01.004
  * [51] F. Pianosi and T. Wagener “Distribution-based sensitivity analysis from a generic input-output sample” In _Environmental Modelling and Software_ 108 Elsevier, 2018, pp. 197–207 DOI: 10.1016/j.envsoft.2018.07.019
  * [52] Arnald Puy “R code of ”A comprehensive comparison of total-order estimators for global sensitivity analysis”” In _Zenodo_ 4946559 Zenodo, 2021 DOI: 10.5281/zenodo.4946559
  * [53] Arnald Puy, Samuele Lo Piano, Andrea Saltelli and Simon A. Levin “sensobol: an R package to compute variance-based sensitivity indices” In _Journal of Statistical Software_ , 2021 arXiv: http://arxiv.org/abs/2101.10103