Tutorials › Quantitative Systems Pharmacology (QSP) › Saltelli Sampling for Sobol Index Estimation

Sensitivity Analysis and Model Reduction · Tutorial 665 of 799

Saltelli Sampling for Sobol Index Estimation

Construct paired Saltelli samples and calculate first- and total-order Sobol indices with explicit estimators and model-evaluation counts.

Advanced 12 min read

What You'll Learn

  • Build independent base matrices and one hybrid matrix per input.
  • Identify which model outputs each Saltelli estimator pairs.
  • Calculate first-order indices with a Saltelli covariance estimator.
  • Calculate total-order indices with the Jansen squared-difference estimator.
  • Count evaluations for one index and for all first- and total-order indices.
  • Recognize assumptions and implementation checks for pick-freeze sampling.

Why the Sample Design Matters

A Sobol index describes how variation in a model output is apportioned across uncertain inputs under a specified input distribution. Estimating these indices by simulation requires more than a large collection of parameter draws: the estimator depends on deliberately paired model evaluations. Saltelli sampling provides that pairing through a set of pick-freeze matrices. Each hybrid row “freezes” most inputs at one draw and replaces a selected input with its value from another draw.

This tutorial focuses on constructing those matrices and calculating first- and total-order indices from their model outputs. The indices depend on the assumed input distributions, so the sample construction must represent those distributions as well as preserve the required independence among input columns.

Assumptions: The inputs \(X_1,\ldots,X_d\) are mutually independent random variables with specified distributions. The model is a scalar-valued function \(Y=f(X_1,\ldots,X_d)\) with finite, positive output variance. The same model definition and numerical settings are used for every evaluation. The standard Sobol decomposition and the estimators below are not generally valid for dependent inputs.

Let \(A\) and \(B\) be independent \(N \times d\) sample matrices. Each row is an independent draw from the joint input distribution, which under the independence assumption is the product of the marginal distributions. For each input \(i\), define \(AB^{(i)}\) to be the matrix that starts as \(A\), but takes column \(i\) from \(B\). In other words, row \(r\) of \(AB^{(i)}\) has the \(i\)th input from row \(r\) of \(B\) and all its other inputs from row \(r\) of \(A\).

$$ AB^{(i)}_{rj} = \begin{cases} B_{ri}, & j=i,\\ A_{rj}, & j\ne i. \end{cases} $$

The input coordinates in each row of \(A\), \(B\), and each hybrid matrix must have the intended marginal distributions. A common implementation starts with independent draws on the unit hypercube and transforms each column using the inverse cumulative distribution function for the corresponding input. The hybrid is then formed by swapping the already-transformed input column. This is also where an independence assumption matters: replacing one column in this way preserves the intended product distribution, but does not preserve an arbitrary joint distribution with correlated inputs.

Estimators and Their Pairings

Write \(f(A_r)\), \(f(B_r)\), and \(f(AB^{(i)}_r)\) for the model outputs from row \(r\) of each matrix. Let \(S_i\) denote the first-order index for input \(i\), and let \(S_{T_i}\) denote its total-order index. The first-order index measures the variance contribution associated with that input alone; the total-order index also includes contributions from interactions involving that input. For their variance-decomposition definitions, see the earlier material on first- and total-order Sobol indices.

Estimator choice: Use the Saltelli covariance estimator for \(S_i\) and the Jansen squared-difference estimator for \(S_{T_i}\). Both can use the same \(A\), \(B\), and \(AB^{(i)}\) evaluations. Estimate output variance from the pooled outputs of \(A\) and \(B\), unless a different variance estimate is explicitly specified.
$$ \widehat{S}_i = \frac{ \displaystyle \frac{1}{N}\sum_{r=1}^{N} f(B_r)\left[f(AB^{(i)}_r)-f(A_r)\right] }{ \widehat{V} } $$
$$ \widehat{S}_{T_i} = \frac{ \displaystyle \frac{1}{2N}\sum_{r=1}^{N} \left[f(A_r)-f(AB^{(i)}_r)\right]^2 }{ \widehat{V} } $$

For the pooled sample-variance choice used here, let \(\bar{f}_{AB}\) be the mean of the \(2N\) outputs from \(A\) and \(B\). Then

$$ \widehat{V} = \frac{1}{2N-1} \sum_{r=1}^{N} \left[ \left(f(A_r)-\bar{f}_{AB}\right)^2+ \left(f(B_r)-\bar{f}_{AB}\right)^2 \right]. $$

The first-order numerator pairs \(f(B_r)\) with the hybrid output and subtracts the output from \(A_r\). The total-order numerator instead compares the two outputs that share all inputs except \(X_i\). It is important not to interchange these pairings: the same matrices support both estimators, but the formulas use them differently. These are finite-sample estimates, so an estimated index can fall slightly below zero or above one; that does not change the theoretical range of the corresponding index.

