5.3. Sensitivity Analysis

Sensitivity analysis functionality of GTSDA allows to perform sensitivity analysis on user-provided data. Sensitivity analysis is defined as the study of how the variation (uncertainty) in the output of a statistical model can be attributed to different variations in the inputs of the model. In other words it is a technique for systematically changing variables (features) in a model to determine the effects of such changes.

5.3.1. Introduction

This functionality of GTSDA allows computing sensitivity indices (scores) to rank inputs according to selected criteria that measures effect of each input \(x\) on output \(y\). The user-provided dependency (also known as function or model) can be represented as data sample or interface to some blackbox (some device, system or object that provides output for a given input).

For a sensitivity analysis technique we wish it to share the following properties:

• If one variable is more important than the other in a technique defined way, we want its score to be higher.
• We want feature scores to be proportional to the corresponding variables influence, so that comparing scores one would get the idea of relative importance of variables.

These properties allow to rank features in the order of importance and give the idea of approximately to what extent one feature is more important than other features. Exact meaning of importance depends on the exact method selected. Based on computed scores, one may understand if dependency is linear, monotonic, additive, understand the strength of feature interactions.

Example of usage:

# in case sample is given
rankResult = gtsda.Analyzer().rank(x=x, y=y, options={}, approx_options={})

# in case blackbox is given
rankResult = gtsda.Analyzer().rank(blackbox=blackbox, budget=budget, options={})

Input parameters are:

• training sample or blackbox;
• GTSDA options (optional);
• GTApprox options (optional) - may be constructed for some techniques as a substitute to blackbox in case sample input is given.

The output parameters are:

• sensitivity scores;
• approximation model constructed (if any);
• generated sample (in case of blackbox mode).

The GTSDA options are described in detail in section Option Reference. Any combination of GTApprox options is supported.

5.3.1.1. Toy example

As an illustration we give the following simple example. Consider the Newton’s law of universal gravitation. Say we know that every point mass attract every other point mass, but don’t know what features affect that and what the form of dependence is.

For example, we selected following features that we expect may affect the force of attraction:

• \(m_1\), \(m_2\) - the masses of the bodies;
• \(r\) - distance between bodies;
• \(T\) - environment temperature;
• \(p\) - atmospheric pressure;
• \(L_1\), \(L_2\) - bodies luminosity.

The Newton’s law of universal gravitation

Imagine, we did \(30\) experiments where we varied input features and measured resulting force of attraction. After that applied GTSDA for the sample obtained (also if we had computational code that calculates attraction force we may have used it directly). The following results would follow:

• All sensitivity analysis techniques of GTSDA would show that attraction force only depends on masses and distance between bodies (distance being more important than masses). Masses have equal importance, while distance has greater effect on output.
• Additionally, computing Sobol Indices would tell us that interaction between features are strong, so dependency is not additive.
• Computing Screening Indices would show that attraction force monotonically depends on distance between bodies.

5.3.1.2. Input definition domain importance

It is important to note that the sensitivity scores returned by GTSDA depend on the variation intervals of the factors. If a factor is restricted to a very narrow interval, then its score might be low even if factor is important. However, the sensitivity scores returned by GTSDA are invariant under changes of units of measurement for individual factors (as long as changes are linear). In such cases the effects of rescaled intervals are compensated by the corresponding changes in the response function.

For example, consider the case when we have a function \(f(x_1,x_2) = x_1 + x_2\), with \(x_1 \in [-1, 1]\) and \(x_2 \in [-1, 1]\). It’s obvious to expect \(x_1\) and \(x_2\) to have equal scores in these conditions. Now, let us expand \(x_1\) to region \([-2, \, 2]\), while keeping \(f(x_1,x_2)\) the same. In this case though in each point local importance of \(x_1\) and \(x_2\) remains similar, on the global scale \(x_2\) would provide \(4\) times more variation to the output, thus rising its feature score. It’s equivalent to the case when we leave \(x_1\) at \([-1, \, 1]\) and change function to \(f(x_1,x_2) = 2*x_1 + x_2\). On the contrary, consider the case when we change the measurement units of the feature. For example, \(x_1\) and \(x_2\) were defined in kilograms and we want to change the measurement units of \(x_1\) to grams. In this case, though new values of rescaled \(x_1\) would become \(1000\) times larger but its feature score would remain the same.

