# Simple binning algorithm and SimDec visualization for comprehensive
sensitivity analysis of complex computational models

Mariia Kozlovaa∗, Antti Aholab, Pamphile Royc, Julian Scott Yeomansd.
aLUT Business School, LUT University, Lappeenranta, Finland
bLaboratory of Steel Structures, LUT University, Lappeenranta, Finland
cQuansight, Austria
dSchulich School of Business, York University, Toronto, Canada

∗Corresponding author: Mariia Kozlova<EMAIL_ADDRESS>
###### Abstract

Models of complex technological systems inherently contain interactions and
dependencies among their input variables that affect their joint influence on
the output. Such models are often computationally expensive and few
sensitivity analysis methods can effectively process such complexities.
Moreover, the sensitivity analysis field as a whole pays limited attention to
the nature of interaction effects, whose understanding can prove to be
critical for the design of safe and reliable systems. In this paper, we
introduce and extensively test a simple binning approach for computing
sensitivity indices and demonstrate how complementing it with the smart
visualization method, simulation decomposition (SimDec), can permit important
insights into the behavior of complex engineering models. The simple binning
approach computes first-, second-order effects, and a combined sensitivity
index, and is considerably more computationally efficient than Sobol’ indices.
The totality of the sensitivity analysis framework provides an efficient and
intuitive way to analyze the behavior of complex systems containing
interactions and dependencies.

##  1 Introduction

Engineering models of complex technological systems are characterized by a
high degree of complexity and specificity. While the researchers developing
them are experts in such specific mathematical approaches, as finite element
methods, structural reliability and integrity modeling, and fitness-for-
service (prEN 1993-1-14, 2023; Mahadevan and Ni, 2003; Shen et al., 2023;
Shittu et al., 2021), their professional expertise is frequently deficient in
the mathematically comprehensive field of sensitivity analysis (SA) (Saltelli
et al., 2019). In many structural elements, potential errors in the underlying
structural analysis models can lead to failures that cause catastrophic
consequences in service. To prevent such outcomes, it is essential to conduct
an effective SA on those models. Furthermore, SA can also be used in model
development to identify the most important parameters and their interactions
and guide the process. SA can be used to explicitly test the models and
exclude those parameters that do not significantly affect the outputs to build
simplified approaches for engineering. The development of simple and intuitive
SA approaches is paramount for incorporating SA into standard practice to
fully complement the processes of building and analyzing a computational model
(Saltelli et al., 2019; Iooss et al., 2022).

Consequently, this paper provides several contributions to SA for assessing
the behavior of complex engineering models. Firstly, we introduce a novel
intuitive procedure that efficiently computes global variance-based
sensitivity indices for a given dataset (Section 2), together with its
motivation presented (Section 1.1). Secondly, in contrast to recent
mathematical developments (Da Veiga, 2015; Lamboni and Kucherenko, 2021; Mara
et al., 2015; Borgonovo et al., 2023; Barr and Rabitz, 2023), the new approach
is shown to capture dependent inputs automatically while simultaneously
revealing several interesting properties (Section 4). Thirdly, the entirety of
modern sensitivity analysis literature is focused solely on the computation of
sensitivity indices (Shi et al., 2023; Barr and Rabitz, 2023; Vuillod et al.,
2023; Wang and Jia, 2023; Jung and Taflanidis, 2023; Shang et al., 2023;
Ballester-Ripoll and Leonelli, 2022; Thapa and Missoum, 2022; Xiong et al.,
2022). In contrast, we advocate for the use of visualizations to investigate
the shape of the effects in addition to their strength (Section 1.2), and
illustrate this with a structural reliability model that possesses nested
heterogeneous interaction effects that could not be uncovered without such
visual representation (Section 5). Overall, this paper introduces a simple but
powerful framework for the SA of engineering models, that computes sensitivity
indices to prioritize the input factors followed by a smart visualization for
communicating the nature of the discovered effects. There are open-source
codes in several programming languages that implement the entire framework to
accompany the article.

###  1.1 Sensitivity indices computation

Understanding the behavior of computational models forms the basis for design
and decision-making. The decision-maker obtains information on: what sources
of uncertainties are most crucial and might require protection against; which
actionable parameters make the most difference and thus represent the perfect
levers for managing the system; and, which model parameters create the most
noise and require clarification to obtain a clearer representation of the
system being modeled.

The sole purpose of the scientific field of SA is to identify which input
variables affect the output(s) the most in the computational models (Saltelli
et al., 2008). For educational purposes, the conceptual idea behind the
sensitivity indices is often explained using scatter plots (Caers, 2018;
Saltelli, 2023) as depicted in Figure 1, borrowed from one of the talks
devoted to SA (Saltelli, 2023).

Figure 1: Conceptual representation of the sensitivity indices, adapted from
(Saltelli, 2023).

Figure 1 displays the intuition behind the variance-based sensitivity indices:
namely, bin the $X$, compute the average of $Y$ in every bin, compute the
variance of these average values of $Y$, and divide the resulting conditional
variance to the total variance of $Y$ to obtain the sensitivity index value.
Quite unexpectedly, however, no common methods have previously adopted this
logic in computing sensitivity indices. Among the widely used methods for
computing sensitivity indices: some involve game theory (Shapley et al., 1953)
and some Fourier analysis (FAST) (Cukier et al., 1973); the classic Sobol’
indices (Sobol, 1993; Homma and Saltelli, 1996) require juggling multiple
simulation matrices; some methods employ polynomial chaos expansion (Sudret,
2008; Crestaux et al., 2009); there is a so-called moment-independent measure
which is based on quantifying the shift in the distribution of the output
caused by an input (Borgonovo, 2007); some are based on derivatives (Sobol and
Kucherenko, 2010; Kucherenko and Song, 2016); there is a variogram-based
method (VARS) (Razavi and Gupta, 2016); while the most recent approach
utilizes discrepancy measures (Puy, Roy, and Saltelli, 2022) to quantify the
difference between the cumulative probability distributions of simulated data
and its running averages. Modern research on sensitivity analysis is actively
developing and predominantly targets an increase in the computational
efficiency of these various methods (Iooss et al., 2022; Saltelli et al.,
2021).

Marzban and Lahmer (2016) introduced a method that follows the logic of Figure
1, though this approach has been largely overlooked by the modeling community
(their article has only nine citations in Scopus). The authors call it
conceptual implementation of the variance-based sensitivity analysis and
demonstrate how the approximation of first-order sensitivity indices can be
computed with this approach (see Section 2.1 for details). Our focus is drawn
to this method because:

  * •

it is intuitive and works precisely as sensitivity indices are conceptually
introduced (Figure 1);

  * •

it is computationally efficient;

  * •

it works with a given dataset. As Plischke (2012) notes, for such methods: (i)
no access to the model is required; (ii) they can work on measured data; and
as a result, (iii) the method is indifferent to distribution types and
sampling methods. Furthermore, it enables smoother passing of information
between the model and sensitivity analysis, if they are performed in different
software or by different people.

This paper extends the concept of Marzban and Lahmer (2016) to capture second-
order effects. This modified new approach, referred to as the Simple binning
approach, introduces additional benefits:

  * •

the choice of the number of bins is no longer the user’s responsibility, it is
automated based on experiments by Marzban and Lahmer (2016);

  * •

it works with dependent inputs (although developments that allow sensitivity
analysis of models with dependent input variables exist (Da Veiga, 2015;
Lamboni and Kucherenko, 2021; Mara et al., 2015; Borgonovo et al., 2023), they
all require certain additional transformations, whereas the simple binning
approach captures dependent inputs ”as is” and preserves the conservation
property);

  * •

it works with categorical/discrete variables;

  * •

it is resistant to different numbers of observations in bins due to the usage
of weighted variance;

  * •

as a result, it can work with empirical datasets, or

  * •

it can be used to quantify the sensitivity of the output to intermediate
outputs in simulation models.

###  1.2 Importance of visualization

The modern practice of sensitivity analysis systematically stops at computing
sensitivity indices (Shi et al., 2023; Barr and Rabitz, 2023; Vuillod et al.,
2023; Wang and Jia, 2023; Jung and Taflanidis, 2023; Shang et al., 2023;
Ballester-Ripoll and Leonelli, 2022; Thapa and Missoum, 2022; Xiong et al.,
2022) and does not take an additional step of investigating the shape of the
effects in a model.