Why the Pairings Estimate the Indices

The core justification is that each pair shares a specified set of inputs. For the first-order estimator, \(B_r\) and \(AB^{(i)}_r\) share \(X_i\), while \(A_r\) is independent of \(B_r\). Conditional on \(X_i\), the two shared-input outputs are independent draws of the model’s remaining inputs. Therefore,

$$ \mathbb{E}\!\left[f(B)f(AB^{(i)})\right] = \mathbb{E}\!\left[\mathbb{E}(Y\mid X_i)^2\right], \qquad \mathbb{E}\!\left[f(B)f(A)\right] = \mathbb{E}(Y)^2. $$

Subtracting the second expectation from the first gives \(\operatorname{Var}(\mathbb{E}[Y\mid X_i])\), the first-order variance contribution. Dividing by \(\operatorname{Var}(Y)\) gives \(S_i\). This explains the product and subtraction in the Saltelli numerator; the expectation of its numerator is the required variance contribution.

For the total-order estimator, \(A\) and \(AB^{(i)}\) share all inputs except \(X_i\). Conditional on the shared inputs \(X_{-i}\), their \(X_i\) values are independent draws from the same distribution. The conditional expected squared difference is therefore twice the conditional variance over \(X_i\):

$$ \mathbb{E}\!\left[ \left(f(A)-f(AB^{(i)})\right)^2 \mid X_{-i} \right] = 2\operatorname{Var}(Y\mid X_{-i}). $$

Taking expectations over \(X_{-i}\) and dividing by \(2\operatorname{Var}(Y)\) gives the total-order index. In practice, the unknown expectations and variance are replaced by the sample averages and variance estimate above. This justification depends on the shared-input structure and on independently sampled input columns; an arbitrary pairing of model outputs does not have the same interpretation.

Worked Matrix Construction and Evaluation Count

Worked Example: Building Hybrids for Two Inputs

Suppose a hypothetical QSP model has two independently sampled inputs, \(X_1\) and \(X_2\), and we use \(N=2\) rows. For clarity, the entries below are illustrative values on the unit scale; in an application they would be draws from the specified input distributions or transformed unit-hypercube draws.

MatrixRow 1: \(X_1\)Row 1: \(X_2\)Row 2: \(X_1\)Row 2: \(X_2\)
\(A\)0.200.800.600.30
\(B\)0.900.100.100.70
\(AB^{(1)}\)0.900.800.100.30
\(AB^{(2)}\)0.200.100.600.70

For example, row 1 of \(AB^{(1)}\) takes \(X_1=0.90\) from \(B\) and \(X_2=0.80\) from \(A\). Row 1 of \(AB^{(2)}\) does the reverse. The same construction applies independently to every row.

To estimate all first- and total-order indices, evaluate the model at the \(N\) rows of \(A\), the \(N\) rows of \(B\), and the \(N\) rows of each of the \(d\) hybrids. The total is

$$ N(d+2)=2(2+2)=8 $$

model evaluations. This count assumes each matrix row requires one model evaluation and that there is no reuse from another calculation. It counts full model runs, not internal solver steps. With \(N\) rows and \(d\) inputs, the corresponding count for all first- and total-order estimates is \(N(d+2)\).

Worked First-Order Calculation

Worked Example: Estimating a First-Order Index

Consider a hypothetical four-row calculation for one input \(X_1\). Suppose the model returns the following outputs. The values are chosen to make the arithmetic transparent; they do not represent a real study.

Row \(r\)\(f(A_r)\)\(f(B_r)\)\(f(AB^{(1)}_r)\)
10.0000\(-0.5000\)\(-0.05785\)
20.00000.50000.05785
30.0000\(-0.9798\)\(-0.1134\)
40.00000.97980.1134

Here the unrounded \(B\) outputs are \(-0.5, 0.5, -\sqrt{0.96}, \sqrt{0.96}\), and the unrounded hybrid outputs equal \(k f(B_r)\), where \(k=0.28/2.42\). Thus \(\sum_r f(B_r)^2=2.42\) and \(\sum_r f(B_r)f(AB^{(1)}_r)=0.28\). The displayed table rounds the values for readability; use the stated unrounded values for the calculation.

The \(A\) outputs are all zero, and the \(B\) outputs sum to zero, so the pooled mean is zero. The pooled sum of squared deviations is \(2.42\), and there are \(2N=8\) outputs. Therefore,

$$ \widehat{V}=\frac{2.42}{8-1}=\frac{2.42}{7}=0.3457. $$

The Saltelli numerator is

$$ \frac{1}{4}\sum_{r=1}^{4} f(B_r) \left[f(AB^{(1)}_r)-f(A_r)\right] = \frac{0.28}{4} = 0.0700. $$