5.3.2. Techniques

There are lots of approaches to the problem of global sensitivity analysis [1], [2], [3], [4], [5], [6], [7], [8], [9], [10]. Technique appropriate for each task depends on the problem conditions and user requirements. We’ve designed the GTSDA sensitivity analysis tool to include the most effective and popular state of the art methods, covering different problem settings. In this section brief overview of the sensitivity analysis techniques used in the GTSDA is provided.

At the moment three types of score indices are implemented:

• Screening Indices - its value may be considered as crude estimation of average partial derivative for each input. In addition to providing feature ranking index allow to check each input feature dependency for monotonicity and linearity.
• Sobol Indices - estimate the amount of variance of output described by each of input variables. In addition to providing feature ranking index allow to check if function is additive and allows to measure the strength of feature interactions.
• Taguchi Indices - estimates the importance of input variable for the either outputs signal-to-noise ratio, minimization or maximization. It requires specific design - Orthogonal Array (see DoE section).

To select index type one should set corresponding value to GTSDA/Ranker/Technique option. By default Screening Indices are computed.

5.3.2.1. Screening Indices

Changed in version 6.42: supports discrete, categorical, and stepped variables (previously supported continuous variables only).

This method is a screening technique able to work with relatively small samples. For linear functions scores returned by method coincide with coefficients for linear regression. For non-linear functions they may be sought as crude estimates of partial derivatives.

Screening indices are calculated using the Morris Screening algorithm. The idea of Morris Screening approach is to generate uniform (in terms of space-filling properties) set of trajectories or points in the design space. In the first case Morris Screening generate random trajectories, on each step of trajectory only one component \(x_i\) of input vector \(X\) is changed, so we get function values in different parts of domain. In the second case we get more uniform distribution of points, because we use LHS method for its generating.

The following function (elementary effect) is estimated:

\[d_i(X) = \frac{Y(x_1, \ldots, x_i+\delta_i, \ldots, x_p) - Y(x_1, \ldots, x_i, \ldots, x_p)}{\delta_i},\]

where \(\delta_i\) is a step size. If GTSDA/Ranker/NormalizeInputs is equal to:

• False, then \(\delta_i\) is absolute step size;
• True, then \(\delta_i\) is normalized to corresponding input range.

Score for \(i\)-th feature is computed as

\[\mu^*_i = \frac{1}{r}\sum_{j=1}^r \Big|d_i(X^j)\Big|, \; i=1, \ldots , p,\]

\[\mu_i = \frac{1}{r}\sum_{j=1}^r d_i(X^j), \; i=1, \ldots , p,\]

\[\sigma^2_i = Var\Big(d_i(X^j)\Big), \; i=1, \ldots , p,\]

where \(r\) is a number of steps changing \(i\)-th feature value on all trajectories, \(X^j\) is the input value at these steps, \(var\Big(\cdot\Big)\) is the variance.

The meaning of these quantities is as follows:

• \(\mu^*_i\) measures the strength of dependency: zero if output does not depend on \(i\)-th feature, in case of linear dependency \(\mu^*_i\) is exactly the value of linear coefficient.
• \(\mu_i\) measures the average direction of dependency of output and \(i\)-th feature: if it is positive, then output generally increases if one increases the value of \(i\)-th feature, if it is negative, then output generally decreases if one increases the value of \(i\)-th feature. If \(\mu_i\) is equal to \(\mu^*_i\) the dependency on the feature is monotonic.
• \(\sigma_i\) measures the degree of nonlinearity of dependency of output and \(i\)-th feature: zero value means linear dependency, \(\sigma_i\) being of magnitude of \(\mu^*_i\) means that dependency is very non-linear).

If sample input is given, the method trains an internal GTApprox model and uses it as a blackbox to generate trajectories.

Pros:

• May work even with very small budgets.
• Very effective at finding features that have no influence.
• Can measure the level of nonlinearity and monotonicity of dependency.
• Supports continuous, discrete, categorical, and stepped variables.

Cons:

• Requires an input-output dependency model (a blackbox). When used in sample-based mode, trains an internal GTApprox model to use it as the blackbox, which may take significant time for large samples.
• Scores do not have a strict physical meaning (though \(d_i\) may be thought as crude estimation of continuous derivatives).

5.3.2.2. Sobol Indices

The idea of Sobol indices is to measure what portion of output variance is described by the variance of the feature.

The score is the estimate of part of variation of output that is described by the variation of considered input.

Method can provide following types of indices:

• Main indices take into account only sole influence of feature while others being fixed. The index tells what portion of output variance would be described by considered input provided all other inputs are fixed at their mean values.
• Total indices take into account interactions of feature with others. The index tells what portion of output variance would be lost if we fix considered feature to its mean value, while still vary other inputs.
• Interaction indices are equal to difference between total and main indices, that represents strength of variables interactions.

For each feature main indices are estimated as

\[S_i = \frac{V_{x_i}[E_{\sim x_i}(Y|x_i)]}{V(Y)}\]

and total indices are estimated as

\[T_i = 1 - \frac{V_{\sim x_i}[E_{x_i}(Y|x_1, \ldots, x_{i-1}, x_{i+1}, \ldots, x_p)]}{V(Y)},\]

where \(E_{x_i}[\cdot]\) and \(V_{x_i}[\cdot]\) are the mean and the variance with respect to \(x_i\), whereas \(E_{\sim x_i}(\cdot|x_i)\) and \(V_{\sim x_i}(\cdot|x_i)\) are the conditional mean and variance with respect to all features except \(x_i\).

There are several methods to compute Sobol indices: FAST, CSTA and EASI. Each of the methods has its pros and cons (see below). In particular, FAST and CSTA methods are suited for computation of all the types of indices: total, main and interactions, while EASI is designed to compute specifically main indices. Additionally, with CSTA method, tool can compute second-order indices \(S_{i, j}\) which correspond to pairwise interactions between features.

5.3.2.2.1. Sobol Indices: FAST Method

Method uses space filling one-dimensional curves of the form

\[x_i(s)=\frac{1}{2}+\frac{1}{\pi}arcsin(sin(v_is+\phi _i))\]

to generate sample points. Here each feature has some frequency \(v_i\) assigned from some incommensurate set \({v_i}\), \(s\) is the coordinate on one-dimensional curve and \(\phi _i\) is a some random constant phase shift. Using Fourier decomposition we may say that

\[f(X) = \sum_{j=-\infty}^{\infty }(A_j \cos (js) + B_j\sin (js)),\]

\[A_j = \frac{1}{2\pi }\int_{-\pi}^{\pi}f(s)cos(js)ds,\]

\[B_j = \frac{1}{2\pi }\int_{-\pi}^{\pi}f(s)sin(js)ds.\]

These integrals can be estimated using points generated on the curve. In this case, e.g. conditional variance can be estimated as

\[V_{x_i}[E_{\sim x_i}(Y|x_i)] = 2 \sum_{j=1}^{K}(A_{jv_i}^2 + B_{jv_i}^2), \; \;jv_i \; \text{is an integer},\]

where \(K\) is some predefined number.

Another appealing property of this approach is its ability to accurately estimate total indices.

To do this estimation unique frequency \(v_i\) is given to \(x_i\) and the same frequency \(v\) is given to all other features, then the same procedure as above is performed.

Notably, both main and total indices can be estimated during one method run.

To generate sample method needs from given training sample FAST method constructs GT Approx model and uses it as a blackbox.

Method is implemented according to [2].

Pros:

• Can measure the strength of non-linear dependency.
• Provided it has enough budget, can measure precise values of main, interaction and total effects.

Cons:

• Requires a high budget to obtain meaningful results.
• Requires an input-output dependency model (a blackbox). When used in sample-based mode, trains an internal GTApprox model to use it as the blackbox, which may take significant time for large samples.
• Accuracy of estimates depends on the true importance of inputs: variables, which have greater effect, have more accurate estimates.
• Supports continuous variables only.

5.3.2.2.2. Sobol Indices: CSTA Method

Changed in version 6.2: now also provides second-order indices.

Changed in version 6.42: supports discrete, categorical, and stepped variables (previously supported continuous variables only).

The method is implemented according to [8]. The main idea of method is based on next theorems:

• For any input variable \(j\) , the main index \(S_j\) equals the expectation value of the correlation coefficient \(C_j\) of the output vectors from two model runs in which the realizations of all inputs in \(j\) are common to both runs, while all other inputs use independent samples in the two runs.
• For any input group \(j\) , the total effect \(T_j\) equals one minus the expectation value of the correlation coefficient \(C_{-j}\) of the output from two model runs which use independently sampled values for inputs in \(j\) , but use the same realizations in both runs for all other inputs.

This method also computes second-order indices:

• For any input variables \(i\), \(j\) , the main index \(S_{i,j}\) equals \(C_{\{i,j\}} - S_i - S_j\), where \(C_{\{i,j\}}\) is the expectation value of the correlation coefficient of the output vectors from two model runs in which the realizations of all inputs in \(i\), \(j\) are common to both runs, while all other inputs use independent samples in the two runs.
• For any input variables \(i\), \(j\) , the total effect \(T_{i, j}\) equals one minus the expectation value of the correlation coefficient \(C_{-\{i,j\}}\) of the output from two model runs which use independently sampled values for inputs in \(i\), \(j\) , but use the same realizations in both runs for all other inputs.

The method proceeds in the following steps:

• Randomly generate two sets of \(N\) sample values for each of \(J\) input variables. Here, these sets are called the primed and unprimed sets.
• Make \((2J + 2)N\) runs of the studied function.
• Standardize outputs.
• Make estimates for \(C_j\) and \(C_{-j}\).
• Make estimates for \(C_{\{i,j\}}\) and \(C_{-\{i,j\}}\).
• Make correction for spurious correlation between primed and unprimed sets.
• Construct confidence intervals for indices (see GTSDA/Ranker/VarianceEstimateRequired).

Pros:

• Can measure the strength of non-linear dependency.
• Provided it has enough budget, can measure precise values of main, interaction and total effects.
• Returns zero for main and total Sobol indices for variable \(X_i\), if \(Y\) does not depend on the variable \(X_i\).
• Accuracy of estimates does not depend on the true importance of inputs.
• Estimates second-order indices.
• Makes it easy to estimate confidence intervals.
• Supports continuous, discrete, categorical, and stepped variables.

Cons:

• Requires a high budget to obtain meaningful results.
• Requires an input-output dependency model (a blackbox). When used in sample-based mode, trains an internal GTApprox model to use it as the blackbox, which may take significant time for large samples.

5.3.2.2.3. Sobol Indices: EASI Method

EASI is a computationally cheap method for which existing data can be used.

Unlike the Sobol Indices: FAST Method method which use a specially generated sample set that contains suitable frequency data for the input factors. In EASI these frequencies are introduced by sorting and shuffling the available input samples.

The sorted input will be a nearly symmetric, periodic signal of frequency 1 (irrespectively of the input distribution).

The output data are sorted accordingly matching the resorted input samples, hence avoiding re-evaluation of the model. These sorted data can then be analysed using the power spectrum of the output. This latter analysis forms the standard back-end procedure of the Fourier-based techniques.

Method is implemented according to [7].

Pros:

• Works very fast even for large samples.
• Computes Sobol indices directly from the sample, with no need to train an internal approximation model, thus avoids introducing additional inaccuracy due to modeling.

Cons:

• Estimates only main indices.
• Does not support blackbox-based mode.
• Supports continuous variables only.

5.3.2.2.4. Noise Correction

New in version 6.6.

GT SDA allows to correct the values of Sobol indices in case of noisy output values. The correction can be switched on by option GTSDA/Ranker/NoiseCorrection and works only for Sobol Indices: FAST Method and Sobol Indices: CSTA Method methods.