Model interactions may arise in many different forms (Kozlova et al., in
press). In one model, they might be the result of multiplication in which the
effect of one input variable linearly increases in the higher values of
another input variable. In another model, the effect of the variable might be
strong in a certain region of another input variable, but negligible in the
others. The direction of the influence of one input factor on the output might
be flipped in different regions of another input variable. In all of these
cases, the sensitivity index will show a positive second-order effect of a
certain degree, but no more. Namely, no indication of the type of interaction
effect can be obtained through sensitivity indices alone.

However, explicit knowledge of the nature of interaction effects is crucial
for effective decision-making and engineering design and, thus, visualization
of the effects is a ”must-have” practice for computational model analysis.
Heat maps (Pleil et al., 2011; M. P. Owen et al., 2019) and response surfaces
(Myers et al., 2016) are inherently three-dimensional visualizations that can
portray the joint effect of a pair of input variables on the model output.
However, these visualization types cannot handle simultaneous variations of
multiple inputs and, thus, fail to depict higher-order interactions.

These visualization limitations are overcome by the simulation decomposition
(SimDec) approach (Kozlova and Yeomans, 2022b). SimDec partitions the output
probability or frequency distribution into sub-distributions comprised of
combinations of specific regions of influential input variables, thereby
exposing the nature of interaction effects. The SimDec approach has
demonstrated value in multiple applications (Kozlova and Yeomans, 2019;
Deviatkin et al., 2021; Liu et al., 2022). In this paper, it will be shown
that the concurrent employment of the simple binning approach to identify the
most influential input variables with SimDec to display the nature of those
effects results in a simple-to-implement, intuitive, and sophisticated
framework for conducting SA.

##  2 Simple binning approach to sensitivity indices

###  2.1 Existing conceptual implementation for first-order effects

Marzban and Lahmer (2016) introduced a conceptual implementation to Sobol’
first-order indices. Unlike the multiple matrices normally needed for
calculating Sobol’ indices, their approach requires only a single
$N_{observations}*K_{inputs}$ matrix generated from one Monte Carlo
simulation. The approach involves binning along the $X$’s range, computing
averages of $Y$s within every bin, and then taking the variance of those
averages (as shown in Figure 1). The first-order sensitivity index is then
computed as the conditional variance divided by the variance of $Y$:

$\ S_{x_{i}}=\frac{\operatorname{Var}(\mathbb{E}(Y\mid
X_{i}))}{\operatorname{Var}(Y)}$ (1)

Furthermore, Marzban and Lahmer (2016) define the optimal number of bins for
different $N_{observations}$ and $K_{inputs}$ as shown in Table 1, and
demonstrate this approach on several benchmark problems with known analytical
solutions for their sensitivity indices.

Table 1: Optimal number of bins for different combinations of the number of input factors and simulation runs according to Marzban and Lahmer (2016). |  Number of input variables, $K_{\text{inputs}}$
---|---
Number of observations, $N_{\text{obs}}$ |  3 |  6 |  12
1000 |  10 |  10 |  10
2500 |  25 |  10 |  10
5000 |  50 |  10 |  10
7500 |  50 |  25 |  10
10000 |  50 |  50 |  10
25000 |  100 |  50 |  25
50000 |  100 |  50 |  50

###  2.2 Extension to second-order effects

Following on from the first-order index calculation approach, we extend the
binning idea into the computation of second-order effects. This implies that
instead of binning a line defined by a single input, we must now bin a two-
dimensional area defined by two inputs. Figure 2.

Figure 2: Simple binning approach for first-order indices (left) and second-
order indices (right).

The calculation steps follow the same procedure for computing averages in the
bins and their variance, except we employ a weighted variance to account for
possibly unequal numbers of observations in bins relevant to categorical
variables with different frequencies for different categories, and for second-
order effects:

$\operatorname{Var}_{w}(Y_{X_{i}})=\frac{\sum_{b=1}^{N_{\text{bins}}}N_{\text{obs}_{b}}(\overline{Y_{b}}-\overline{Y})^{2}}{\sum_{b=1}^{N_{\text{bins}}}N_{\text{obs}_{b}}}.$
(2)

where the binning happens over $X_{i}$ and the averages are taken of $Y$
values as illustrated in Figure 2.

The second-order effect is defined in the classic way.