Hence the estimated first-order index is

$$ \widehat{S}_1 = \frac{0.0700}{2.42/7} = 0.202479\ldots \approx 0.2025. $$

The estimate is dimensionless. Under the specified input distributions, it estimates the fraction of output variance attributable to the first-order contribution of \(X_1\); it is not a causal effect or a statement that the input explains exactly that fraction in every sample.

Worked Total-Order Calculation

Worked Example: Estimating a Total-Order Index

For another hypothetical calculation, suppose the outputs for one input \(X_i\) are as follows:

Row \(r\)\(f(A_r)\)\(f(B_r)\)\(f(AB^{(i)}_r)\)
1102
2232
3011
4321

The eight pooled \(A\) and \(B\) outputs have mean \(1.5\). Their squared deviations sum to \(10\), so the specified pooled sample variance is \(\widehat{V}=10/7=1.4286\). The four paired squared differences are \(1,0,1,4\), giving

$$ \widehat{S}_{T_i} = \frac{(1+0+1+4)/8}{10/7} = \frac{0.7500}{1.4286} = 0.5250. $$

The estimate includes the first-order contribution of \(X_i\) and any interaction contributions involving it. For \(d=3\) inputs and \(N=4\) rows, estimating all first- and total-order indices with this shared design requires \(4(3+2)=20\) model evaluations, even though this example displays the outputs for only one hybrid.

Evaluation Counts, Alternatives, and Checks

The main computational advantage of the shared Saltelli design is that the same base matrices support estimates for every input. The count depends on what is being estimated and which outputs are needed:

Target using the estimators aboveRequired model evaluations
First- and total-order indices for all \(d\) inputs\(N(d+2)\): \(A\), \(B\), and \(d\) hybrid matrices
One input’s first- and total-order indices\(3N\): \(A\), \(B\), and that input’s hybrid
First-order indices for all \(d\) inputs\(N(d+2)\) with this shared design; some alternative designs can reduce evaluations if total-order estimates are not required

A nested conditional Monte Carlo approach is a useful conceptual alternative: it approximates conditional expectations by repeatedly sampling the remaining inputs at fixed values of a selected input. It can be straightforward to explain, but nested sampling may require many evaluations and requires careful choices of inner and outer sample sizes. Pick-freeze estimation replaces that nested structure with paired matrices and simple averages. Other pick-freeze estimators rearrange which base or hybrid outputs appear in the numerator; their formulas and sample pairings must be used together rather than mixed.

Before interpreting results, verify that each hybrid really differs from \(A\) only in its designated input column, that \(A\) and \(B\) are independently generated, and that every column transformation matches the intended input distribution. Check that all evaluations use the same output definition, units, and solver settings. If the model output is nearly constant, the variance denominator is near zero and normalized indices can be unstable or undefined. If the model is stochastic, independent simulation noise can contaminate the paired differences; use a deliberate noise-handling strategy rather than treating noisy outputs as deterministic model values.

A common mistake is to form \(AB^{(i)}\) by replacing every column except \(i\), then apply the formulas above. That hybrid shares \(X_i\) with \(A\), not with \(B\), so it does not match the first-order covariance pairing defined here. The correction is to write down the intended shared coordinates first: for \(AB^{(i)}\), take column \(i\) from \(B\) and all other columns from \(A\), and then use the corresponding Saltelli and Jansen formulas consistently.

Finite-sample estimates need not lie in the theoretical interval from zero to one, and indices need not sum to one because total-order indices include interactions. These facts are reasons to check the calculations and sampling design, not reasons to silently truncate an estimate or renormalize the indices. The design and estimators describe variance contributions under the specified distributions; changing those distributions changes the sensitivity question itself.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For \(N=500\) and \(d=6\), how many model evaluations are needed for all first- and total-order estimates using the shared design?
  2. In \(AB^{(3)}\), which matrix supplies column 3, and which matrix supplies the other columns?
  3. Why does the first-order numerator use \(f(B_r)[f(AB^{(i)}_r)-f(A_r)]\) rather than a squared difference?
  4. A total-order estimate is 1.04 in a finite run. Does this prove the theoretical total-order index exceeds one?
  5. Why can the standard hybrid construction be inappropriate for correlated inputs?

Solutions. (1) The count is \(N(d+2)=500(6+2)=4{,}000\) evaluations. (2) Column 3 comes from \(B\); all other columns come from \(A\). (3) The paired product estimates the conditional-mean variance contribution for the shared input, whereas the squared-difference pairing targets the total-order contribution. (4) No. The finite-sample estimate can exceed one because of sampling variability or a noisy denominator; the theoretical index is bounded by one under the stated framework. (5) Replacing one coordinate while leaving the others from another draw generally fails to preserve the joint dependence structure, so the hybrid rows need not have the intended joint distribution.

Clinical Biostats