The correction is done by adding a dummy variable that does not influence the output function. The estimates of main \(S_N\) and total \(T_N\) Sobol indices for this dummy variable allow to correct the values of Sobol indices of input variables in the following way:

\[ \begin{align}\begin{aligned}S_i^* = \frac{S_i}{1 - S_N},\\T_i^* = \frac{T_i - T_N}{1 - T_N},\end{aligned}\end{align} \]

where \(S_i\) and \(T_i\) are the values of main and total Sobol indices computed from noisy observations. It should be noted, that this particular type of correction implies the additive noise in output values.

5.3.2.3. Taguchi Indices

New in version 6.1.

Changed in version 6.42: supports discrete, categorical, and stepped variables (previously supported continuous variables only).

This method is an implementation of Taguchi techniques for engineering. It is constructed to work with noisy outputs and is normally used with small sample size. Taguchi indices require the design to be formed by discrete variables. The method shows the best quality if design forms so-called orthogonal array (see section Orthogonal Array for details).

In case of sample-based analysis, this technique checks whether the sample is a nearly orthogonal array. If it is not, results may be flawed, and in such cases GTSDA issues a warning. In case of blackbox analysis, an orthogonal array sample is generated by GTSDA (and you can specify the number of levels for continuous variables — see GTSDA/Ranker/Taguchi/LevelsNumber). In the both cases, you should define the analysis goal with the GTSDA/Ranker/Taguchi/Method option:

• maximum - search for the most significant variable for maximization of the mean output value;
• minimum - search for the most significant variable for minimization of the mean output value;
• signal-to-noise ratio - search for the most significant variable for maximization of the mean output value divided by its standard deviation.

Let \(\{x_1, \dots, x_N\}\) be the set of unique input points in sample \(X\). Our tool performs analysis for each output component separately. That’s why in what follows we will consider just one output \(y\). For Taguchi analysis it is important that at each unique point several evaluations of output are made. Let \(\{y_{i, 1}, \dots, y_{i, N_i}\}\) be the set of \(N_i\) output evaluations for unique input point \(x_i\).

The results of Taguchi analysis are scale invariant, but are not shift invariant. The canonical Taguchi analysis requires all output values to be greater than \(0\). We recommend to perform Taguchi analysis for such type of output.

The analysis starts with aggregating values for each unique point \(i\):

• maximum

  \[S_i = -10 \cdot \log_{10}\left(\frac{1}{N_i} \sum_{k=1}^{N_i} \frac{1}{y_{i, k}^2}\right);\]

• minimum

  \[S_i = -10 \cdot \log_{10}\left(\frac{1}{N_i} \sum_{k=1}^{N_i} y_{i, k}^2\right);\]

• signal-to-noise

  \[\overline{y}_i = \frac{1}{N_i} \sum_{k=1}^{N_i} y_{i, k};\]
  \[\sigma_i^2 = \frac{1}{N_i - 1} \sum_{k=1}^{N_i} \left(y_{i, k} - \overline{y}_i\right)^2;\]
  \[S_i = 10 \cdot \log_{10}\left(\frac{\overline{y}_i^2}{\sigma_i^2}\right).\]

It is important that if for some \(i\) all \(y_{i, k}\) for different \(k\) are equal, then the signal-to-noise statistic \(S_i\) is not defined. In this case the \(i\)-th unique point is ignored when calculating score (see below). Also if for some \(i\) all \(y_{i, k} = 0\), then the maximum and minimum statistics will have NaN values.

As a next step for each feature \(j=1, \dots, d\) and all its values from \(\{X_{j, 1}, \dots, X_{j, s_j}\}\), where \(s_j\) is the number of levels for \(j\)-th input (see section Orthogonal Array for details). For each possible value \(X_{j, k}\) of \(j\)-th input variable the following quantity is calculated:

\[S_{k, j} = \frac{s_j}{N} \sum_{i \colon x_{i, j}=X_{j, k}} S_i,\]