$\ S_{{X_{i}}{X_{j}}}=\frac{V(E(Y|X_{i}X_{j})}{V(Y)}-S_{X_{i}}-S_{X_{j}}$ (3)

or

$\
S_{{X_{i}}{X_{j}}}=\frac{V(E(Y|X_{i}X_{j})}{V(Y)}-\frac{V(E(Y|X_{i})}{V(Y)}-\frac{V(E(Y|X_{j})}{V(Y)}$
(4)

The number of bins is defined based on the optimal number of observations per
bin. For the first-order effects, the number of bins is a linear approximation
of the optimal experimental number of bins as a function of the number of
simulation runs and the number of inputs, Table 1, floored at 10. For the
calculation of second-order effects, to preserve approximately the same number
of observations in two-dimensional bins, the range of each input variable is
broken down into a number of intervals equal to a square root (rounded to the
nearest integer) of the previously computed number of bins for the first-order
effects. The same number of bins should be applied to all terms in Equation 4.
For example, a stylized scheme in Figure 2 depicts four bins for $X_{i}$ for
the first-order effect, which translates into two bins for $X_{i}$ and two for
$X_{j}$ for second-order effects, thus, the bin size remains the same.

An aggregate or combined sensitivity index can now be calculated which, for
each input variable, sums up its first-order effect and halves second-order
indices with all other input variables:

$\ S_{\text{combined}_{X_{i}}}=S_{X_{i}}+\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{K_{\text{inputs}}}\frac{S_{{X_{i}}{X_{j}}}}{2}$ (5)

This combined sensitivity index aggregates the individual and the interaction
effects of each input variable. The sum of combined indices for all model
input variables would equal $1$ if there are no effects higher than second-
order in the model and no unaccounted variation. However, due to the
approximation nature of the method, potentially noisy results due to binning
(Marzban and Lahmer, 2016), and the imperfections of sampling, the sum of
combined sensitivity indices can slightly deviate from $1$ even in the absence
of higher-order effects. Irrespectively, we advocate for the usage of this
approach in combination with SimDec to determine the selection of the prime
input variables for decomposition. As long as the ranking of inputs’
importance is preserved, noise does not prove to be an impediment.

##  3 Testing the simple binning approach

###  3.1 Capturing interactions – The toy model

To demonstrate the efficacy of the simple binning approach in capturing
second-order effects, we examine the simple toy model presented in the seminal
sensitivity analysis textbook by Saltelli et al. (2004). At the heart of the
toy model are summations and multiplications which are basic operations in any
computational model. The model equation and the input distributions are as
follows:

$\ Y=C_{s}P_{s}+C_{t}P_{t}+C_{j}P_{j}\\\ $ (6)

where $P_{s}\sim\ N(0;4)$, $P_{t}\sim\ N(0;2)$, $P_{j}\sim\ N(0;1)$,
$C_{s}\sim\ N(250;200)$, $C_{t}\sim\ N(400;300)$, and $C_{j}\sim\ N(500;400)$.

A side-by-side comparison of sensitivity index results obtained by Sobol’
approach as reported in (Saltelli et al., 2004) and the new simple binning
method are presented in Table 2.

Table 2:  Comparison of the Sobol’ indices and simple binning approach on a toy model. Effect | Sobol’ indices (Saltelli et al., 2004) | Simple binning method | Delta
---|---|---|---
First-order effects
$P_{s}$ | 36 % | 35 % |  $-1$ %
$C_{s}$ | 0 % | 1 % | 1 %
$P_{t}$ | 22 % | 20 % |  $-2$ %
$C_{t}$ | 0 % | 1 % | 1 %
$P_{j}$ | 8 % | 8 % | 0 %
$C_{j}$ | 0 % | 2 % | 2 %
Sum of first-order effects | 66 % | 67 % | 1 %
Second-order effects (of selected pairs of variables)
$P_{s}C_{s}$ | 18 % | 16 % |  $-2$ %
$P_{t}C_{t}$ | 11 % | 10 % |  $-1$ %
$P_{j}C_{j}$ | 5 % | 6 % | 1 %
Sum of second-order effects | 34 % | 32 % |  $-2$ %
Sum of all effects | 100 % | 99 % |  $-1$ %
Model evaluations | 21000 | 1000 |

Table 2 demonstrates that the simple binning method produces meaningful
sensitivity indices close to Sobol’ indices. However, the simple binning
method used only 1000 model evaluations, while 21000 data points were required
to compute the Sobol’ indices. The number of bins used was ten for first-order
effects and three for second-order effects.

The indices produced by both approaches in Table 2 are based on limited model
evaluations with simple random sampling. Because different simulation runs can
generate to some degree different datasets, we repeated the exercise 1000
times to ascertain how stable the index estimations were. Figure 3 provides a
comparison of the results obtained with 10K and 1K model evaluations with the
simple binning approach, respectively, contrasted by Sobol’ indices computed
on an efficient quasi-random (Sobol’ sequence) sampling of the size 1792 (14
different simulations required for 6 inputs variables, 128 points each). The
simple binning approach even based on a smaller 1K sample produces more
accurate estimates than Sobol’ indices. A larger sample narrows the deviations
in sensitivity indices estimation even further.

Figure 3: Sobol’ indices obtained with quasi-random sampling (QMC) of the size
1792 and sensitivity indices obtained with the simple binning approach and
simple random sampling (MC) for 1K and 10K model evaluations of the toy model.

In addition, we compared how the simple binning approach performs with (i)
simple random sampling (Olken and Rotem, 1986; Singh and Singh, 2003), which
is often used for Monte Carlo simulation, (ii) quasi-random sampling (Sobol’,
1967; A. B. Owen, 2023) with the Sobol’ low discrepancy sequence, which is
reportedly a more efficient sampling strategy that fills the space more
uniformly, and (iii) full factorial design, frequently used in engineering
contexts (Alidoosti et al., 2013; Suard et al., 2013; Tong, 2006), Figure 4.
For this experiment, we change the distribution of inputs from normal to
uniform, selecting as minimum and maximum values $\pm$ two standard deviations
of the corresponding normal distributions.

Figure 4: Toy model sensitivity indices estimation by the simple binning
approach with simple random sampling Monte Carlo simulation (MC), quasi-random
sampling or quasi-Monte Carlo (QMC), and full factorial design (FFD) sampling.

From Figure 4, one can observe that simple random sampling results in noisier
estimates and that this noise increases substantially the smaller the sample.
Full factorial design sampling results in deterministic estimates for
sensitivity indices, but the simple binning approach performed on such a
sample underestimates first-order indices and overestimates second-order
effects. Quasi-Monte Carlo (QMC), done with scrambling, gives more accurate
and clear estimates than other sampling strategies, even with a smaller
sample, and, thus, is recommended for sensitivity indices computation. QMC is
easy to implement and its coded functions are available in Python (SciPy,
2023), R (Chalabi et al., 2023), Julia (Robbe, 2018), and Matlab (MathWorks,
2013). However, when using any of these packages, it is important to have the
first point in the sequence sampled, otherwise, its performance drops (A. B.
Owen, 2020).

###  3.2 Capturing cyclic behavior – The Ishigami function

The Ishigami function, equation 7, is one of the benchmark functions often
used in sensitivity analysis studies for the validation of different methods,
since analytical solutions for its sensitivity indices can be readily
determined.

$\ Y=\sin(X_{1})+a\sin(X_{2})^{2}+bX_{3}^{4}\sin(X_{1})\\\ $ (7)

The function is periodic in nature, Figure 5, which represents a significant
challenge for approximate methods (Ziehn and Tomlin, 2009).

Figure 5: Ishigami function, equation 7, with $a=7$ and $b=0.1$ (the third
factor is fixed at 0 on each plot).

The simple binning algorithm is tested on the Ishigami function with
parameters $a=7$ and $b=0.1$, quasi-random sampling of size $10^{5}$,
$10^{4}$, and $10^{3}$, Figure 6.

Figure 6: Analytical and approximate sensitivity indices for Ishigami
function.

Figure 6 indicates that the simple binning algorithm with QMC sample size
$10^{4}$ or larger captures first-order effects very well, but underestimates
the highly-curved (see Figure 5 (center)) interaction effect between $X_{1}$
and $X_{3}$. A smaller sample translates into fewer bins for estimating the
first-order effects and the square root of that for second-order effects (a
sample size of 1000 is analyzed with 10 bins for first-order effects and only
3 bins for the second-order effects), so the high-frequency relationships are
downplayed. Thus, if a model is known to possess periodic or cyclic effects, a
larger sample size becomes essential in order to produce reliable estimates
for sensitivity analysis.

###  3.3 Capturing correlation – The mechanical engineering model

In this section, we reproduce and then modify the mechanical engineering
problem presented in Marzban and Lahmer (2016) by introducing dependent inputs
into their model and examining the resulting performance capabilities of our
sensitivity analysis algorithm in processing this additional complication.

####  3.3.1 Case background

Civil engineering structures play an essential role in most modern
infrastructures (e.g., in buildings, process industry, and shells). A wide
spectrum of steel frames are required within these structures. From an
engineering perspective, such frames must be designed to withstand numerous
different kinds of load actions due to such things as payloads, self-weight,
and environmental impacts. In statically loaded structures (characterized by
permanent load actions), it is necessary to analyze stresses and
displacements. Stress analyses are usually related to the ultimate limit state
(ULS), while analyses on displacement cover the serviceability limit states
(SLS) that are affiliated with the operational conditions and functionalities
(service) of structures. It is crucial to understand the effects of different
parameters on the mechanical system behavior. For some simple cases, the
effects can be determined analytically.

Frequently, however, nonlinear becomes impractical or impossible to solve
analytically in even the most straightforward systems. In such cases,
numerical modeling becomes a requisite tool for assessing structural behavior
under mechanical loads. From an engineering viewpoint, it is usually easy to
determine the most influential geometric and material parameters in the
underlying linear analyses (LAs) of mechanical behavior of the system, such as
deflections and stresses. However, geometrically nonlinear analyses (GNAs) may
uncover different behavior of the system with a different set of influential
parameters. Thus, the decision, which analysis to use is often critical for
engineering design and structural analysis.

To support the decision as to whether a nonlinear approach is required, a
simplified model that computes the ratio between the results of two different
analysis types (e.g., LA and GNA) can be employed. The sensitivity indices for
different parameters can then provide insight for decision-making.

A two-dimensional frame structure is employed to demonstrate such an approach
(Figure 7), as in (Marzban and Lahmer, 2016). The frame structure is comprised
of two vertical columns joined with a horizontal beam. The beam is loaded by a
uniformly distributed shear load along the beam and the top left corner of the
frame is loaded by a horizontal load. The structural behavior, focusing on the
displacements, of the frame structure is numerically solved via LA and GNA.

Figure 7: Shape and dimensions of the studied frame structure.

As the original CALFEM finite element (FE) code (Austrell et al., 2004)
applied in Marzban and Lahmer (2016) was not accessible, a similar frame
structure compiled by three beam elements was recreated (Figure 7). The
numerical FE models were parameterized in an API tool using the FEMAP 2022.2
(Siemens PLM) software. The system was studied employing nine different input
parameters. Namely, the height $H$ and width $W$ of the frame; the modulus of
elasticity $E$; the cross-sectional area of columns and beam, $A_{c}$ and
$A_{b}$, respectively; the second moment of areas of columns and beam, $I_{c}$
and $I_{b}$, respectively; lateral force, P; and uniformly distributed shear
load $q_{0}$. Each parameter was varied by $\pm 25\%$ of its mean value and
Table 3 displays all of the corresponding minimum and maximum values. A total
number of 10,000 models were simulated and analyzed using both LAs and GNAs.
As an output value, the ratio between lateral displacements in the top left
corner of the frame, obtained using GNA and LA, was used.

Table 3: Model parameter ranges in the frame model. | $H$ | $W$ | $E$ | $A_{c}$ | $I_{c}$ | $A_{b}$ | $I_{b}$ | $P$ | $q_{0}$
---|---|---|---|---|---|---|---|---|---
| (m) | (m) | (GPa) | (e-3 m2) | (e-5 m4) | (e-3 m2) | (e-5 m4) | (kN) | (kN/m)
Minimum | 2.55 | 3.0 | 150 | 1.5 | 1.2 | 4.5 | 4.05 | 7.5 | 37.5
Maximum | 4.25 | 5.0 | 250 | 2.5 | 2.0 | 7.5 | 6.75 | 12.5 | 62.5

####  3.3.2 Estimation of sensitivity indices

Table 4 displays the first- and second-order sensitivity indices computed for
the original model used in Marzban and Lahmer (2016). To clarify the visual
perception, all index values less than 2% have been grayed out. As can be seen
from Table 4, the frame height is the mode influencing parameter. The relative
magnitudes shown for the first-order effects coincide with the those obtained
by Marzban and Lahmer (2016). The negligible second-order indices signify an
absence of any interaction effects.

Table 4:  Sensitivity indices for the engineering model. | $H$ | $W$ | $E$ | $A_{c}$ | $I_{c}$ | $A_{b}$ | $I_{b}$ | $P$ | $q_{0}$
---|---|---|---|---|---|---|---|---|---
| (m) | (m) | (N/$m^{2}$) | (m2) | (m4) | (m2) | (m4) | (N) | (N/m)
First-order indices | 41 % | 16 % | 13 % | 0 % | 10 % | 0 % | 0 % | 0 % | 13 %
Second-order indices |  |  |  |  |  |  |  |  |
$H$ (m) |  | 1 % | 1 % | 0 % | 0 % | 0 % | 0 % | 0 % | 1 %
$W$ (m) |  |  | 0 % | 0 % | 0 % | 0 % | 0 % | 0 % | 0 %
$E$ (N/m2) |  |  |  | 0 % | 0 % | 0 % | 0 % | 0 % | 0 %
$A_{c}$ (m2) |  |  |  |  | 0 % | 0 % | 0 % | 0 % | 0 %
$I_{c}$ (m4) |  |  |  |  |  | 0 % | 0 % | 0 % | 0 %
$A_{b}$ (m2) |  |  |  |  |  |  | 0 % | 0 % | 0 %
$I_{b}$ (m4) |  |  |  |  |  |  |  | 0 % | 0 %
$P$ (N) |  |  |  |  |  |  |  |  | 0 %

To introduce dependency into the system, the $H$/$I_{c}$ ratio was fixed as
per the dimensions of the case structure (i.e., $I_{c}$ = 4.706e-6$\cdot$$H$).
With this change in effect, to investigate the second-order effects of input
parameters with correlated input parameters, a subsequent computational
experiment was conducted with a total of 10,000 simulations. The new
sensitivity indices calculated for this correlated system are presented in
Table 5.

Table 5:  Sensitivity indices for the modified engineering model with correlated inputs. | $H$ | $W$ | $E$ | $A_{c}$ | $I_{c}$ | $A_{b}$ | $I_{b}$ | $P$ | $q_{0}$
---|---|---|---|---|---|---|---|---|---
| (m) | (m) | (N/m2) | (m2) | (m4) | (m2) | (m4) | (N) | (N/m)
First-order indices | 21 % | 27 % | 23 % | 0 % | 21 % | 0 % | 1 % | 0 % | 22 %
Second-order indices |  |  |  |  |  |  |  |  |
$H$ (m) |  | 1 % | 0 % | 0 % | -20 % | 0 % | 0 % | 0 % | 2 %
$W$ (m) |  |  | 0 % | 0 % | 1 % | 0 % | 0 % | 0 % | 0 %
$E$ (N/m2) |  |  |  | 0 % | 0 % | 0 % | 0 % | 0 % | 0 %
$A_{c}$ (m2) |  |  |  |  | 0 % | 0 % | 0 % | 0 % | 0 %
$I_{c}$ (m4) |  |  |  |  |  | 0 % | 0 % | 0 % | 2 %
$A_{b}$ (m2) |  |  |  |  |  |  | 0 % | 0 % | 0 %
$I_{b}$ (m4) |  |  |  |  |  |  |  | 0 % | 0 %
$P$ (N) |  |  |  |  |  |  |  |  | 0 %

Comparing Tables 4 and 5, we can clearly observe the striking difference
arising in the second-order effects. While for the original model, all second-
order effects were negligible, in the correlated model, we can see a negative
effect of 20%. The negative second-order effect reflects the introduced
dependency between the inputs. It should be noted that this negative second-
order effect essentially zeroes out the first-order effect of the $I_{c}$ if
summing up all effects. Thus, the negative second-order effect represents a
correction of an overlapping first-order effect.

It can be further observed that the first-order indices have changed their
values as well. By effectively eliminating one input variable via correlation,
the overall output variance of the output has changed so that corresponding
portions of the explained variance have also changed accordingly.

This engineering case demonstrates that the binning approach for calculating
sensitivity indices is capable of capturing and identifying the impacts from
dependent inputs.

Even though such sensitivity analyses can be time-consuming for complex
systems (e.g., for a high number of elements), this approach could be used for
sub-models to highlight the influencing factors and importance of GNAs. While
in the presented example, the correlation was introduced artificially, it
occurs naturally in more complex systems (e.g., strength properties are
dependent on the material chosen, or, if giving an example outside engineering
applications, different price ranges lead to different demand levels, etc.).
Projecting these findings to a larger scope of applications,an increase in the
accuracy of structural analyses for SLS estimates would improve the overall
reliability and integrity of the structures. Consequently, costly required
changes to structural elements potentially found post hoc during
service/operation could be avoided via proper a priori engineering analysis
and design.

##  4 Relationships between second-order indices and correlation

The mechanical engineering model demonstrated that in the presence of positive
correlation, the second-order sensitivity index turns negative (Table 5). But
how does the index behave when the correlation is negative? Furthermore, is
the relationship the same in the presence of an interaction of inputs that
affect the output?

To explore these questions, we set up a series of simulation experiments for
two simple two-factor models. The first model is additive, $Y=A+B$, where the
variables do not interact. The second model is multiplicative, $Y=AB$, which
possesses a positive second-order effect when the variables are independent,
thereby manifesting interaction. We also examine two different types of
correlation, (i) joint multivariate uniform distribution copula and (ii)
equating values of $B$ to $A$ over a portion of the sample, which might be
considered as a structural change in the system. Both correlation types are
modeled for different correlation strengths. For validation purposes, Pearson
and Spearman correlation coefficients are reported in each case. Both $A$ and
$B$ assume uniform distribution between $0$ and $5$. The indices are
calculated based on a sample size of $10^{6}$ to ensure a high level of
precision.

The sensitivity indices estimates are presented in Table 6, and the
relationship between the second-order indices and the Pearson correlation
coefficient is visualized in Figure 8. Hart and Gremaud (2018) have indicated
that theoretical second-order indices can become negative in the presence of
dependent inputs and such an outcome can be observed for the approximated
indices computed by our simple binning algorithm (Table 6).

Table 6:  Sensitivity indices and correlation coefficients for a simple $Y=AB$ model with $A$ and $B$ dependent to different degrees. Correlation type | Model | $S$ | Correlation
---|---|---|---
negative |  | positive
-100 % | … | -50 % | … | 0 % | … | 50 % | … | 100 %
| $A+B$ | $S_{A}$ | - | 0.14 | 0.27 | 0.38 | 0.50 | 0.62 | 0.74 | 0.87 | 1.00
$S_{B}$ | - | 0.14 | 0.27 | 0.38 | 0.50 | 0.62 | 0.74 | 0.87 | 1.00
$S_{AB}$ | - | 0.71 | 0.47 | 0.23 | 0.00 | -0.24 | -0.49 | -0.74 | -1.00
$\sum$ | - | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
$AB$ | $S_{A}$ | 1.00 | 0.31 | 0.27 | 0.33 | 0.43 | 0.56 | 0.70 | 0.85 | 1.00
$S_{B}$ | 1.00 | 0.31 | 0.27 | 0.33 | 0.43 | 0.56 | 0.70 | 0.85 | 1.00
$S_{AB}$ | -1.00 | 0.39 | 0.47 | 0.35 | 0.14 | -0.11 | -0.40 | -0.70 | -1.00
$\sum$ | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Pearson | -1.00 | -0.73 | -0.48 | -0.24 | 0.00 | 0.24 | 0.48 | 0.73 | 1.00
Spearman | -1.00 | -0.73 | -0.48 | -0.24 | 0.00 | 0.24 | 0.48 | 0.73 | 1.00
| $A+B$ | $S_{A}$ | - | 0.61 | 0.36 | 0.45 | 0.50 | 0.65 | 0.79 | 0.93 | 1.00
$S_{B}$ | - | 0.04 | 0.26 | 0.47 | 0.50 | 0.67 | 0.76 | 0.83 | 1.00
$S_{AB}$ | - | 0.35 | 0.39 | 0.09 | 0.00 | -0.32 | -0.56 | -0.76 | -1.00
$\sum$ | - | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
$AB$ | $S_{A}$ | 1.00 | 0.52 | 0.39 | 0.42 | 0.43 | 0.44 | 0.51 | 0.74 | 1.00
$S_{B}$ | 1.00 | 0.09 | 0.17 | 0.35 | 0.43 | 0.59 | 0.86 | 0.93 | 1.00
$S_{AB}$ | -1.00 | 0.39 | 0.44 | 0.24 | 0.14 | -0.03 | -0.36 | -0.67 | -1.00
$\sum$ | 1.00 | 1.00 | 1.00 | 1.01 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Pearson | -1.00 | -0.71 | -0.52 | -0.23 | 0.00 | 0.23 | 0.52 | 0.71 | 1.00
| Spearman | -1,00 | -0.71 | -0.53 | -0.24 | 0.00 | 0.24 | 0.53 | 0.71 | 1.00
Figure 8: Second-order effect as a function of correlation for additive (0%
interaction) and multiplicative (14% interaction) functions with two
correlation types, (1) - portion of equal values, (2) - copula.

Additive models cross the origin (see Figure 8) and possess zero second-order
effects under a no-dependency case. No interaction is anticipated from the
additive model. Copula variables show the second-order effect values very
close to the correlation coefficients and form a linear relationship between
the two. Positive correlation results in negative second-order effects,
conveying overlapping effects of the variables. Conversely, negative
correlation results in a positive second-order effect for additive models as
if the combined effect of the variables is greater than the sum of their
individual effects (see equation 4). The additive model, in which dependence
is modeled through a portion of equal values, has largely the same pattern
except for lower positive second-order effects in the presence of negative
correlation. The calculation of significance indices for both additive models
at -100% correlation fails, because all $Y$ values turn to $0$. Instead, the
values for the negative 95% correlation are computed and displayed in Figure
8.

In the multiplicative model, the second-order effect index behaves
asymmetrically for positive and negative correlation (Figure 8). With
increasing positive correlation, the second-order index becomes increasingly
negative (as in the additive models), signifying the overlapping effects of
the input variables on the output. For negative correlations, the second-order
index is initially positive and increasing, peaking at 50% correlation. This
translates into a synergistic effect where the negatively correlated inputs
magnify the impact on the output. But when the negative correlation exceeds
50%, the second-order index decreases and approaches $-1$. One could speculate
that with a strong negative correlation for interacting variables, the joint
overlapping effect becomes stronger than the synergistic one, driving the
second-order index to $-1$.

For all cases, the conservation property (the sum of all indices is equal to
1) holds. The boundedness property for each index extends from $[0,1]$, as for
the classic Sobol’ indices, to $[-1,1]$. Because the conservation property
holds under all cases, this implies that correlation and interaction effects
are actually both combined into the single second-order effect, and,
simultaneously affecting the first-order index behavior, perfectly fitting
into the variance decomposition of the output. The presence of interaction can
be observed in Figure 8 as the vertical distance between the additive and
multiplicative models’ lines on the correlation range between $-0.25$ and
$0.25$. In other words, the different degrees of dependency move the second-
order effect up or down, but the presence of interaction reveals itself in
higher second-order effect values.

##  5 Visualizing second-order effects – Beyond sensitivity indices

The sensitivity analysis field has focused on computing sensitivity indices to
be able to rank the input variables by their relative importance (Saltelli et
al., 2008; Da Veiga et al., 2021). And yet it becomes apparent that
sensitivity indices are not able to convey the full relationship story. As
observed in our earlier analysis with dependent input variables (Table 6), the
second-order indices for the same model, but operating under different
assumptions, appear to be very similar. However, the scatter plots (Table 6)
reveal how different these relationships actually are, and this knowledge can
lead to different structural design decisions.

A fact that has been largely neglected throughout modern sensitivity analysis
studies is that although specific summary values can appear near-identical in
magnitude, an examination of the underlying data can reveal very different
shapes and perspectives (Puy, Beneventano, et al., 2022; Saltelli et al.,
2004; Kozlova et al., in press).

Consequently, different kinds of visualizations can be used to augment and
improve the understanding of the underlying model. These data visualizations
can include such formats as basic scatter plots, response surfaces, heat maps,
and parallel coordinate plots. However, all of these representations possess
dimensionality limitations. SimDec is a recent approach that is able to
project multidimensional systems onto a two-dimensional graph (Kozlova and
Yeomans, 2022b) and, in combination with the simple binning sensitivity
analysis approach, can be used to convey the nature of interactions in complex
engineering models to produce a much more comprehensive sensitivity analysis
framework (Kozlova et al., in press).

To illustrate the power of this framework, this section showcases the
application of the simple binning approach together with SimDec visualization
for investigating the behavior of a complex structural reliability model
containing both correlations and interactions.

###  5.1 Case background: structural reliability model for fatigue assessments

Many structural applications experience cyclic or fluctuating load conditions
during service life (such as ship and offshore structures as well as various
vehicles, trucks and trains). In these complex systems, fatigue is amongst the
most important design criteria to control and maintain the lifespan of ageing
components without concerns of structural failures. Due to the random nature
of acting cyclic loads in many applications and the complexity of physical
modeling of fatigue phenomenon, the vast majority of existing models for
predicting fatigue life are necessarily based on probabilistic or data-driven
approaches (Bai et al., 2023; Ye et al., 2023; Ruiz Muñoz and Sorensen, 2020;
Leonetti et al., 2020). However, current trends for improving material usage
and optimization in mechanical components have led to an ever-increasing need
for creating more sophisticated and more accurate deterministic models to
assess the fatigue life of components. This is particularly the case for
modern infrastructures that utilize high-strength steels, due to the fact that
an increase in material strength does not contribute to higher fatigue
performance in welded components (Lieurade et al., 2008).

To account for the effects of material strength and welding quality, a multi-
parametric model – named the 4R method – was introduced for fatigue assessment
of welded components. The 4R model incorporates four parameters - material
strength, welding residual stresses, applied stress ratio, and welding quality
(Ahola et al., 2021, 2020). These four key parameters affect the fatigue
strength of welded components depending on the conditions (Hobbacher, 2016).
As the 4R model employs more parameters than conventional models, a key
question still rests on a determination of the sensitivity of each parameter
on the 4R output value.

The basis of the 4R method is to utilize the effective stress concept, usually
known as the effective notch stress (ENS) approach (Sonsino et al., 2012), in
association with the mean stress correction using the requisite four
parameters. In the 4R model, material elastic-plastic behavior at the notch
root is analytically or numerically computed, after which the local hysteresis
loops of cyclic stress are determined. The output reference stress applied in
fatigue assessments (in the form of S-N curve) is computed based on the well-
known Smith-Watson-Topper (Smith et al., 1970) mean stress correction on the
linear-elastic notch stress (equation 8). Figure 9 demonstrates the workflow
of the 4R method.

$\ \Delta\sigma_{k,ref}=\frac{\Delta\sigma_{k}}{\sqrt{1-R_{local}}}\\\ $ (8)

Figure 9: Description and workflow of the 4R model for fatigue assessment.

###  5.2 Sensitivity indices for the structural reliability model

The 4R model is simulated $10^{4}$ times with inputs uniformly distributed
within the ranges specified in Table 7.

Table 7:  Assumptions on input variation for the structural reliability 4R model. Case material |  | Weld toe radius,
---
$r_{true}$ (mm)
| Steel grade,
---
$R_{p0.2}$ (MPa)
| Residual stress,
---
$\sigma_{res}$ (MPa)
| Stress ratio,
---
$R$ (-)
S355 | 0.5 ± 0.5 | 355 ± 100 | 250 ± 100 | 0.5 ± 0.2
0 ± 100 | 0 ± 0.2
-100 ± 100 | -1 ± 0.2
S960 | 0.5 ± 0.5 | 960 ± 100 | 850 ± 100 | 0.5 ± 0.2
0 ± 100 | 0 ± 0.2
-300 ± 100 | -1 ± 0.2

The three levels of residual stress and stress ratio for each case material
are formed based on actual conditions in structural components. Usually
welding causes high tensile residual stresses, up to the yield strength of
parent material. On the other hand, residual stresses can be reduced by, for
example, post-heat treatments and/or mechanical post-weld treatments or
residual stress relaxation via the mechanical load. The applied stress ratio
is specific to load conditions in service, and the selected three cases
represent available cases in engineering components. The weld toe radius is
converted in the model to fatigue-effective stress, $K_{f}$, the sensitivity
to which is analyzed. The three other inputs used in the sensitivity analysis,
as is. The sensitivity indices to the model output, $\Delta\sigma_{k,ref}$,
computed with the simple binning method, are presented in Table 8.

Table 8:  Sensitivity indices for the structural reliability 4R model. Variable |  First-order effect | Second-order effect |  Combined effect
---|---|---|---
| Residual
---
stress,
$\sigma_{res}$
| Stress
---
ratio,
$R$
| Steel
---
grade,
$R_{p0.2}$
| Fatigue-
---
effective
stress, $K_{f}$
Residual stress, $\sigma_{res}$ | 50% |  | 11% | $-6\%$ | 0% | 51%
Stress ratio, $R$ | 28% |  |  | 4% | 0% | 35%
Steel grade, $R_{p0.2}$ | 11% |  |  |  | 0% | 10%
Fatigue-effective stress, $K_{f}$ | 4% |  |  |  |  | 4%
$\Sigma$ | 92% |  |  |  |  | 100%

The sensitivity indices in Table 8 show that the residual stress
$\sigma_{res}$, alone, explains half of the variance of the output. Residual
stress $\sigma_{res}$ also has $11\%$ interaction with stress ratio $R$, and a
negative $-6\%$ second-order effect with steel grade $R_{p0.2}$. In addition,
stress ratio $R$ and steel grade $R_{p0.2}$ have a $4\%$ interaction. The
remaining second-order effects are zero. Thus, there are three input variables
in the model that have pair-wise interaction effects with each other: residual
stress $\sigma_{res}$, stress ratio $R$, and steel grade $R_{p0.2}$. Together
these variables explain $96\%$ of the output variability. However, as
mentioned above, the sensitivity indices alone cannot adequately explain the
nature of those interactions and how they should be accounted for in decision-
making processes.

###  5.3 Exploring the nature of second-order effects

To investigate how the model output is affected simultaneously by the three
input variables identified as possessing pairwise interaction effects, we
employ SimDec to provide a multidimensional visualization perspective of their
impacts (Kozlova and Yeomans, 2022a; Kozlova et al., in press). Since SimDec
uses a frequency distribution, the same simulation data that was used to
compute sensitivity indices with the simple binning method, can be used for
the visualization. The idea behind SimDec is to decompose the simulation data
into regions formed by the combinations of different ranges (or states) of
influential input variables. An intelligent color-coding of these regions in
the frequency distribution communicates a clear visual representation of the
inherent input-output relationships within the model (for the algorithm
details, see (Kozlova et al., in press)). The open-source SimDec code is
available for Python, R, Julia, and Matlab (Kozlova et al., 2023).

For the structural reliability model, we choose to break down the residual
stress $\sigma_{res}$ into three states (low, medium, high), the stress ratio
$R$ into two states (reversed, pulsating), and the steel grade $R_{p0.2}$ into
two states (mild, UHSS) (Table 9).

Table 9:  Decomposition set-up for the structural reliability 4R model (UHSS stands for ultra-high steel strength). Residual stress |  | Stress ratio |  | Steel grade
---|---|---|---|---
State | min | max |  | State | min | max |  | State | min | max
Low | $-400$ | 100 |  | Reversed | $-1.2$ | $-0.25$ |  | Mild | 255 | 657
Medium | 100 | 650 |  | Pulsating | $-0.25$ | 0.7 |  | UHSS | 657 | 1060
High | 650 | 950 |  | - | - | - |  | - | - | -

Altogether, all possible combinations of these state settings for the input
variables (Table 9) form 12 separate scenarios. The probability distribution
of the model output, $\Delta\sigma_{k,ref}$, is partitioned using these
scenarios and color-coded. The coloring logic is important to facilitate the
overall interpretability of the visual perceptions. The states of the most
influential input variable assume distinct main colors, and each of these main
colors is sub-shaded to further highlight the partitions, Figure 10.

Figure 10: Simulation decomposition of the structural reliability 4R model.
---
Output Summary Color Scenario Residual stress, $\sigma_{res}$ Stress ratio,
$R$ Steel grade, $R_{p0.2}$ min mean max probability sc1 Mild $107$ $332$
$475$ $11\%$ sc2 Reversed UHSS $11$ $264$ $515$ $11\%$ sc3 Mild $320$ $467$
$591$ $22\%$ sc4 Low Pulsating UHSS $68$ $543$ $819$ $22\%$ sc5 Mild $383$
$456$ $524$ $6\%$ sc6 Reversed UHSS NaN NaN NaN NaN sc7 Mild $422$ $499$ $622$
$12\%$ sc8 Medium Pulsating UHSS NaN NaN NaN NaN sc9 Mild NaN NaN NaN NaN sc10
Reversed UHSS $630$ $704$ $795$ $6\%$ sc11 Mild NaN NaN NaN NaN sc12 High
Pulsating UHSS $656$ $746$ $851$ $12\%$ Table 10: Summary of simulation
decomposition of the structural reliability 4R model.

Figure 10 demonstrates how non-monotonic the effects of the input variables on
the model output are, and that the effect of one input is conditioned to the
states of another one. First, the residual stress $\sigma_{res}$ divides the
distribution into two narrow scenarios, medium and high, whereas the low
scenario sub-distribution spreads out through the entire range of the output
values. Stress ratio $R$ has only a minor effect on the output in medium and
high states of residual stress (small horizontal shift of the respective sub-
distributions) but plays an important role in the low residual stress, almost
dividing the output into halves. In low residual stress, steel grade
$R_{p0.2}$ has a minor effect when the stress ratio is reversed but creates a
substantial rightward shift if the stress ratio is pulsating. The combinations
of medium residual stress & UHSS steel grade and high residual stress & mild
steel grade are non-existent, which reflects physical material constraints and
thereby explains the negative second-order effect. These visceral
relationships would not have been revealed or apparent without the SimDec
visualization.

##  6 Conclusions

This paper introduces and extensively tests the simple binning approach for
computing first- and second-order sensitivity indices, which provide more
accurate results than classic in global sensitivity analysis Sobol’ indices.
Further, the paper demonstrates the importance of visualizing input-output
relationships on a structural reliability model using a SimDec approach. The
resulting methodological framework is streamlined (i.e., it works with a given
data of as little as 1000 data points), open-sourced, and provides an
intuitive way of ’looking’ into a computational model and grasping its
behavior. The presented framework can be recommended for analyzing disparate
complex technological models for engineering design and decision-making.

## Acknowledgements

The authors are grateful to Prof. Art Owen for his comments in the early
stages of this work. The work was supported by grant 220178 from the Finnish
Foundation for Economic Education, by project CaNeLis, 3406/31/2022, funded by
Business Finland, and by grant OGP0155871 from the Natural Sciences and
Engineering Research Council.

## References

  *   * Ahola et al. (2021) Ahola, A., Muikku, A., Braun, M., and Björk, T.  (2021).  Fatigue strength assessment of ground fillet-welded joints using 4r method.  _International Journal of Fatigue_ , _142_ , 105916.
  * Ahola et al. (2020) Ahola, A., Skriko, T., and Björk, T.  (2020).  Fatigue strength assessment of ground fillet-welded joints using 4r method.  _Journal of Constructional Steel Research_ , _167_ , 105861.
  * Alidoosti et al. (2013) Alidoosti, A., Ghafari-Nazari, A., Moztarzadeh, F., Jalali, N., Moztarzadeh, S., and Mozafari, M.  (2013).  Electrical discharge machining characteristics of nickel–titanium shape memory alloy based on full factorial design.  _Journal of intelligent Material systems and Structures_ , _24_(13), 1546–1556.
  * Austrell et al. (2004) Austrell, P.-E., Dahlbom, O., Lindemann, J., Olsson, A., Olsson, K.-G., et al.  (2004).  Calfem - a finite element toolbox, version 3.4.  _Lund, Sweden_.
  * Bai et al. (2023) Bai, S., Huang, T., Li, Y.-F., Lu, N., and Huang, H.-Z.  (2023).  A probabilistic fatigue life prediction method under random combined high and low cycle fatigue load history.  _Reliability Engineering & System Safety_, _238_ , 109452.
  * Ballester-Ripoll and Leonelli (2022) Ballester-Ripoll, R., and Leonelli, M.  (2022).  Computing sobol indices in probabilistic graphical models.  _Reliability Engineering & System Safety_, _225_ , 108573.
  * Barr and Rabitz (2023) Barr, J., and Rabitz, H.  (2023).  Kernel-based global sensitivity analysis obtained from a single data set.  _Reliability Engineering & System Safety_, _235_ , 109173.
  * Borgonovo (2007) Borgonovo, E.  (2007).  A new uncertainty importance measure.  _Reliability Engineering & System Safety_, _92_(6), 771–784.
  * Borgonovo et al. (2023) Borgonovo, E., Plischke, E., and Prieur, C.  (2023).  Total effects with constrained features.
  * Caers (2018) Caers, J.  (2018).  _04-2 sensitivity analysis global._ Lecture, Stanford University course ”QUSS GS 260”.  (Available at https://youtu.be/vBuWB9WuFhA)
  * Chalabi et al. (2023) Chalabi, Y., Dutang, C., Savicky, P., Wuertz, D., Knuth, D., Matsumoto, M., and Saito, M.  (2023).  _Toolbox for pseudo and quasi random number generation and random generator tests._ https://cran.r-project.org/web/packages/randtoolbox/randtoolbox.pdf.
  * Crestaux et al. (2009) Crestaux, T., Le Maıtre, O., and Martinez, J.-M.  (2009).  Polynomial chaos expansion for sensitivity analysis.  _Reliability Engineering & System Safety_, _94_(7), 1161–1172.
  * Cukier et al. (1973) Cukier, R., Fortuin, C., Shuler, K. E., Petschek, A., and Schaibly, J. H.  (1973).  Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. i theory.  _The Journal of chemical physics_ , _59_(8), 3873–3878.
  * Da Veiga (2015) Da Veiga, S.  (2015).  Global sensitivity analysis with dependence measures.  _Journal of Statistical Computation and Simulation_ , _85_(7), 1283–1305.
  * Da Veiga et al. (2021) Da Veiga, S., Gamboa, F., Iooss, B., and Prieur, C.  (2021).  _Basics and trends in sensitivity analysis: Theory and practice in r_.  SIAM.
  * Deviatkin et al. (2021) Deviatkin, I., Kozlova, M., and Yeomans, J. S.  (2021).  Simulation decomposition for environmental sustainability: Enhanced decision-making in carbon footprint analysis.  _Socio-Economic Planning Sciences_ , _75_ , 100837.
  * Hart and Gremaud (2018) Hart, J., and Gremaud, P. A.  (2018).  An approximation theoretic perspective of sobol’indices with dependent variables.  _International Journal for Uncertainty Quantification_ , _8_(6).
  * Hobbacher (2016) Hobbacher, A. F.  (2016).  _Recommendations for fatigue design of welded joints and components_.  Springer: Cham, Switzerland.
  * Homma and Saltelli (1996) Homma, T., and Saltelli, A.  (1996).  Importance measures in global sensitivity analysis of nonlinear models.  _Reliability Engineering & System Safety_, _52_(1), 1–17.
  * Iooss et al. (2022) Iooss, B., Sudret, B., Piano, L., and Prieur, C.  (2022).  Editorial for the special issue on “sensitivity analysis of model outputs” reliability engineering and system safety.  _Reliability Engineering & System Safety_, _223_ , 108477.
  * Jung and Taflanidis (2023) Jung, W., and Taflanidis, A. A.  (2023).  Efficient global sensitivity analysis for high-dimensional outputs combining data-driven probability models and dimensionality reduction.  _Reliability Engineering & System Safety_, _231_ , 108805.
  * Kozlova et al. (in press) Kozlova, M., Moss, R. J., Yeomans, J. S., and Caers, J.  (in press).  Uncovering heterogeneous effects in computational models for sustainable decision-making.  _Environmental Modelling and Software_.  DOI: 10.2139/ssrn.4550911
  * Kozlova et al. (2023) Kozlova, M., Roy, P., Moss, R. J., and Alam, A.  (2023).  _Simdec repositories in python, r, julia, and matlab._ https://github.com/Simulation-Decomposition.
  * Kozlova and Yeomans (2019) Kozlova, M., and Yeomans, J.  (2019).  Multi-variable simulation decomposition in environmental planning: An application to carbon capture and storage.  _J. Environ. Inform. Lett_ , _1_ , 20–26.
  * Kozlova and Yeomans (2022a) Kozlova, M., and Yeomans, J.  (2022a).  Extending simulation decomposition analysis into systemic risk planning for domino-like cascading effects in environmental systems.  _J. Environ. Informatics Lett_ , _7_ , 64–68.
  * Kozlova and Yeomans (2022b) Kozlova, M., and Yeomans, J. S.  (2022b).  Monte carlo enhancement via simulation decomposition: A “must-have” inclusion for many disciplines.  _INFORMS Transactions on Education_ , _22_(3), 147–159.
  * Kucherenko and Song (2016) Kucherenko, S., and Song, S.  (2016).  Derivative-based global sensitivity measures and their link with sobol’sensitivity indices.  In _Monte carlo and quasi-monte carlo methods: Mcqmc, leuven, belgium, april 2014_ (pp. 455–469).
  * Lamboni and Kucherenko (2021) Lamboni, M., and Kucherenko, S.  (2021).  Multivariate sensitivity analysis and derivative-based global sensitivity measures with dependent variables.  _Reliability Engineering & System Safety_, _212_ , 107519.
  * Leonetti et al. (2020) Leonetti, D., Majlaars, J., and Snijder, H. H.  (2020).  Probabilistic fatigue resistance model for steel welded details under variable amplitude loading – inference and uncertainty estimation.  _International Journal of Fatigue_ , _135_ , 105515.
  * Lieurade et al. (2008) Lieurade, H. P., Huther, I., and Lefebvre, F.  (2008).  Effect of weld quality and postweld improvement techniques on the fatigue resistance of extra high strength steels.  _Welding in the World_ , _52_ , 106–115.
  * Liu et al. (2022) Liu, Y.-C., Leifsson, L., Pietrenko-Dabrowska, A., and Koziel, S.  (2022).  Analysis of agricultural and engineering systems using simulation decomposition.  In _International conference on computational science_ (pp. 435–444).
  * Mahadevan and Ni (2003) Mahadevan, S., and Ni, K.  (2003).  Damage tolerance reliability analysis of automotive spot-welded joints.  _Reliability Engineering & System Safety_, _81_ , 9–21.
  * Mara et al. (2015) Mara, T. A., Tarantola, S., and Annoni, P.  (2015).  Non-parametric methods for global sensitivity analysis of model output with dependent inputs.  _Environmental modelling & software_, _72_ , 173–183.
  * Marzban and Lahmer (2016) Marzban, S., and Lahmer, T.  (2016).  Conceptual implementation of the variance-based sensitivity analysis for the calculation of the first-order effects.  _Journal of Statistical Theory and Practice_ , _10_(4), 589–611.
  * MathWorks (2013) MathWorks.  (2013).  _Generating quasi-random numbers._ https://se.mathworks.com/help/stats/generating-quasi-random-numbers.html.
  * Myers et al. (2016) Myers, R. H., Montgomery, D. C., and Anderson-Cook, C. M.  (2016).  _Response surface methodology: process and product optimization using designed experiments_.  John Wiley & Sons.
  * Olken and Rotem (1986) Olken, F., and Rotem, D.  (1986).  Simple random sampling from relational databases.
  * A. B. Owen (2020) Owen, A. B.  (2020).  On dropping the first sobol’point.  In _International conference on monte carlo and quasi-monte carlo methods in scientific computing_ (pp. 71–86).
  * A. B. Owen (2023) Owen, A. B.  (2023).  _Practical Quasi-Monte Carlo_.  https://artowen.su.domains/mc/practicalqmc.pdf.
  * M. P. Owen et al. (2019) Owen, M. P., Panken, A., Moss, R., Alvarez, L., and Leeper, C.  (2019).  Acas xu: Integrated collision avoidance and detect and avoid capability for uas.  In _2019 ieee/aiaa 38th digital avionics systems conference (dasc)_ (pp. 1–10).
  * Pleil et al. (2011) Pleil, J. D., Stiegel, M. A., Madden, M. C., and Sobus, J. R.  (2011).  Heat map visualization of complex environmental and biomarker measurements.  _Chemosphere_ , _84_(5), 716–723.
  * Plischke (2012) Plischke, E.  (2012).  How to compute variance-based sensitivity indicators with your spreadsheet software.  _Environmental Modelling & Software_, _35_ , 188–191.
  * prEN 1993-1-14 (2023) prEN 1993-1-14.  (2023).  _Eurocode 3 - design of steel structures - part -1-14: Design assisted by finite element analysis_.  European Committee for Standardization.
  * Puy, Beneventano, et al. (2022) Puy, A., Beneventano, P., Levin, S. A., Lo Piano, S., Portaluri, T., and Saltelli, A.  (2022).  Models with higher effective dimensions tend to produce more uncertain estimates.  _Science Advances_ , _8_(42), eabn9450.
  * Puy, Roy, and Saltelli (2022) Puy, A., Roy, P., and Saltelli, A.  (2022).  Discrepancy measures for sensitivity analysis in mathematical modeling.  _arXiv preprint arXiv:2206.13470_.
  * Razavi and Gupta (2016) Razavi, S., and Gupta, H. V.  (2016).  A new framework for comprehensive, robust, and efficient global sensitivity analysis: 1. theory.  _Water Resources Research_ , _52_(1), 423–439.
  * Robbe (2018) Robbe, P.  (2018).  _The qmc module for julia._ https://github.com/PieterjanRobbe/QMC.jl.
  * Ruiz Muñoz and Sorensen (2020) Ruiz Muñoz, S., and Sorensen, J.  (2020).  Probabilistic inspection planning of offshore welds subject to the transition from protected to corrosive environment.  _Reliability Engineering & System Safety_, _202_ , 107009.
  * Saltelli (2023) Saltelli, A.  (2023, June).  _Sensitivity analysis made easy._ Presentation at Workshop on Modelling.  (Available at http://www.andreasaltelli.eu/file/repository/Talk_SA_June_28_BIS.pdf)
  * Saltelli et al. (2019) Saltelli, A., Aleksankina, K., Becker, W., Fennell, P., Ferretti, F., Holst, N., … Wu, Q.  (2019).  Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices.  _Environmental modelling & software_, _114_ , 29–39.
  * Saltelli et al. (2021) Saltelli, A., Jakeman, A., Razavi, S., and Wu, Q.  (2021).  Sensitivity analysis: A discipline coming of age.  _Environmental Modelling & Software_, _146_ , 105226.
  * Saltelli et al. (2008) Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., … Tarantola, S.  (2008).  _Global sensitivity analysis: the primer_.  John Wiley & Sons.
  * Saltelli et al. (2004) Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M., et al.  (2004).  Sensitivity analysis in practice: a guide to assessing scientific models.  _Chichester, England_.
  * SciPy (2023) SciPy.  (2023).  _Quasi-monte carlo submodule (scipy.stats.qmc)._ https://docs.scipy.org/doc/scipy/reference/stats.qmc.html.
  * Shang et al. (2023) Shang, X., Su, L., Fang, H., Zeng, B., and Zhang, Z.  (2023).  An efficient multi-fidelity kriging surrogate model-based method for global sensitivity analysis.  _Reliability Engineering & System Safety_, _229_ , 108858.
  * Shapley et al. (1953) Shapley, L. S., et al. (1953).  A value for n-person games.
  * Shen et al. (2023) Shen, X., Feng, K., Xu, H., Wang, G., Zhang, Y., Dai, Y., and Yun, W.  (2023).  Reliability analysis of bending fatigue life of hydraulic pipeline.  _Reliability Engineering & System Safety_, _231_ , 109019.
  * Shi et al. (2023) Shi, W., Zhou, Q., and Zhou, Y.  (2023).  An efficient elementary effect-based method for sensitivity analysis in identifying main and two-factor interaction effects.  _Reliability Engineering & System Safety_, _237_ , 109365.
  * Shittu et al. (2021) Shittu, A. A., Mehmanparast, A., and Hart, P.  (2021).  Comparative study between s-n and fracture mechanics approach on reliability assessment of offshore wind turbine jacket foundations.  _Reliability Engineering & System Safety_, _215_ , 107838.
  * Singh and Singh (2003) Singh, S., and Singh, S.  (2003).  Simple random sampling.  _Advanced Sampling Theory with Applications: How Michael ‘selected’Amy Volume I_ , 71–136.
  * Smith et al. (1970) Smith, K. N., Watson, P., and Topper, T. H.  (1970).  A stress-strain function for the fatigue of metals.  _Journal of Materials_ , _5_ , 767–778.
  * Sobol and Kucherenko (2010) Sobol, I., and Kucherenko, S.  (2010).  Derivative based global sensitivity measures.  _Procedia-Social and Behavioral Sciences_ , _2_(6), 7745–7746.
  * Sobol’ (1967) Sobol’, I. M.  (1967).  On the distribution of points in a cube and the approximate evaluation of integrals.  _Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki_ , _7_(4), 784–802.
  * Sobol (1993) Sobol, I. M.  (1993).  Sensitivity analysis for non-linear mathematical models, originally “sensitivity estimates for non-linear mathematical models”.  _Math Model Comput Exp_ , _1_ , 407–414.
  * Sonsino et al. (2012) Sonsino, C. M., Fricke, W., de Bruyne, F., Hoppe, A., Ahmadi, A., and Zhang, G.  (2012).  Notch stress concepts for the fatigue assessment of welded joints – background and applications.  _International Journal of Fatigue_ , _34_ , 2–16.
  * Suard et al. (2013) Suard, S., Hostikka, S., and Baccou, J.  (2013).  Sensitivity analysis of fire models using a fractional factorial design.  _Fire safety journal_ , _62_ , 115–124.
  * Sudret (2008) Sudret, B.  (2008).  Global sensitivity analysis using polynomial chaos expansions.  _Reliability engineering & system safety_, _93_(7), 964–979.
  * Thapa and Missoum (2022) Thapa, M., and Missoum, S.  (2022).  Uncertainty quantification and global sensitivity analysis of composite wind turbine blades.  _Reliability Engineering & System Safety_, _222_ , 108354.
  * Tong (2006) Tong, C.  (2006).  Refinement strategies for stratified sampling methods.  _Reliability Engineering & System Safety_, _91_(10-11), 1257–1265.
  * Vuillod et al. (2023) Vuillod, B., Montemurro, M., Panettieri, E., and Hallo, L.  (2023).  A comparison between sobol’s indices and shapley’s effect for global sensitivity analysis of systems with independent input variables.  _Reliability Engineering & System Safety_, _234_ , 109177.
  * Wang and Jia (2023) Wang, Z., and Jia, G.  (2023).  Extended sample-based approach for efficient sensitivity analysis of group of random variables.  _Reliability Engineering & System Safety_, _231_ , 108991.
  * Xiong et al. (2022) Xiong, Q., Du, P., Deng, J., Huang, D., Song, G., Qian, L., … Luo, Y.  (2022).  Global sensitivity analysis for nuclear reactor lbloca with time-dependent outputs.  _Reliability Engineering & System Safety_, _221_ , 108337.
  * Ye et al. (2023) Ye, Y., Yang, Q., Zhang, J., Meng, S., and Wang, J.  (2023).  A dynamic data driven reliability prognosis method for structural digital twin and experimental validation.  _Reliability Engineering & System Safety_, _240_ , 109543.
  * Ziehn and Tomlin (2009) Ziehn, T., and Tomlin, A. S.  (2009).  Gui–hdmr–a software tool for global sensitivity analysis of complex models.  _Environmental Modelling & Software_, _24_(7), 775–785.