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

HAL
open science

# Making best use of permutations to compute sensitivity indices with replicated orthogonal arrays

Laurent Gilquin, Élise Arnaud, Clémentine Prieur, Alexandre Janon

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

## ► To cite this version:

Laurent Gilquin, Élise Arnaud, Clémentine Prieur, Alexandre Janon. Making best use of permutations to compute sensitivity indices with replicated orthogonal arrays. Reliability Engineering and System Safety, 2019, 187, pp.28-39. 10.1016/j.ress.2018.09.010. hal-01558915v3

HAL Id: hal-01558915

https://inria.hal.science/hal-01558915v3

Submitted on 5 Sep 2019

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Making the best use of permutations to compute
sensitivity indices with replicated orthogonal arrays

Laurent Gilquin

Inria Grenoble - Rhône-Alpes, LJK, 655 avenue de l'Europe, 38330 Montbonnot, France

Elise Arnaud

Université Grenoble Alpes, CNRS, LJK, F-38000 Grenoble, France

Inria project/team AIRSEA

Clémentine Prieur

Université Grenoble Alpes, CNRS, LJK, F-38000 Grenoble, France

Inria project/team AIRSEA

Alexandre Janon

Laboratoire de Mathématiques dOrsay, Université Paris-Sud, CNRS, Université Paris-Saclay, 91405 Orsay, FRANCE

Abstract

Among practitioners, the importance of inputs to a model output is com-
monly measured via the computation of Sobol' sensitivity indices. Various es-
timation strategies exist in the literature, most of them requiring a very high
number of model evaluations. Designing methods that compete favorably both
in terms of computational cost and accuracy is therefore an issue of crucial
importance. In this paper, an efficient replication-based strategy is proposed
to estimate the full set of first- and second-order Sobol' indices. It relies on a
Sobol' pick-freeze estimation scheme and requires only two replicated designs
based on randomized orthogonal arrays of strength two. The precision of this
procedure is assessed with bootstrap confidence intervals, presented for the first
time in the replication framework. Our developments are compared to known
approaches and validated on numerical test cases. A way to estimate the full set
of first-, second-order but also total-effect Sobol' indices at a very competitive
cost is also described, as a combination of our procedure and the one introduced
by Saltelli in [1].

Keywords: Bootstrap confidence intervals, Computer experiments, Replicated designs, Sensitivity Analysis, Sobol' indices

# 1. Introduction

Many mathematical models encountered in applied sciences involve numer-
ous poorly-known inputs. It is important for the practitioner to understand
how the output uncertainty can be apportioned to the uncertainty in the in-
puts. One way to do so is to perform a global sensitivity analysis in which
statistical methods allow one to calculate importance measures (see, for exam-
ple, [2] and references therein). In this framework, the input vector is treated
as a random vector, whose joint probability distribution reflects the modeler's
knowledge about the inputs uncertainty. This turns the output into a random
variable. In this paper, it is assumed that the model output has a finite vari-
ance. When inputs are mutually independent, the estimation of Sobol' indices
introduced in [3] provides a way to measure the importance of individual (or set
of) component(s) of the input random vector. The derivation of these indices
is based on a functional analysis of variance (ANOVA) decomposition of the
model (for details and references on functional ANOVA, see [4]).

The Sobol' index can be loosely defined as the Pearson correlation coeffi-
cient of two output vectors evaluated at two random input vectors sharing an
identical subset of their components ([5, 6]). When the subset contains only
one component, a first-order Sobol' index is estimated, if the subset contains
two elements, then a closed second-order Sobol' index is calculated and so forth.
Estimators based on such a strategy have been introduced in [3], and are known
as pick-freeze estimators (see, for example, [7]).

This paper focuses on the Sobol' pick-freeze (SPF) estimation procedure.
This procedure can be applied with various sampling strategies to obtain the
pairs of output vector. The SPF scheme can be used to estimate any Sobol'
indices. However, it suffers from a number of model evaluations that is propor-
tional to the dimension *d* of the input vector. In typical engineering applications,
this dimension can be prohibitive. Indeed, the computation time required for a
single model evaluation can drastically restrict the total number of affordable
model runs.

30

This motivates the use of *replicated* designs for the estimation of Sobol' indices. Replicated designs based on Latin hypercube sampling (LHS) have been introduced in [8]. They are used in [9, 10] in which an arbitrary number *r* of replicated LHS is used to define an estimator of first-order indices, that does not belong to the SPF family. Numerical simulations reveal that this estimator suffers from non negligible biases. Therefore, a new sampling strategy is proposed in [10], constraining the concatenation of the *r* replicated LHS to be an orthogonal array of strength two. This improved strategy estimates efficiently first-order indices. Replicated designs have been used in an SPF scheme in [11], in which the authors suggest the use of only two replicated LHS to estimate the overall set of first-order Sobol' indices. Asymptotic properties of this so-called replication procedure have been studied in [12]. The authors in [12] have also extended the replication procedure to the estimation of closed second-order Sobol' indices using two orthogonal arrays.

35

40

The first aim of this paper is to improve the efficiency of the replication
procedure proposed in [12] for the estimation of first- and second-order indices.
Improving the efficiency has to be understood here as decreasing the number
of model evaluations required to reach a given accuracy. The latter accuracy
is measured by confidence intervals whose computation constitutes the second
objective of this paper. Indeed, when performing a global sensitivity analysis,
practitioners are not only interested in a point estimation of each Sobol' index,
but also in assessing the accuracy of the estimates. While asymptotic confidence
intervals (derived from the central limit theorem) for the replication procedure
are detailed in [12], and bootstrap procedure for classical Sobol' index estima-
tion in [13], bootstrap confidence intervals for the replication framework are
presented here for the first time, to the best of our knowledge.

45

50

55
The paper is organized as follows: in Section 2, the definitions of Sobol' indices and of their pick-freeze estimators are recalled. Section 3 is devoted to the replication procedure. First, the classical version [12] is presented. Secondly, our procedure is introduced, that uses only two replicated designs based on randomized orthogonal arrays of strength two to estimate the full set of first- and second-order Sobol' indices. Then, a way to estimate the full set of first-, second-order but also total-effect Sobol' indices at a very competitive cost is also described, as a combination of our procedure with the one described in [1]. The end of this section provides formula of the bootstrap confidence intervals. Section 4 presents numerical simulations and highlights the efficiency of our procedures through comparisons with other methods, in particular with [12], and through the study of a high-dimensional problem ($d = 50$). In parallel, our bootstrap confidence intervals are compared to asymptotic ones for the estimation of first-order indices.

## 2. Background on Sobol' indices

### 2.1. Definition of Sobol' indices

Denote by $\boldsymbol{x} = (x_1, \dots, x_d)$ the vector of inputs of a model response $f$. We assume without loss of generality that $\boldsymbol{x}$ is a random vector uniformly distributed over the unit hypercube $\mathbb{H}_d = [0, 1]^d$. We further assume that $f$ is square integrable and denote by $\mathcal{D} = \{1, \dots, d\}$ the set of indices. Let $\boldsymbol{u}$ be a non-empty subset of $\mathcal{D}$, $-\boldsymbol{u}$ its complement and $|\boldsymbol{u}|$ its cardinality. Then, $\boldsymbol{x}_\boldsymbol{u}$ represents a point in $\mathbb{H}_{|\boldsymbol{u}|}$ with components $x_j, j \in \boldsymbol{u}$. Given two points $\boldsymbol{x}$ and $\boldsymbol{x}'$, the hybrid point $(\boldsymbol{x}_\boldsymbol{u} : \boldsymbol{x}'_{-\boldsymbol{u}})$ is defined as $x_j$ if $j \in \boldsymbol{u}$ and $x'_j$ if $j \notin \boldsymbol{u}$.

Consider the Sobol-Hoeffding decomposition [14, 3, 15] of $f$:

$$f(\boldsymbol{x}) = f_{\emptyset} + \sum_{\boldsymbol{u} \subseteq \mathcal{D}} f_{\boldsymbol{u}}(\boldsymbol{x}_{\boldsymbol{u}}). \qquad (1)$$

Under the assumptions that $f_{\emptyset}$ is constant and each term $f_{\boldsymbol{u}}$ satisfies:

$$\int_{\mathbb{H}} f_{\boldsymbol{u}}(\boldsymbol{x}_{\boldsymbol{u}})dx_j = 0, \quad \forall j \in \boldsymbol{u},$$
the Hoeffding decomposition becomes unique and,

$$f_{\emptyset} = \mathbb{E}[f(\boldsymbol{x})] = \mu, \\ f_u(\boldsymbol{x}_u) = \int_{\mathbb{H}_{d-|u|}} f(\boldsymbol{x})d\boldsymbol{x}_{-u} - \sum_{\boldsymbol{v} \preceq \boldsymbol{u}} f_v(\boldsymbol{x}_v).$$

As a consequence, taking the variance of both sides of Eq. (1) leads to the variance decomposition of $f$:

$$\sigma^2 = \mathrm{Var}[f(\boldsymbol{x})] = \sum_{\boldsymbol{u} \subseteq \mathcal{D}, \boldsymbol{u} \neq \emptyset} \sigma_{\boldsymbol{u}}^2, \quad \text{where } \sigma_{\boldsymbol{u}}^2 = \int_{\mathbb{H}_{|\boldsymbol{u}|}} f_{\boldsymbol{u}}(\boldsymbol{x}_{\boldsymbol{u}})^2 d\boldsymbol{x}_{\boldsymbol{u}}.$$

From this latter decomposition, one can define the following two quantities:

$$\underline{\tau}_{\boldsymbol{u}}^2 = \sum_{\boldsymbol{v} \subseteq \boldsymbol{u}} \sigma_{\boldsymbol{v}}^2, \quad \bar{\tau}_{\boldsymbol{u}}^2 = \sum_{\boldsymbol{v} \cap \boldsymbol{u} \neq \emptyset} \sigma_{\boldsymbol{v}}^2, \quad \boldsymbol{u} \subseteq \mathcal{D}.$$

These two quantities $\underline{\tau}_{\boldsymbol{u}}^2$ and $\bar{\tau}_{\boldsymbol{u}}^2$ measure the importance of variables $\boldsymbol{x}_{\boldsymbol{u}}$: $\underline{\tau}_{\boldsymbol{u}}^2$ quantifies the main effect of $\boldsymbol{x}_{\boldsymbol{u}}$, that is the effect of all the variables within $\boldsymbol{x}_{\boldsymbol{u}}$ including interactions, and $\bar{\tau}_{\boldsymbol{u}}^2$ quantifies the main effect of $\boldsymbol{x}_{\boldsymbol{u}}$ plus the effect of all interactions between variables in $\boldsymbol{x}_{\boldsymbol{u}}$ and variables in $\boldsymbol{x}_{-\boldsymbol{u}}$.

Therefore, we have the relation $0 \le \underline{\tau}_{\boldsymbol{u}}^2 \le \bar{\tau}_{\boldsymbol{u}}^2 \le 1$. These two measures are commonly found in the literature in their normalized form: $\underline{S}_{\boldsymbol{u}} = \underline{\tau}_{\boldsymbol{u}}^2/\sigma^2$ is the closed $|\boldsymbol{u}|^{\text{th}}$-order Sobol' index for inputs $\boldsymbol{x}_{\boldsymbol{u}}$, while $\bar{S}_{\boldsymbol{u}} = \bar{\tau}_{\boldsymbol{u}}^2/\sigma^2$ is the total-effect Sobol' index of order $|\boldsymbol{u}|$. In this paper, we focus on first-order and total-effect Sobol' indices, namely $S_j = \underline{S}_{\{j\}}$ and $\bar{S}_{\{j\}}$, $j \in \mathcal{D}$. We also consider second-order Sobol' indices $S_{k,l} = \underline{S}_{\{k,l\}} - \underline{S}_{\{k\}} - \bar{S}_{\{l\}}$, $\{k,l\} \in \mathcal{D}^2; k \ne l$.

The computation of the normalized indices is performed based on the following integral formulas:

$$\underline{\tau}_{\boldsymbol{u}}^2 = \int_{\mathbb{H}_{d-|\boldsymbol{u}|}} f(\boldsymbol{x}_{\boldsymbol{u}} : \boldsymbol{x}'_{-\boldsymbol{u}}) f(\boldsymbol{x}) d\boldsymbol{x} d\boldsymbol{x}' - \mu^2, \quad (2)$$

$$\sigma^2 = \int_{\mathbb{H}_d} f(\boldsymbol{x})^2 d\boldsymbol{x} - \mu^2, \\ \mu = \int_{\mathbb{H}_d} f(\boldsymbol{x}) d\boldsymbol{x}. \quad (3)$$
According to the law of total variance, the numerator of the total-effect Sobol' index can be written:

$$ \bar{\tau}_{\{j\}}^2 = \sigma^2 - \underline{\tau}_{\{-j\}}^2. \qquad (4) $$

Usually the complexity of $f$ causes the solution of integrals (2 - 4) to be analytically intractable. In such cases, one can instead estimate these quantities.

## 2.2. Estimation of Sobol' indices

While estimating Sobol' indices, one has to choose both an estimator and a design of experiments (simply called design hereafter). In this paper, we focus on the estimator introduced in [16], as it was proven in [7] to have optimal asymptotic variance properties. A design is a point set $\mathcal{P} = \{x_i\}_{i=1}^n$ in which each point is obtained by sampling $n$ times each input variable $x_j$, $j = 1, \dots, d$. Each row of the design is a point $x_i$ in $\mathbb{H}_d$ and the j-th column of the design refers to a sample of $x_j$. For $u \subset \mathcal{D}$, $x_{i,u}$ is a point in $\mathbb{H}_{|u|}$ with components $x_{i,j}, j \in u$. Consider $\mathcal{P} = \{x_i\}_{i=1}^n$ and $\mathcal{P}' = \{x'_i\}_{i=1}^n$ two designs uniformly distributed over $\mathbb{H}_d$. One way to estimate the quantity in (2) and (4) is as follows,

$$ \hat{\tau}_{u}^{2} = \frac{1}{n} \sum_{i=1}^{n} f(x_i) f(x_{i,u} : x'_{i,-u}) - \hat{\mu}^{2}, \qquad (5) $$

with:

$$ \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^{n} \frac{f^2(x_i) + f^2(x_{i,u} : x'_{i,-u})}{2} - \hat{\mu}^2, \qquad (6) $$

$$ \text{and } \hat{\mu} = \frac{1}{n} \sum_{i=1}^{n} \frac{f(x_i) + f(x_{i,u} : x'_{i,-u})}{2}. $$

The Sobol' indices estimators are then:

$$ \hat{S}_u = \hat{\tau}_u^2 / \hat{\sigma}^2, \quad u \subset \mathcal{D}, \qquad (7) $$

and

$$ \hat{S}_u = 1 - \frac{\frac{1}{n} \sum_{i=1}^{n} f(x'_i) f(x_{i,u} : x'_{i,-u}) - \hat{\mu}^2}{\hat{\sigma}^2}. \qquad (8) $$
By using Eqs. (5 - 8) and the sampling strategy proposed in [1], one can compute the overall sets of first-order and total-effect Sobol' indices with only $d+2$ designs, namely the designs $\{x_{i,u} : x'_{i,-u}\}_{i=1}^n$ constructed for $u \in \emptyset, \{1\}, \dots, \{d\}, \mathcal{D}$.$^{105}$ This is a fairly common strategy in variance-based sensitivity analysis and it is referred to as SAL02 in the present work.

While ingenious, this last approach requires a number of model evaluations that grows linearly with respect to the input space dimension, which may be unaffordable for some real applications with limited budget. A possible solution to this issue lies in the use of the replication procedure (see [17] for other alternatives involving additional regularity properties). The replication procedure allows one to estimate the full set of first- and second-order effect indices, at a cost independent from the input space dimension. The replication procedure is detailed in the next section.

### 3. Replication procedure

#### 3.1. TIS15, the current best replication procedure

The replication procedure has been introduced in [11] with the aim of evaluating all first-order indices at a reduced cost, namely $2n$ model evaluations. This procedure relies on the construction of two replicated designs. The notion of replicated designs was first introduced in [8] through the notion of replicated Latin Hypercubes. To extend this definition to other types of designs, one may use the generalization from [18]:

**Definition 1.** Let $\mathcal{P} = \{x_i\}_{i=1}^n$ and $\mathcal{P}' = \{x'_i\}_{i=1}^n$ be two non-identical designs in $\mathbb{H}_d$. $\mathcal{P}$ and $\mathcal{P}'$ are two replicated designs of order $p$, if for any $u \subset \mathcal{D}$ such that $|u| = p$, there exists a permutation $\pi_u$ of $\{1, \dots, n\}$ such that $\forall i \in \{1, \dots, n\}, x_{i,u} = x'_{\pi_u(i),u}$.

The procedure introduced in [11] allows one to estimate the full set of first-order Sobol' indices with only two replicated designs of order 1. This procedure has been deeply studied and generalized in Tissot and Prieur [12] to the estimation of the full set of closed second-order indices.
Let $\mathcal{P} = \{\mathbf{x}_i\}_{i=1}^n$ and $\mathcal{P}' = \{\mathbf{x}'_i\}_{i=1}^n$ be two replicated designs of order $|\mathbf{u}|$. The key point of the replication procedure is to use the permutation $\pi_{\mathbf{u}}$ defined in Definition 1 to produce the hybrid design $\{\mathbf{x}_{i,\mathbf{u}} : \mathbf{x}'_{i,-\mathbf{u}}\}_{i=1}^n$ used in Eq. (5). Denote by $\{y_i\}_{i=1}^n = \{f(\mathbf{x}_i)\}_{i=1}^n$ and $\{y'_i\}_{i=1}^n = \{f(\mathbf{x}'_i)\}_{i=1}^n$ the two sets of model responses obtained with $\mathcal{P}$ and $\mathcal{P}'$ respectively. From Definition 1, it results that,

$$ 
\begin{aligned}
y'_{\pi_{\mathbf{u}}(i)} &= f(\mathbf{x}'_{\pi_{\mathbf{u}}(i),\mathbf{u}} : \mathbf{x}'_{\pi_{\mathbf{u}}(i),-\mathbf{u}}), \\
&= f(\mathbf{x}_{i,\mathbf{u}} : \mathbf{x}'_{\pi_{\mathbf{u}}(i),-\mathbf{u}}).
\end{aligned}
$$

Hence, each $\underline{\tau}_{\mathbf{u}}^2$ can be estimated via formula (5) by using $y'_{\pi_{\mathbf{u}}(i)}$ instead of $f(\mathbf{x}_{i,\mathbf{u}} : \mathbf{x}'_{i,-\mathbf{u}})$ without requiring further model responses but $\{y_i\}_{i=1}^n$ and $\{y'_i\}_{i=1}^n$.

**Estimation of first-order indices.** Consider first the case $|\mathbf{u}| = 1$. As noted hereinbefore, the full set of first-order Sobol' indices can be estimated from two replicated designs of order 1, e.g. two replicated Latin hypercubes of size $n$.

**Estimation of closed second-order indices.** Consider now the case $|\mathbf{u}| = 2$. The full set of closed second-order Sobol' indices can be estimated from two replicated designs of order 2. Such designs are, for instance, designs based on orthogonal arrays of strength 2 (see [12] for more details). The structure of orthogonal arrays has been introduced by Kishen [19] and further extended by Rao [20]. It is defined as follows:

**Definition 2.** A $t - (q, d, \lambda)$ orthogonal array ($t \le d$) is a $\lambda q^t \times d$ matrix whose entries are chosen from a $q$-set of $N$ such that in every subset of $t$ columns of the array, every $t$-subset of points of this $q$-set appears in exactly $\lambda$ rows.

From this definition, one can construct, by setting $t = 2$, a structure consisting of points in $\{1, \dots, q\}^{\lambda q^2}$ in which each 2-set of columns have the same 2-set of points $\lambda$ times.

The so-called method of differences [21] is used to construct $2 - (d, q, 1)$ orthogonal arrays, with $q$ a prime number greater or equal to $d - 1$. This last
constraint is not particularly restrictive, as discussed in [12, 22], and can be relaxed when using other constructions (see, for example, [23]). This design is called an orthogonal array of strength two. Note that the orthogonal array should be transformed to fit the desired probability density function of each input. This requires two successive transformations. First, the values {1, ..., q} of the orthogonal array are normalized between [0, 1]. Then, each column of the normalized orthogonal array are mapped back with the corresponding inverse cumulative distribution function of each input.

With this construction, the replication procedure (called TIS15 in the following) allows one to estimate the full set of first- and closed second-order indices at the cost of $2n + 2q^2$, which depends on the input space dimension only through the constraint on $q$ (a prime number satisfying $q \ge d-1$). Note that, as soon as $n > 2(d-1)$, the number of model evaluations required for implementing procedure TIS15 is smaller than the one required to implement procedure SAL02. Thus, for time-consuming applications with limited budget, TIS15 should be preferred even if it does not allow to estimate total-effect Sobol' indices. Indeed, in many applications main effects and interactions of order 2 are often leading the sensitivity of the model, higher-order interactions being negligible.

## 3.2. REP18A, a new efficient strategy

This section introduces a more efficient strategy that estimates the full set of first- and second-order Sobol' indices from two replicated designs of order 2. This estimation strategy is referred to as REP18A and is explained below.

Consider two replicated designs of order 2: $\mathcal{P} = \{x_i\}_{i=1}^n$ and $\mathcal{P}' = \{x'_i\}_{i=1}^n$ each based on a $2 - (d, q, 1)$ orthogonal array of size $n = q^2$ (see [12] for more details). REP18A stems from the following observations:

(a) for any subset $u$ of $\mathcal{D}$ such that $|u| = 2$, there exists a unique permutation $\pi_u$ satisfying:

$$x'_{\pi_u(i),u} = x_{i,u}, \quad \forall i \in \{1, \dots, q^2\};$$
(b) for any $k \in \mathcal{D}$, there exists a set $\mathcal{H}_k = \{\pi_k^{(\gamma)}, \gamma \in \{1, \dots, (q!)^q\}\}$ of $(q!)^q$ mappings (or permutations) $\pi_k^{(\gamma)}$, satisfying:

$$ \mathbf{x}'_{\pi_k^{(\gamma)}(i),k} = \mathbf{x}_{i,k}, \quad \forall \gamma \in \{1, \dots, (q!)^q\}, \forall i \in \{1, \dots, q^2\}. $$

Based on these observations, it is possible to estimate the full set of first- and second-order Sobol' indices from the single pair of designs $(\mathcal{P}, \mathcal{P}')$. First, for any $\mathbf{u}$ subset of $\mathcal{D}$ with cardinality 2, for any $i \in \{1, \dots, q^2\}$, let

$$ y_i = f(\mathbf{x}_i), $$

$$ y'_i = f(\mathbf{x}'_i), $$

$$ \begin{aligned} y'_{\pi_u(i)} &= f(\mathbf{x}'_{\pi_u(i),u} : \mathbf{x}'_{\pi_u(i),-u}), \\ &= f(\mathbf{x}_{i,u} : \mathbf{x}'_{\pi_u(i),-u}). \end{aligned} $$

Additionally, for any $k \in \mathcal{D}$, for any $\gamma \in \{1, \dots, (q!)^q\}$, let

$$ \begin{aligned} y'_{\pi_k^{(\gamma)}(i)} &= f(\mathbf{x}'_{\pi_k^{(\gamma)}(i),k} : \mathbf{x}'_{\pi_k^{(\gamma)}(i),-k}), \\ &= f(\mathbf{x}_{i,k} : \mathbf{x}'_{\pi_k^{(\gamma)}(i),-k}). \end{aligned} $$

Note that the quantities $\{y'_{\pi_u(i)}\}_{i=1}^{q^2}$ and $\{y'_{\pi_k^{(\gamma)}(i)}\}_{i=1}^{q^2}$ require no additional evaluations of the model $f$, since the permutations $\pi_u$ and $\pi_k^{(\gamma)}$ act as permutations of rows on the design $\mathcal{P}'$.

Then, the main results of REP18A are outlined in Proposition 1.

**Proposition 1.** Let $\mathcal{P} = \{\mathbf{x}_i\}_{i=1}^{n=q^2}$ and $\mathcal{P}' = \{\mathbf{x}'_i\}_{i=1}^{n=q^2}$ be two replicated designs based on 2 – (d, q, 1) orthogonal arrays as in [12], where q is a prime number satisfying the constraint $q \ge d-1$. From model evaluations on $\mathcal{P} = \{\mathbf{x}_i\}_{i=1}^{q^2}$ and $\mathcal{P}' = \{\mathbf{x}'_i\}_{i=1}^{q^2}$, on can obtain:

* a single estimate for each $\underline{S}_{\mathbf{u}}$, $|\mathbf{u}| = 2$,

* $(q!)^q$ estimates for each $S_k$, $k \in \mathcal{D}$.

*Proof of Proposition 1.* Each $\underline{S}_{\mathbf{u}}$, $|\mathbf{u}| = 2$, can be estimated via formula (5, 6) by using $y'_{\pi_u(i)}$ in place of $f(\mathbf{x}_{i,u} : \mathbf{x}'_{i,-u})$. Each $S_k$, $k \in \mathcal{D}$, can be estimated via formula (5, 6) by using $y'_{\pi_k^{(\gamma)}(i)}$ in place of $f(\mathbf{x}_{i,k} : \mathbf{x}'_{i,-k})$ and there exists $(q!)^q$ choices for $\pi_k^{(\gamma)}$. $\square$
*Example.* Proposition 1 is illustrated on a simple example. Consider two replicated orthogonal arrays of strength two with $q = 3$ and $d = 4$:

$$ \mathcal{P} = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 2 & 1 & 3 & 2 \\ 3 & 1 & 2 & 3 \\ 2 & 2 & 2 & 1 \\ 3 & 2 & 1 & 2 \\ 1 & 2 & 3 & 3 \\ 3 & 3 & 3 & 1 \\ 1 & 3 & 2 & 2 \\ 2 & 3 & 1 & 3 \end{pmatrix} \qquad \mathcal{P}' = \begin{pmatrix} 2 & 1 & 2 & 1 \\ 3 & 3 & 2 & 2 \\ 2 & 3 & 3 & 3 \\ 3 & 2 & 3 & 1 \\ 1 & 1 & 3 & 2 \\ 3 & 1 & 1 & 3 \\ 1 & 3 & 1 & 1 \\ 2 & 2 & 1 & 2 \\ 1 & 2 & 2 & 3 \end{pmatrix} $$

Each column contains the set of values {1, 2, 3}, each value being repeated three times. Let us focus on the estimation of $S_{\{1,2\}}$. Item (a) in Section 3.2 states that there exists a unique permutation that re-arranges the first two columns of $\mathcal{P}'$ identically to the two first columns of $\mathcal{P}$. In our example, this permutation is $\pi_{\{1,2\}} = (5, 1, 6, 8, 4, 9, 2, 7, 3)$.

Let us now focus on the estimation of $S_1$. According to Item (b), there exist $(3!)^3$ permutations $\pi_1^{(\gamma)}$ that re-arrange the rows of the first column of $\mathcal{P}$ identically to the first column of $\mathcal{P}'$. Let denote by $\mathcal{H}_1 = \{\pi_1^{(\gamma)}, \gamma \in \{1, \dots, (3!)^3\}\}$ the set of possible re-arrangements. The following permutation $\pi_1^{(1)}$ is one out of the 81 possible re-arrangements: $\pi_1^{(1)} = (5, 1, 2, 3, 4, 5, 6, 9, 8)$.

Therefore, for each input $x_k$, one may pick one element of $\mathcal{H}_k$ to estimate the first-order index $S_k$. By making use of multiple mappings $\pi_k^{(\gamma)}$, one can enhance the estimation of the first-order indices.

*Remark.* For each $S_k$, several estimations can be obtained, one for each permutation in $\mathcal{H}_k$. It is rather natural to consider the mean of these estimates to increase the accuracy of the procedure. In practice, one should avoid computing the full set of $(q!)^q$ permutations $\pi_k^{(\gamma)}$ and rather select randomly a predetermined amount $\kappa$. In the numerical experiments, $\kappa$ has been chosen equal to 100 which is enough to yield satisfactory results. An optimal choice for $\kappa$, and an
optimal selection procedure for the $\kappa$ permutations in $\mathcal{H}_k$ is out of the scope of that paper, but is an interesting perspective.

For any $k \in \mathcal{D}$, denote by $\pi_k^{(\gamma_1)}, \dots, \pi_k^{(\gamma_\kappa)}$ the $\kappa$ permutations sampled from $\mathcal{H}_k$. $S_k$ is then estimated by:

$$ \hat{S}_k^{(\gamma_1, \dots, \gamma_\kappa)} = \frac{1}{\kappa} \sum_{s=1}^{\kappa} \hat{S}_k^{(\gamma_s)}, \qquad (9) $$

where $\hat{S}_k^{(\gamma_s)}$ is the estimator associated with $\pi_k^{(\gamma_s)}$. Then:

$$ \hat{S}_{\{k, \ell\}}^{(\gamma_1, \dots, \gamma_\kappa)} = \underline{\hat{S}}_{\{k, \ell\}} - \hat{S}_k^{(\gamma_1, \dots, \gamma_\kappa)} - \hat{S}_\ell^{(\gamma_1, \dots, \gamma_\kappa)} \qquad (10) $$

is the estimate of the second-order Sobol' index $S_{\{k, \ell\}}$.

### 3.3. REP18B, an extension to total-effect indices

Procedure REP18A does not provide estimates of total-effect Sobol' indices. This is not convenient when interactions of order higher than two prevail. REP18A can be adapted so as to incorporate the estimation of total-effect indices. This extension, named REP18B afterward, exploits the sampling strategy SAL02 [1].

REP18B proceeds as follows: first, all first-order and second-order indices are estimated with two replicated designs (see Section 3.1). Then, all total-effect Sobol' indices are estimated using one of the two replicated designs plus $d$ additional designs obtained following Saltelli's scheme [1]. More specifically, let $\mathcal{P}$ and $\mathcal{P}'$ denote the two replicated orthogonal arrays. The $j$-th additional design is constructed by substituting the $j$-th column of $\mathcal{P}$ for a $j$-th MC sample of equal length. Then, the $j$-th additional design and $\mathcal{P}$ are used together to estimate $\bar{S}_{\{j\}}$ according to Saltelli's scheme.

For $n = q^2$, the cost of REP18B reads $q^2(d + 2)$, which equals SAL02 cost. Note that REP18B outperforms SAL02 in the sense that it also provides second-order indices. However, due to constraints on the construction of strength two OAs, one has to select $q \ge d - 1$. If strong interactions are expected, of order greater or equal to 3, one should prefer SAL02 for models with high dimension $d$.
### 3.4. Bootstrap intervals

Sampling error in the estimation of Sobol' indices can be classically estimated, at a moderate cost, by using bootstrap resamplings such as the bias-corrected (BC) percentile method [24, 25]. This method is reviewed hereafter.

Let $\hat{\theta}_n(w)$ be an estimator of an unknown statistic $\theta$ based on $n$ realizations $w_1, \dots, w_n$ of a random vector $w$. Then, $B$ subsets $\{w_1[b], \dots, w_n[b]\}$, $b = 1, \dots, B$, are randomly drawn from $\{w_1, \dots, w_n\}$ with replacement. For each subset, a replication of $\hat{\theta}$ is obtained by computing $\hat{\theta}_b = \hat{\theta}(w_1[b], \dots, w_n[b])$. As a result, a set $\mathcal{R} = \{\hat{\theta}_1, \dots, \hat{\theta}_B\}$ of $B$ replications of $\hat{\theta}$ is obtained.

Denote by $\Phi$ the standard normal cdf, and by $\Phi^{-1}$ its inverse:

$$ \Phi(z) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{z} \exp\left(-\frac{t^2}{2}\right) dt. $$

Using the set $\mathcal{R}$ and the point estimate $\hat{\theta} = \hat{\theta}_n(w)$, one can estimate a “bias correction constant” $z_0$:

$$ \hat{z}_0 = \Phi^{-1} \left( \frac{\#\{\hat{\theta}_b \in \mathcal{R} \text{ s.t. } \hat{\theta}_b \le \hat{\theta}\}}{B} \right). $$

Then, for $\beta \in ]0; 1[$, the “corrected quantile estimate” $\hat{q}(\beta)$ is defined as:

$$ \hat{q}(\beta) = \Phi(2\hat{z}_0 + z_{\beta}), $$

where $z_\beta$ satisfies $\Phi(z_\beta) = \beta$. Finally, the central BC bootstrap confidence interval of level $1 - \alpha$ is estimated by the interval whose endpoints are the $\hat{q}(\alpha/2)$ and $\hat{q}(1 - \alpha/2)$ quantiles of $\mathcal{R}$.

These bootstrap confidence intervals are used to assess the precision of the Sobol' indices estimators. Let $k \in \mathcal{D}$. For the first-order Sobol' index estimate $\hat{S}_k^{(\gamma_1, \dots, \gamma_\kappa)}$, defined in Eq. (9), the bootstrap procedure is applied to $w = (y, y', y'_{\pi_k^{(\gamma_1)}}, \dots, y'_{\pi_k^{(\gamma_\kappa)}})$. For the second-order Sobol' index estimate $\hat{S}_{j_1, j_2}^{(\gamma_1, \dots, \gamma_\kappa)}$, defined in Eq. (10) with $\{j_1, j_2\} \in \mathcal{D}^2$, the bootstrap procedure is applied to $w = (y, y', y'_{\pi_{j_1}^{(\gamma_1)}}, \dots, y'_{\pi_{j_1}^{(\gamma_\kappa)}}, y'_{\pi_{j_2}^{(\gamma_1)}}, \dots, y'_{\pi_{j_2}^{(\gamma_\kappa)}})$.
# 4. Simulation study and application to benchmarks

The procedure REP18A is applied on two academic examples and compared with procedure TIS15. In parallel, a comparison between bootstrap and asymptotic confidence intervals is provided for the first academic example. Then, procedure REP18B is compared with procedure SAL02 on two academic examples, one of which is a test case for high-dimensional problems. Finally, REP18B is illustrated on an engineering application.

For clarity purposes, Table 1 lists the diverse procedures that will be applied or referred to in the following sections, along with the Sobol' indices they estimate, their computational cost $\nu$, and the nature of the sampling applied. Notations LHS, OA, MC, QMC stand for: Latin Hypercube Sampling, Orthogonal Array, Monte Carlo, quasi-Monte Carlo.

## 4.1. REP18A vs TIS15

This section compares REP18A with TIS15 on two academic examples: the Ishigami function introduced in [26] and the function introduced by Bratley et al. [27]. To reach a fair comparison, a common computational cost $\nu$ is used for REP18A and TIS15. As the computational cost for TIS15 equals $2n + 2q^2$, a balance between $n$ and $q^2$ is chosen.

The precision of REP18A is gauged by computing bootstrap confidence intervals according to Section 3.4. The number of permutations $\kappa$ (Eqs. (9) and

<table><caption>Table 1: Procedures for the estimation of Sobol' indices. The first three columns indicate the category of Sobol' indices estimated by each procedure. LHS, OA, MC, QMC stand for: Latin Hypercube Sampling, Orthogonal Array, Monte Carlo, quasi-Monte Carlo. The "+" symbol indicates that both designs are needed.</caption><thead><tr><th>name</th><th>1st</th><th>2nd</th><th>total</th><th>cost (ν)</th><th>sampling</th></tr></thead><tbody><tr><td>REP18A</td><td>✓</td><td>✓</td><td></td><td>2q<sup>2</sup></td><td>OA</td></tr><tr><td>REP18B</td><td>✓</td><td>✓</td><td>✓</td><td>q<sup>2</sup>(d + 2)</td><td>OA + MC/QMC</td></tr><tr><td>SAL02</td><td>✓</td><td></td><td>✓</td><td>n(d + 2)</td><td>MC/QMC</td></tr><tr><td>TIS15</td><td>✓</td><td>✓</td><td></td><td>2n + 2q<sup>2</sup></td><td>LHS + OA</td></tr></tbody></table>
(10) in Section 3.2) and of bootstrap replications B (see Section 3.4) are both fixed to 100. Besides, 100 replicate estimates of the Sobol' indices are calculated and the mean values are considered in the studies. The results are undertaken for the following computational costs $\nu \in \{162, 1058, 1922\}$.

### 4.1.1. Ishigami function

The Ishigami function is defined as follows:

$$f: \quad [-\pi, \pi]^3 \quad \rightarrow \quad \mathbb{R} \\ \boldsymbol{x} = (x_1, x_2, x_3) \quad \mapsto \quad \sin(x_1) + 7 \sin^2(x_2) + 0.1 x_3^4 \sin(x_1)$$

The main peculiarity of this test function is its strong nonlinearity and non-monotonicity. The input vector $\boldsymbol{x}$ is uniformly distributed over $[-\pi, \pi]^3$.

Fig. 1 provides means and bootstrap confidence intervals of first- and second-order Sobol' indices estimated with REP18A. Means are represented as black dots with their corresponding bootstrap confidence intervals (vertical bars). To highlight the efficiency of REP18A, Fig. 1 also displays in gray color the means and bootstrap confidence intervals obtained with TIS15.

REP18A yields good results for the first-order indices. From $\nu = 1058$ onward, the precision of the bootstrap confidence intervals allows one to distinguish the main effect of each input without any confusion. For $\nu = 1922$, the radii of the bootstrap confidence intervals are lower than $3 \times 10^{-2}$. Fig. 1 - (a),(c),(e) display the drastic gain obtained by using REP18A instead of TIS15, in particular for $S_3$.

The results for the second-order Sobol' indices are slightly less accurate. This was expected: in procedure REP18A, the estimation of the first-order Sobol' indices is obtained as an average over $\kappa$ different permutations (see Eq. (9)), but this is not the case for the estimation of the second-order Sobol' indices. As such and given the formula (10), the accuracy of the second-order Sobol' estimate $\hat{S}_{k,\ell}^{(\gamma_1, \dots, \gamma_\kappa)}$ is mainly driven by the precision of the corresponding closed second-order Sobol' estimate $\hat{S}_{\{k,\ell\}}$.

Nonetheless, these results are still satisfactory considering the low values selected for the computational cost $\nu$. Fig. 1 - (b),(d),(f) highlight once again
<figure><img src="image_3.png" /><figcaption>Figure 1: Ishigami function - Means (dots) and bootstrap confidence intervals (vertical bars) of first- (a, c and e) and second-order (b, d and f) Sobol' estimates. Grey refers to procedure TIS15, black to procedure REP18A. The black crosses mark the true values of the indices.</figcaption></figure>
<table><caption>Table 2: Ishigami function - Asymptotic and bootstrap confidence intervals (IC) of first-order Sobol' indices. True values are: $S_1 = 0.3139$, $S_2 = 0.4424$, $S_3 = 0$.</caption><thead><tr><th>ν</th><th>IC</th><th>S<sub>1</sub></th><th>S<sub>2</sub></th><th>S<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="2">162</td><td>bootstrap</td><td>[0.25, 0.39]</td><td>[0.35, 0.52]</td><td>[-0.02, 0.03]</td></tr><tr><td>asymptotic</td><td>[0.14, 0.49]</td><td>[0.29, 0.58]</td><td>[-0.28, 0.31]</td></tr><tr><td rowspan="2">1058</td><td>bootstrap</td><td>[0.28, 0.34]</td><td>[0.41, 0.48]</td><td>[-0.01, 0.01]</td></tr><tr><td>asymptotic</td><td>[0.25, 0.39]</td><td>[0.38, 0.49]</td><td>[-0.12, 0.11]</td></tr><tr><td rowspan="2">1922</td><td>bootstrap</td><td>[0.29, 0.34]</td><td>[0.42, 0.47]</td><td>[-0.01, 0.00]</td></tr><tr><td>asymptotic</td><td>[0.26, 0.37]</td><td>[0.40, 0.49]</td><td>[-0.09, 0.09]</td></tr></tbody></table>

the efficiency of REP18A compared to TIS15. The confidence intervals obtained with TIS15 are twice larger than those obtained with REP18A.

Table 2 provides a comparison between bootstrap and asymptotic confidence intervals for the estimation of first-order Sobol' indices with REP18A. This table shows a gain in precision when using bootstrap ones. Moreover, unlike bootstrap confidence intervals, asymptotic ones need additional evaluations of the model. The asymptotic theory for the estimation of second-order Sobol' indices via the replication method also requires additional regularity assumptions on the model (see, for example, [12]). We refer to the appendix for a detailed description on the way the asymptotic confidence intervals are constructed.

### 4.1.2. Bratley et al. function

The Bratley et al. function is defined by:

$$f: \quad \mathbb{H}_6 \to \mathbb{R} \\ \boldsymbol{x} = (x_1, \dots, x_6) \mapsto \sum_{i=1}^{6} (-1)^i \prod_{j=1}^{i} x_j$$

$x$ is assumed to be uniformly distributed over $\mathbb{H}_6$. The importance of each input $x_j$ depends on its own rank. More explicitly, $x_1$ is more influential than $x_2$, which is more influential than $x_3$, and so on. Fig. 2 displays the means and bootstrap confidence intervals of the first- and second-order Sobol' indices estimated with REP18A (black color) and TIS15 (grey color).
<figure><img src="image_4.png" /><figcaption>Figure 2: Bratley et al. function - Means (black dots) and bootstrap confidence intervals (vertical bars) of first- and second-order Sobol' estimates obtained with REP18A. Grey refers to procedure TIS15, black to procedure REP18A. The black crosses mark the true values of the indices.</figcaption></figure>
The results of REP18A for the first-order indices are similar to those obtained for the Ishigami function. For $\nu = 1058$, the radii of the confidence intervals are lower than $3 \times 10^{-2}$ and the non-influential inputs are well identified. The results of REP18A for the second-order indices are slightly less accurate, with radii ranging from 0.025 to 0.06 for $\nu = 1058$. The origin of this discrepancy is the same as the one evoked in Section 4.1.1.

For that example, most of second-order Sobol' indices are close to zero, which make their estimation challenging. Nonetheless, it can still be observed that REP18A provides a far better precision than TIS15.

## 4.2. REP18B vs SAL02

This section compares REP18B with SAL02 on both the Ishigami function and a high-dimensional test case. For REP18B, the number of permutations $\kappa$ and of bootstrap replications $B$ are once again both fixed to 100. For both procedures, 100 replicate estimates of the first-order and total-effect Sobol' indices are calculated and the mean values are considered.

### 4.2.1. Ishigami function

The improvement brought by REP18B over SAL02 is illustrated on the Ishigami function. The results are computed for the following costs: $\nu \in \{405, 845, 1445, 2645, 3645, 4805\}$. The precision of each procedure is assessed by drawing boxplots of the estimation errors of first-order and total-effect indices. These errors correspond to the absolute difference between true values of Sobol' indices and their estimates:

$$ \varepsilon_{S_u} = |S_u - \hat{S}_u|, \quad u \subset \mathcal{D}, \qquad (11) $$

where $S_u$ is the true value and $\hat{S}_u$ the estimate.

Fig. 3 shows the results obtained with both REP18B and SAL02. It is clear that REP18B performs the best for the first-order indices, while matching SAL02 precision for the total-effect indices.

It is worth mentioning that SAL02 has been extended to estimate the second-order indices for a total computational cost $\nu = n(2d + 2)$ (see Theorem 2 in
<figure><img src="image_5.png" /><figcaption>Figure 3: Ishigami function - Boxplots of estimation errors for first- and total-effect indices. White boxplots refer to REP18B, grey boxplots refer to variant B of SAL02. The x-axis indicates the size *n* of the designs.</figcaption></figure>
[1]), which is far more costly than REP18B for moderate dimension $d$. This extension outperforms REP18B only for high dimensional models from around 100 input parameters due to the construction constraint $q \ge d - 1$.

### 4.2.2. High-dimensional example

This section treats the case of a high-dimensional model by considering a modified version of the Morokoff & Caflisch function [28]:

$$f: \quad \mathbb{H}_d \quad \rightarrow \quad \mathbb{R} \\ \boldsymbol{x} = (x_1, \dots, x_d) \quad \mapsto \quad \prod_{j=1}^{d} x_j^{\frac{1}{\alpha_j}}$$

$\boldsymbol{x}$ is assumed to be uniformly distributed over $\mathbb{H}_d$. The importance of each input $x_j$ depends on its assigned coefficient. The higher the coefficient $\alpha_j$, the less influential the input $x_j$. A total of 50 inputs is considered. Values assigned to coefficients $\alpha_j$ are drawn from the set $\{2, 5, 6, 8, 10\}$ with the following approximation: four (*resp.* six, eight, twelve, twenty) among the fifty $\alpha_j$s are equal to 2 (*resp.* 5, 6, 8, 10). In that example, the sum of the first-order indices only explains 56% of the model variance and the second-order interactions have negligible effects. This makes this example interesting for procedures REP18B and SAL02.

The results are computed for $\nu = 53^2 \times (50 + 2) = 146 068$. The size of the designs, $n = 53^2$, are chosen to satisfy the construction constraint $q \ge d - 1$ in REP18B. Fig. 4 displays the means and bootstrap confidence intervals of the first-order Sobol' indices estimated with each procedure. The dashed lines represent the true value of the Sobol' indices with respect to the coefficients $\alpha_j$s.

The four influential inputs, corresponding to $\alpha_j = 2$, are identified by both procedures. For the remaining inputs, REP18B allows one to distinguish the sets of indices according to their coefficient $\alpha_j$. At the opposite, the estimates obtained with SAL02 are not distinguishable from one another. Overall, a net gain in precision is observed when using REP18B.

Fig. 5 displays the means and bootstrap confidence intervals of the total-effect Sobol' indices estimated with each procedure. As for the Ishigami func-
<figure><img src="image_6.png" /><figcaption>Figure 4: Morokoff & Caflish *et al.* function - Means and bootstrap confidence intervals of first-order Sobol' estimates obtained with REP18A (black color) and SAL02 (grey color) for $\nu = 146\,068$. The dashed lines represent the true values of the indices.</figcaption></figure>

tion, both procedures give similar results. Note that SAL02 could have been applied with a lower cost $\nu$ but at the expense of a less accurate estimation.

### 4.3. Application: Level E model

This section illustrates procedure REP18B on an engineering example. The Level E model has been used as a benchmark for sensitivity analysis by several authors (see [2] for a review). The model predicts the radiological dose to humans over geological time scales due to the underground migration of four radionuclides from a nuclear waste disposal site through two geosphere layers characterized by different hydro-geological properties.

The Level E model is mathematically represented by a set of partial differential equations, modeling the different chemical processes that govern the underground migration of the nuclides (a detailed description of the equations can be found in [2][Section 3.4]). The quantity of interest $Y(t)$ corresponds to
<figure><img src="image_7.png" /><figcaption>Figure 5: Morokoff & Caflish et al. function - Means and bootstrap confidence intervals of total-effect Sobol' estimates obtained with REP18A (black color) and SAL02 (grey color) for $v = 146\,068$. The dashed lines represent the true values of the indices.</figcaption></figure>

the annual radiological dose due to the four radionuclides. $Y(t)$ is evaluated at $26$ time frames ranging from $2 \times 10^4 = 20K$ to $9 \times 10^6 = 9M$ years. Following the simplification proposed in [2], we consider the twelve independent parameters listed in Table 3.

Procedure REP18B is applied to estimate first-order, second-order and total-effect Sobol' indices of the twelve input parameters. Two replicated designs of size $n = 5041$ are constructed first to estimate first-order and second-order Sobol' indices. Then, twelve additional pick-freeze samples are generated to estimate the total Sobol' indices. The corresponding computation cost reads $v = 70\,574$. Since the quantity of interest $(Y(t_i), i = 1, \dots, 26)$ is a vector, results are displayed as cumulative area plots in Figs. 6, 7 and 8.

For the first-order indices, the results obtained with REP18B match those presented in [2][Fig. 5.3 page 143]. The observations are the following: only
<table><caption>Table 3: Level E model - Parameters and associated distributions.</caption><thead><tr><th>parameters</th><th>distribution</th><th>range of variation</th><th>description</th></tr></thead><tbody><tr><td>T</td><td>uniform</td><td>[100, 1000]</td><td>containment time (yr)</td></tr><tr><td>k<sub>I</sub></td><td>log-uniform</td><td>[10<sup>-3</sup>, 10<sup>-2</sup>]</td><td>leach rate for iodine (mols/yr)</td></tr><tr><td>k<sub>C</sub></td><td>log-uniform</td><td>[10<sup>-6</sup>, 10<sup>-5</sup>]</td><td>leach rate for Np chain (mols/yr)</td></tr><tr><td>v<sup>(1)</sup></td><td>log-uniform</td><td>[10<sup>-3</sup>, 10<sup>-1</sup>]</td><td>water velocity in the first geosphere layer (m/yr)</td></tr><tr><td>l<sup>(1)</sup></td><td>uniform</td><td>[100, 500]</td><td>length of the first geosphere layer (m)</td></tr><tr><td>R<sub>I</sub><sup>(1)</sup></td><td>uniform</td><td>[1, 5]</td><td>retention factor for iodine in the first layer</td></tr><tr><td>R<sub>C</sub><sup>(1)</sup></td><td>uniform</td><td>[3, 30]</td><td>retention factor for the chain elements in the first layer</td></tr><tr><td>v<sup>(2)</sup></td><td>log-uniform</td><td>[10<sup>-2</sup>, 10<sup>-1</sup>]</td><td>water velocity in the second geosphere layer (m/yr)</td></tr><tr><td>l<sup>(2)</sup></td><td>uniform</td><td>[50, 200]</td><td>length of the second geosphere layer (m)</td></tr><tr><td>R<sub>I</sub><sup>(2)</sup></td><td>uniform</td><td>[1, 5]</td><td>retention factor for iodine in the first layer</td></tr><tr><td>R<sub>C</sub><sup>(2)</sup></td><td>uniform</td><td>[3, 30]</td><td>retention factor for the chain elements in the second layer</td></tr><tr><td>W</td><td>log-uniform</td><td>[10<sup>5</sup>, 10<sup>7</sup>]</td><td>stream flow rate (m<sup>3</sup>/yr)</td></tr></tbody></table>

parameters $v^{(1)}$ and $W$ are influential for the first nine time frames, a drop is observed at $t = 200K$ and the parameters $l^{(1)}$, $v^{(2)}$, $R_C^{(1)}$ become influential only after the drop. The main effects of the remaining parameters are negligible. The sum of the first-order indices being always lower than 0.25, it is relevant to study the second-order indices.

Fig. 7 shows the 19 most influent second-order interactions. One can observe that the model is driven by second-order interactions from $t = 300K$ to $t = 2M$, where approximately 75% of the model variance is explained. Since the sum of first- and second-order indices is lower than 0.7 on the two time frames $[20K, 200K]$ and $[4M, 9M]$, it is relevant to estimate the total-effect indices.

Fig. 8 displays the results for the total-effect Sobol' indices. Once again, the results obtained with REP18B match those presented in [2][Fig. 5.4 page 143]. Correspondingly to the first-order estimates, a strong increase of the total-effect estimates is observed at $t = 200K$. The parameters $R_C^{(2)}$, $l^{(2)}$, $R_I^{(1)}$, $R_I^{(2)}$ have significant interactions effects of order higher than two. On the other hand, the parameters $k_C$, $T$, $k_I$ are not influential at all.
Figure 6: Level E model - Cumulative area plot of first-order estimates obtained with REP18B. The y-axis indicates the sum of the estimates. The x-axis refers to the time frame in log scale.

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

Figure 7: Level E model - Cumulative area plot of the 19 most influent second-order estimates obtained with REP18B. The y-axis indicates the sum of the estimates. The x-axis refers to the time frame in log scale.

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Level E model - Cumulative area plot of total-effect estimates obtained with REP18B. The y-axis indicates the sum of the estimates. The x-axis refers to the time frame in log scale.</figcaption></figure>

5. Conclusion

When estimating Sobol' indices, the question, how many model evaluations must be performed to reach a decent precision is often raised by practitioners. This question is all the more critical as the number of evaluations available is often limited (whether by time or budget constraints). The extension of the replication procedure proposed in this article offers a practical solution to the estimation of the full set of first- and second-order Sobol' indices. Aside from halving the cost of the original replication procedure [12], our approach was shown to drastically enhance the estimates precision. The assessment of this precision was made by computing bootstrap confidence intervals. Compared to asymptotic intervals, bootstrap confidence intervals were found more reliable to gauge the quality of the estimation.

Additionally, we also discuss the extension of the procedure to the estimation
of the total-effect Sobol' indices. This was achieved by applying Saltelli's scheme
420

[2] to one of the orthogonal arrays of strength two already at hand. The extended approach allows one to evaluate the overall set of first-order, second-order and total-effect Sobol' indices.

425

References

[1] A. Saltelli, Making best use of model evaluations to compute sensitivity indices, Computer Physics Communications 145 (2002) 280-297.

[2] A. Saltelli, K. Chan, E. M. Scott, Sensitivity Analysis in Practice, Wiley, 2004.

[3] I. M. Sobol', Sensitivity analysis for nonlinear mathematical models, Mathematical Modeling and Computational Experiment 1 (1993) 407-414.

[4] A. B. Owen, Variance components and generalized sobol'indices, SIAM/ASA Journal on Uncertainty Quantification 1 (1) (2013) 19-41.

430

[5] I. Sobol', Y. Levitan, On the use of variance reducing multipliers in monte carlo computations of a global sensitivity index, Comput. Phys. Commun. 117 (1) (1999) 52-61.

[6] G. Glen, K. Isaacs, Estimating sobol sensitivity indices using correlations, Environ. Model. Softw. 37 (2012) 157-166.

435

[7] A. Janon, T. Klein, A. Lagnoux, M. Nodet, C. Prieur, Asymptotic normality and efficiency of two sobol index estimators, ESAIM: Probability and Statistics 18 (2014) 342-364.

[8] M. D. McKay, Evaluating prediction uncertainty, Technical Report NUREG/CR-6311, US Nuclear Regulatory Commission and Los Alamos National Laboratory (1995) 1-79.

440

[9] M. D. Morris, L. M. Moore, M. D. McKay, Sampling plans based on balanced incomplete block designs for evaluating the importance of computer
model inputs, Journal of Statistical Planning and Inference 136 (9) (2006) 3203-3220.

445

[10] M. D. Morris, L. M. Moore, M. D. McKay, Using orthogonal arrays in the sensitivity analysis of computer models, Technometrics 50 (2) (2008) 205-215.

450

[11] T. A. Mara, O. Rakoto Joseph, Comparison of some efficient methods to evaluate the main effect of computer model factors, Journal of Statistical Computation and Simulation 78(2) (2008) 167-178.

[12] J.-Y. Tissot, C. Prieur, A randomized orthogonal array-based procedure for the estimation of first-and second-order sobol'indices, Journal of Statistical Computation and Simulation 85 (7) (2015) 1358-1381.

455

[13] G. Archer, A. Saltelli, I. Sobol, Sensitivity measures, ANOVA-like techniques and the use of bootstrap, Journal of Statistical Computation and Simulation 58 (2) (1997) 99-120.

[14] W. F. Hoeffding, A class of statistics with asymptotically normal distributions, Annals of Mathematical Statistics 19 (1948) 293-325.

460

[15] A. Owen, Lattice sampling revisited: Monte carlo variance of means over randomized orthogonal arrays, The Annals of Statistics (1994) 930-945.

[16] H. Monod, C. Naud, D. Makowski, Uncertainty and sensitivity analysis for crop models, in: D. Wallach, D. Makowski, J. W. Jones (Eds.), Working with Dynamic Crop Models: Evaluation, Analysis, Parameterization, and Applications, Elsevier, Amsterdam, 2006, Ch. 4, pp. 55-99.

465

[17] C. Prieur, S. Tarantola, Variance-based sensitivity analysis: Theory and estimation algorithms, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Handbook of Uncertainty Quantification, Springer, 2016, pp. 1-23.

[18] L. Gilquin, L. A. J. Rugama, E. Arnaud, F. J. Hickernell, H. Monod, C. Prieur, Iterative construction of replicated designs based on sobol'sequences, Comptes Rendus Mathematique 355 (1) (2017) 10-14.

470
[19] K. Kishen, On latin and hyper-graeco-latin cubes and hyper cubes, Current Science 11 (3) (1942) 98-99.

[20] C. R. Rao, Hypercubes of strength d leading to confounded designs in factorial experiments, Bull. Calcutta Math. Soc 38 (3) (1946) 67-78.

[21] R. C. Bose, K. A. Bush, Orthogonal arrays of strength two and three, The Annals of Mathematical Statistics (1952) 508-524.

[22] L. Gilquin, C. Prieur, E. Arnaud, Replication procedure for grouped sobol'indices estimation in dependent uncertainty spaces, Information and Inference (2015) iav010.

[23] A. S. Hedayat, N. J. A. Sloane, J. Stufken, Orthogonal arrays, Springer-Verlag, New York, 1999, theory and applications, With a foreword by C. R. Rao.

[24] B. Efron, Nonparametric standard errors and confidence intervals, Canadian Journal of Statistics 9 (2) (1981) 139-158.

[25] B. Efron, R. Tibshirani, Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy, Statistical science 1 (1) (1986) 54-75.

[26] T. Ishigami, T. Homma, An importance quantification technique in uncertainty analysis for computer models, in: Uncertainty Modeling and Analysis, 1990. Proceedings., First International Symposium on, 1990, pp. 398-403.

[27] P. Bratley, B. L. Fox, H. Niederreiter, Implementation and tests of low-discrepancy sequences, ACM Trans. Model. Comput. Simul. 2 (3) (1992) 195-213.

[28] W. J. Morokoff, R. E. Caflisch, Quasi-monte carlo integration, J. Comput. Phys. 122 (2) (1995) 218-230.
[29] C. J. Stone, Additive regression and other nonparametric models, The Annals of Statistics 13 (2) (1985) 689-705.

[30] T. Hastie, gam: Generalized Additive Models, r package version 1.14-4 (2017).

URL https://CRAN.R-project.org/package=gam

## 6. Appendix

This appendix describes the construction of asymptotic confidence intervals for first-order Sobol' indices, estimated with the replication procedure.

Let $\ell \in \{1, \dots, d\}$; $S_\ell$ is estimated with the replication procedure and formulas (5), (6) and (7) of Section 2.2. Let $\mathbf{x} = (x_\ell, \mathbf{x}_{-\ell}) = (X, Z) \in \mathbb{R} \times \mathbb{R}^{d-1}$. Denote by $X', Z', Z'', Z'''$ independent copies of $X$ (resp. $Z$). Asymptotic confidence intervals are constructed as follows:

let $g_T : \mathbb{R}^{2d-1} \to \mathbb{R}^3$ defined by

$$g_T(x, z, z') = \begin{pmatrix} (f(x, z) - \mu)(f(x, z') - \mu) \\ f(x, z) + f(x, z') - 2\mu \\ (f(x, z) - \mu)^2 + (f(x, z') - \mu)^2 \end{pmatrix},$$

where $x \in \mathbb{R}$, $z, z' \in \mathbb{R}^{d-1}$ and $\mu = \mathbb{E}(f(\mathbf{x}))$.

Denote by $g_{Tadd}$ the best additive (in $X, Z, Z'$) approximation as defined in [29], and set $g_{Trem} = g_T - g_{Tadd}$. Additionally, we define:

$$\Phi_T(x, y, z) = \frac{x - \left(\frac{y}{2}\right)^2}{\frac{z}{2} - \left(\frac{y}{2}\right)^2},$$

and:

$$\mathbf{a}_T = (\mathrm{Cov}(Y, Y'), 0, 2\mathrm{Var}(Y))^T,$$

where Cov and Var denote the covariance and variance operators, respectively, $Y = f(X, Z)$ and $Y' = f(X, Z')$. Then, the asymptotic variance in the central limit theorem for $\hat{S}_\ell$ (see Eq. (7)) is defined as

$$\sigma_T^2 = (\nabla \Phi_T(\mathbf{a}_T))^T \Sigma_T \nabla \Phi_T(\mathbf{a}_T),$$
where $\Sigma_T$ is the $(3 \times 3)$ covariance matrix of $g_{Trem}$. A proof of these expressions is given in (ii) of Proposition 3.2. of [12].

The last remaining step is to explain how the matrix $\Sigma_T$ can be estimated. To do so, $g_T$ is decomposed as:

$$g_T(X, Z, Z') = m + g_{Tadd}(X, Z, Z') + g_{Trem}(X, Z, Z'),$$

where $m = \mathbb{E}(g_T(X, Z, Z'))$ and

$$g_{Tadd}(X, Z, Z') = g_{Ta1}(X) + g_{Ta2}(Z) + g_{Ta3}(Z')$$

is the best additive approximation of $g_T$.

We want to estimate the covariance matrix of the $\mathbb{R}^m$-valued random variable $g_{Trem}(X, Z, Z')$. For $i, j \in \{1, \dots, m\}$, we denote by $C_{i,j}$ the $(i, j)$-coefficient of this matrix.

Let $(X, Z, Z', Z'', Z''')_k, k = 1, \dots, n$ be a $n$-sample of $(X, Z, Z', Z'', Z''')$. Let us define $Y_i, Y_j, Y'_j$ and $Y''_j$ as:

$$Y_i = g_T((X, Z, Z')_i) = m_i + g_{Ta1}(X)_i + g_{Ta2}(Z)_i + g_{Ta3}(Z')_i + g_{Trem}(X, Z, Z')_i,$$

$$Y_j = g_T((X, Z, Z')_j) = m_j + g_{Ta1}(X)_j + g_{Ta2}(Z)_j + g_{Ta3}(Z')_j + g_{Trem}(X, Z, Z')_j,$$

$$Y'_j = g_T((X, Z'', Z''')_j) = m_j + g_{Ta1}(X)_j + g_{Ta2}(Z'')_j + g_{Ta3}(Z''')_j + g_{Trem}(X, Z'', Z''')_j,$$

$$Y''_j = g_T((X', Z, Z''')_j) = m_j + g_{Ta1}(X')_j + g_{Ta2}(Z)_j + g_{Ta3}(Z'')_j + g_{Trem}(X', Z, Z'')_j,$$

$$Y'''_j = g_T((X', Z'', Z')_j) = m_j + g_{Ta1}(X')_j + g_{Ta2}(Z'')_j + g_{Ta3}(Z')_j + g_{Trem}(X', Z'', Z')_j.$$

We also define:

$$A_{i,j} = \mathrm{Cov}(g_{Ta1}(X)_i, g_{Ta1}(X)_j),$$

$$B_{i,j} = \mathrm{Cov}(g_{Ta2}(Z)_i, g_{Ta2}(Z)_j),$$

$$D_{i,j} = \mathrm{Cov}(g_{Ta3}(Z')_i, g_{Ta3}(Z')_j).$$

Thanks to independence, and $L^2$-orthogonality between $g_{Trem}(X, Z, Z')$ and functions of $X$ (resp. $Z, Z'$) alone, we have:

$$\mathrm{Cov}(Y_i, Y_j) = A_{i,j} + B_{i,j} + D_{i,j} + C_{i,j},$$
$$ \mathrm{Cov}(Y_i, Y_j') = A_{i,j}, $$

$$ \mathrm{Cov}(Y_i, Y_j'') = B_{i,j}, $$

$$ \mathrm{Cov}(Y_i, Y_j''') = D_{i,j}. $$

Hence:

$$ C_{i,j} = \mathrm{Cov}(Y_i, Y_j) - \mathrm{Cov}(Y_i, Y_j') - \mathrm{Cov}(Y_i, Y_j'') - \mathrm{Cov}(Y_i, Y_j''') $$

which give rise to a natural empirical estimator $\hat{C}_{i,j}$ of $C_{i,j}$ using a sample of $Y_i, Y_j, Y_j'$ and $Y_j''$.

In practice, to estimate the asymptotic variances without additional model evaluations, one can reuse the samples of input variables used during sensitivity indices estimation (both first-order and second-order). However, the numerical experiments carried out in Section 4.1.1 are made using new Monte-Carlo samples, as the results are expected to be quite similar.

It is also worth noting that $g_{add}$ could be estimated using, for instance, the semi-parametric methods in the R package GAM [30]. However, experiments show that this choice induces some significant bias in the estimation of $\Sigma$, leading to a bias in the estimation of the asymptotic variance of $\hat{T}$.