where \(X_{i, j}\) is the value of \(j\)-th component of \(X_i\).

Finally, the scores are calculated for each feature:

\[\max_{k \in \overline{1, s_j}} S_{k, j} - \min_{k \in \overline{1, s_j}} S_{k, j}.\]

For better understanding of Taguchi analysis please see Using Taguchi Analysis.

Pros:

• Handles noisy outputs.
• Can work with a low budget provided it is sufficient to generate an orthogonal array input sample.
• Supports continuous, discrete, categorical, and stepped variables.

Cons:

• In sample-based mode, requires a sample of a nearly orthogonal array design, otherwise results may be flawed.
• In blackbox-based mode, requires generating an orthogonal array design. That generation itself may take noticeable time depending on the array size (budget) and properties of variables.

5.3.3. Sample size and budget requirements

The maximum size of the training sample, which can be processed by GTSDA, is primarily determined by the user’s hardware. Necessary hardware resources depend significantly on the specific technique — see descriptions of individual techniques. Accuracy of estimation tends to improve as the sample size increases. But the maximum training sample size for the Taguchi Indices is limited by \(150\) unique input points and the problems input dimensionality is upper bounded by \(50\).

Contrary to the maximum size, there is a certain minimum for the size of the training set (or for the available number of blackbox calls), which depends on the technique used. This condition refers to the effective values, i.e. the ones obtained after preprocessing. An error with the corresponding error code will be returned if this condition is violated.

The requirements on minimum sample size (budget) are summarized in Table Minimum sample size (blackbox budget) for GTSDA sensitivity analysis techniques. It should be noted that all the techniques except Taguchi preprocess sample by removing of duplicate points. For most techniques there are two different limits, depending on whether we use blackbox or sample mode.

Table Minimum sample size (blackbox budget) for GTSDA sensitivity analysis techniques denotes the following:

• \(p\): the effective input dimension after the sample preprocessing.
• \(NR\): the GTSDA/Ranker/Sobol/FAST/NumberCurves option value. Corresponding limit is in effect only if the option is set by user.
• \(\left\lceil x \right\rceil\) is the value of \(x\) rounded up (to the next integer).

Minimum sample size (blackbox budget) for GTSDA sensitivity analysis techniques

Technique	minimum sample*	minimum blackbox budget
MS	\(2p+3\)	\(p+1\)
Sobol (FAST)	\(2p+3\)	\(17p\) or \(17p \cdot NR\)
Sobol (CSTA)	\(2p+3\)	\(4p+4\)
Sobol (EASI)	\(150\)
Taguchi	not specified, depends on proper Orthogonal Array existence

* Duplicate rows are removed at preprocessing stage (except Taguchi technique), so the requirement is on unique rows.

Note that even if sample size provided is enough tool still may throw GTException in case the dependency is too complex to estimate indices on given sample.

5.3.4. References

[1]	Saltelli. Global Sensitivity Analysis The Primer. Wiley, 2008.

[2]	(1, 2) Saltelli. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics, 41:39-56, 1999.

[3]	Hastie, R. Tibshirani, and J. Friedman. The elements of statistical learning: data mining, inference, and prediction. Springer, 2008.

[4]	Kononenko. An adaptation of relief for attribute estimation in regression. 1997.

[5]	Schulkopf. Measuring statistical dependence with hilbert-schmidt norms. In Algorithmic Learning Theory, 3734:63-74, 2005.

[6]	Kraskov. Estimating mutual information. Physical review. E, Statistical, nonlinear, and soft matter physics, 69:40-79, 2004.

[7]	(1, 2) Plischke. An Effective Algorithm for Computing Global Sensitivity Indices (EASI), 2010.

[8]	(1, 2) Pukelsheim. Optimal Design of Experiments. SIAM, 2006.

[9]	Atkinson, A. N. Donev, R. D. Tobias. Optimum Experimental Designs, with SAS. Oxford University Press, 2007.

[10]	Fraley, M. Oom, B. Terrien, J. Zalewski. Design of experiments via taguchi methods: orthogonal arrays. controls.engin.umich.edu, 2007