Sensitivity analysis aims at assessing the relative contributions of the different sources of uncertainty to the variability in the output of a model. There are several mathematical techniques available in the literature, with variance-based approaches being the most popular, and variance-based indices the most widely used, especially “total effects” indices. Literature also provides several estimators/approximators for these indices (reviews can be found here [1]), which typically need *N* = *n* × (*M* + 2) model evaluations (and resulting outputs), where *M* is the number of uncertain inputs and *n* is some factor of proportionality that is selected ad hoc, depending on the complexity of the model (e.g. linear or not, monotonic or not, continuous or not). [Note: Literature typically refers to *n* as the *number of samples* and to *N* as the *number of model evaluations*, and this is how I’ll be using them also.]

The generation of this set of model runs of size N is by far the most computationally demanding step in the calculation of variance-based indices, and a major limitation to their applicability, as it’s typically in the order of magnitude of thousands [2] (*n* typically >1000). For example, consider a model with 5 uncertain inputs that takes 1 minute to run. To appropriately estimate sensitivity indices for such a model, we would potentially need about N=7000 model evaluations, which would take almost 5 days on a personal computer, excluding the time for the estimator calculation itself (which is typically very short).

The aim is therefore to pick the minimum *n* needed to ensure our index calculation is reliable. Unfortunately, there is no hard and fast rule on how to do that, but the aim of this post is to illuminate that question a little bit and provide some guidance. I am going off the work presented here [3] and elsewhere, and **the aim is to perform the estimation of sensitivity indices repeatedly, using an increasing number of n until the index values converge.**

I will be doing this using a fishery model, which is a nonlinear and nonmonotonic model with 9 parameters. Based on previous results suggesting that 3 of these parameters are largely inconsequential to the variability in the output, I’ll be fixing them to their nominal values. I’ll be using SALib to perform the analysis. My full script can be found here, but I’ll only be posting the most important snippets of code in this post.

First, I set up my SALib ‘problem’ and create arrays to store my results:

```
# Set up dictionary with system parameters
problem = {
'num_vars': 6,
'names': ['a', 'b', 'c','h',
'K','m'],
'bounds': [[ 0.002, 2],
[0.005, 1],
[0.2, 1],
[0.001, 1],
[100, 5000],
[0.1, 1.5]]
}
# Array with n's to use
nsamples = np.arange(50, 4050, 50)
# Arrays to store the index estimates
S1_estimates = np.zeros([problem['num_vars'],len(nsamples)])
ST_estimates = np.zeros([problem['num_vars'],len(nsamples)])
```

I then loop through all possible *n* values and perform the sensitivity analysis:

```
# Loop through all n values, create sample, evaluate model and estimate S1 & ST
for i in range(len(nsamples)):
print('n= '+ str(nsamples[i]))
# Generate samples
sampleset = saltelli.sample(problem, nsamples[i],calc_second_order=False)
# Run model for all samples
output = [fish_game(*sampleset[j,:]) for j in range(len(sampleset))]
# Perform analysis
results = sobol.analyze(problem, np.asarray(output), calc_second_order=False,print_to_console=False)
# Store estimates
ST_estimates[:,i]=results['ST']
S1_estimates[:,i]=results['S1']
```

I can then plot the evolution of these estimates as *n* increases:

What these figures tell us is that choosing an *n* below 1000 for this model would potentially misestimate our indices, especially the first order ones (S1). As *n* increases, we see the estimates becoming less noisy and converging to their values. For more complex models, say, with more interactive effects, the minimum *n* before convergence could actually be a lot higher. A similar experiment by Nossent et al. [3] found that convergence was reached only after *n=*12,000*. *

An observation here is that the values of the total indices (ST) are higher than those of the the first order indices (S1), which makes sense, as ST includes both first order effects (captured by S1) and second order effects (i.e. interactions between the factors). Another observation here is that the parameters with the most significant effects (*m *and *K*) converge much faster than parameters with less impact on the output (*a *and *b*). This was also observed by Nossent et al. [3].

Finally, sensitivity analysis is often performed for the purposes of factor prioritization, i.e. determining (often rank-ordering) the most important parameters for the purposes of, for example, deciding where to devote most calibration efforts in the model or most further analysis to reduce the uncertainty in a parameter. The figures below show the evolution of that rank-ordering as we increase *n*.

These figures show that with a number of samples that is too small, we could potentially misclassify a factor as important or unimportant when it actually is not.

Now, one might ask: *how is this useful? I’m trying to minimize my n, but you’ve just made me run way too many model evaluations, multiple times, just to determine how I could have done it faster? Isn’t that backwards?*

Well, yes and no. I’ve devoted this time to run this bigger experiment to get insight on the behavior of my model. I have established confidence in my index values and factor prioritization. Further, I now know that *n*>1500 would probably be unnecessary for this system and even if the model itself or my parameter ranges change. **As long as the parameter interactions, and model complexity remain relatively the same**, I can leverage this information to perform future sensitivity analyses, with a known minimum *n* needed.

References:

[1] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola, “Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index,” *Computer Physics Communications*, vol. 181, no. 2, pp. 259–270, Feb. 2010, doi: 10.1016/j.cpc.2009.09.018.

[2] F. Pianosi and T. Wagener, “A simple and efficient method for global sensitivity analysis based on cumulative distribution functions,” *Environmental Modelling & Software*, vol. 67, pp. 1–11, May 2015, doi: 10.1016/j.envsoft.2015.01.004.

[3] J. Nossent, P. Elsen, and W. Bauwens, “Sobol’ sensitivity analysis of a complex environmental model,” *Environmental Modelling & Software*, vol. 26, no. 12, pp. 1515–1525, Dec. 2011, doi: 10.1016/j.envsoft.2011.08.010.

Hi,

How do you plot these figures?

I am having problem when I am trying to plot them.

Could you show me the plots?

Thank you in advance.

Hi Alan, thanks for reading! The script to plot these figures can be found here: https://github.com/antonia-had/SA_index_convergence/blob/master/SAverification.py (lines 46 and 75).

Thank you so much, for the prompt reply.

The plots are working! 🙂

Hi, I have other questions.

1. Is it a good idea to compute the second order?

2. What can this show us?

3. If I try to plot the second order S2, how could I represent it?

Thank you again in advance.

The second order index gives you information about the interactions between the different parameters, specifically how much does a interact with b (for example). The total order does include the first order plus all interactive effects, but does not distinguish where the interactive effects are. So you might know that parameter a has a lot of interactions with other parameters but you won’t know which specifically. The second order index can tell you that.

Whether or not to compute it depends on what your questions about the model are and what you’re trying learn using a sensitivity analysis. For the purposes of illustration here it wasn’t important to compute.

You can plot S2 the same way as I plot the others here if you was to visualize how it changes with more n.

If you’d like to plot all three indices together, one common plot is the radial convergence or chord plot, I wrote a blogpost about them here: https://waterprogramming.wordpress.com/2019/05/06/radial-convergence-diagram-aka-chord-diagram-for-sensitivity-analysis-results-and-other-inter-relationships-between-data/

Another option is a Circos plot. Nice use of it can be found in this paper: https://doi.org/10.1002/wrcr.20413 (figures 7 and 8).

Dear Antonia,

Thank you very much for your help. I am still learning about sensitivity analysis. And your blog was essential to understand it.

If you have any further questions, I’ll be back here! 🙂

Kind regards!

Hi Antonia,

I have another question. I am testing the Sobol method for the SIR model. I followed in your steps, but I’m not sure if I’m doing it right.

Could you check if what I did is right? My question is mainly in the loop where S1 and ST are estimated.

Thank you again for your help.

Kind Regards.

#############

# Code

#############

#The ODE is:

def myModel(INP, t, beta0, gamma):

”’The main set of equations”’

Y=np.zeros((3))

V = INP

Y[0] = – beta0 * V[0] * V[1] #Susceptible

Y[1] = beta0 * V[0] * V[1] – gamma * V[1] #Infectious

Y[2] = gamma * V[1] #Removed

return Y # For odeint

# Initial condition

beta0 = 0.4

gamma = 1/6.

INPUT = (1, 1, 0.0)

# Solving ODE

t_range = np.arange(0, 30, 1)

solution = spi.odeint(myModel, INPUT, t_range, args=(beta0, gamma))

# Sensitivity analysis

problem = {

‘num_vars’: 2,

‘names’: [‘beta0’, ‘gamma’],

‘bounds’: np.column_stack((np.array([0.0, 0.0]),np.array([1, 1])))

}

# Array with n’s to use

nsamples = np.arange(50, 4050, 50)

# Arrays to store the index estimates

S1_estimates = np.zeros([problem[‘num_vars’],len(nsamples)])

ST_estimates = np.zeros([problem[‘num_vars’],len(nsamples)])

# Loop through all n values, create sample, evaluate model and estimate S1 & ST

for i in range(len(nsamples)):

print(‘n= ‘+ str(nsamples[i]))

# Generate samples

sampleset = saltelli.sample(problem, nsamples[i],calc_second_order=False)

# Run model for all samples

output = [sp.odeint(diff_eqs2,[1,1,0.0],t_range,\

args=(sampleset[j][0],sampleset[j][1])) for j in range(len(sampleset))]

sol_ = np.asarray(output)

# Perform analysis

#Susceptible

result= sobol.analyze(problem, sol_[:,i,0], calc_second_order=False, print_to_console=False)

# Store estimates

ST_estimates_S[:,i]=results[‘ST’]

S1_estimates_S[:,i]=results[‘S1’]

#Infectious

result= sobol.analyze(problem, sol_[:,i,1], calc_second_order=False, print_to_console=False)

# Store estimates

ST_estimates_I[:,i]=results[‘ST’]

S1_estimates_I[:,i]=results[‘S1’]

#Removed

result= sobol.analyze(problem, sol_[:,i,2], calc_second_order=False, print_to_console=False)

# Store estimates

ST_estimates_R[:,i]=results[‘ST’]

S1_estimates_R[:,i]=results[‘S1’]

The indentation that I included in the comment has disappeared, sorry about that.

If you have any difficulty reading this, I send you another comment.

Sorry, there is an error on step “# Run model for all samples”. The correct one is that:

# Run model for all samples

output = [sp.odeint(myModel, INPUT,t_range,\

args=(sampleset[j][0],sampleset[j][1])) for j in range(len(sampleset))]

sol_ = np.asarray(output)

Hi Alan, I have a hard time running this without the indentations and the imported packages. If you still need help, email me your script at ah986[at]cornell.edu and I can take a look.