Tutorials › Physiologically-Based Pharmacokinetics (PBPK) › Sobol Indices for Global Sensitivity

Sensitivity Analysis and Parameter Reduction · Tutorial 730 of 799

Sobol Indices for Global Sensitivity

Learn how Sobol indices partition output variance into individual parameter effects and contributions that involve interactions.

Advanced 10 min read

What You'll Learn

  • Define first-order and total-effect Sobol indices for a scalar PBPK output.
  • Derive how the indices correspond to variance components under independent inputs.
  • Calculate Sobol indices for additive, interacting, and nonlinear hypothetical PBPK endpoints.
  • Interpret gaps between first-order and total-effect indices as evidence of interaction contributions.
  • Recognize assumptions and limitations that affect what Sobol indices can establish.
  • Distinguish global variance-based indices from local sensitivities and Morris screening.

Why Use Variance-Based Global Indices?

A local sensitivity describes how an output changes near a chosen parameter value. Morris screening summarizes changes across parameter perturbations. Both can help identify influential parameters, but neither directly answers a variance-allocation question: how much of the output variation over specified parameter uncertainty is associated with each parameter, including its interactions with others? Sobol indices address that question for a defined output and a defined joint distribution of inputs.

For a PBPK model, the output might be an AUC, a concentration at a specified time, or another scalar endpoint. The indices depend on both the model and the uncertainty distributions assigned to its parameters. They are therefore not intrinsic rankings of parameters independent of the analysis conditions.

Definition: Let \(Y=f(X_1,\ldots,X_d)\) be a scalar output with finite, nonzero variance, where the uncertain inputs \(X_1,\ldots,X_d\) are mutually independent and have specified probability distributions. The first-order Sobol index for input \(X_i\) is $$ S_i=\frac{\operatorname{Var}\!\left(\operatorname{E}[Y\mid X_i]\right)}{\operatorname{Var}(Y)}. $$ The total-effect Sobol index for \(X_i\) is $$ S_{T_i}=\frac{\operatorname{E}\!\left[\operatorname{Var}(Y\mid X_{-i})\right]}{\operatorname{Var}(Y)}, $$ where \(X_{-i}\) denotes all inputs except \(X_i\).

The first-order index measures the fraction of output variance explained by the conditional mean as \(X_i\) varies over its distribution, averaging over the other inputs. The total-effect index measures the fraction associated with every variance component that involves \(X_i\): its main effect and all interactions containing it. Neither index is a probability that a parameter is important, nor does either establish that changing a parameter causes a real-world change outside the model and input assumptions.

Variance Decomposition and Interpretation

Under independent inputs and a square-integrable model output, the Sobol (functional ANOVA) decomposition expresses the centered output as a sum of components: single-input effects, pairwise interactions, three-way interactions, and so on. These components are orthogonal under the input distribution, so their variances add. Write \(V=\operatorname{Var}(Y)\), and let \(V_u\) denote the variance component associated with a nonempty subset \(u\) of input indices. Then

$$ V=\sum_{\varnothing\ne u\subseteq\{1,\ldots,d\}} V_u,\qquad S_i=\frac{V_{\{i\}}}{V},\qquad S_{T_i}=\frac{\sum_{u:\,i\in u}V_u}{V}. $$

Thus \(S_i\) contains only the component unique to input \(i\), while \(S_{T_i}\) also includes pairwise and higher-order components involving that input. For independent inputs with this decomposition, \(0\le S_i\le S_{T_i}\le1\). The sum of first-order indices is at most one; it is less than one when variance is contributed by interactions. The sum of total-effect indices can exceed one because an interaction is counted in the total-effect index of every input involved in it.

Why the Definitions Match the Variance Components

The connection follows from conditional expectation and the orthogonality of the decomposition. Let \(f_i(X_i)\) be the single-input component, with mean zero. When conditioning on \(X_i\), terms involving other independent inputs average to zero, leaving

$$ \operatorname{E}[Y\mid X_i]=\operatorname{E}[Y]+f_i(X_i). $$

Taking the variance gives \(\operatorname{Var}(\operatorname{E}[Y\mid X_i])=\operatorname{Var}(f_i(X_i))=V_{\{i\}}\). This proves that the numerator of \(S_i\) is precisely the single-input variance component.

For the total effect, condition instead on all inputs except \(X_i\). Components that do not involve \(i\) are fixed by this conditioning, whereas components that do involve \(i\) retain variation when \(X_i\) varies. Consequently, \(\operatorname{Var}(Y\mid X_{-i})\) contains the conditional variation from those components, and averaging over \(X_{-i}\) sums their variance components:

$$ \operatorname{E}[\operatorname{Var}(Y\mid X_{-i})] =\sum_{u:\,i\in u}V_u. $$

Dividing by \(V\) gives \(S_{T_i}\). This derivation relies on independence: it allows the orthogonal variance decomposition and the clean separation of terms by input subset. It also assumes finite output variance. If \(V=0\), the indices are undefined because there is no output variance to apportion.

Worked Example: An Additive PBPK Endpoint

Suppose a hypothetical PBPK analysis considers an AUC endpoint \(Y\), measured in mg·h/L, with two independent, dimensionless uncertain inputs \(X_1\) and \(X_2\). Each is uniformly distributed on \([-1,1]\), and the model is

$$ Y=10+2X_1+X_2. $$

The input means are zero and their variances are \(1/3\). The output mean is 10 mg·h/L. Because this model is additive and the inputs are independent, the variance contributions are \(\operatorname{Var}(2X_1)=4/3\) and \(\operatorname{Var}(X_2)=1/3\), in squared output units. Therefore,

$$ \operatorname{Var}(Y)=\frac{4}{3}+\frac{1}{3}=\frac{5}{3}. $$

For \(X_1\), \(\operatorname{E}[Y\mid X_1]=10+2X_1\), so $$ S_1=\frac{4/3}{5/3}=0.800. $$ For \(X_2\), \(\operatorname{E}[Y\mid X_2]=10+X_2\), giving $$ S_2=\frac{1/3}{5/3}=0.200. $$ There are no interaction components, so \(S_{T_1}=S_1=0.800\) and \(S_{T_2}=S_2=0.200\). Over these specified input distributions, \(X_1\) accounts for 80.0% of the AUC variance and \(X_2\) for 20.0%. This statement concerns the model output’s variance under these assumptions; it is not a universal ranking for other parameter ranges or endpoints.

Main Effects and Interaction Contributions

A difference between a total-effect and first-order index indicates variance contributions involving that input and at least one other input. The difference is not attributable to one particular partner unless the interaction components are examined separately. In particular, \(S_{T_i}-S_i\) collects all interaction orders involving \(i\).

Worked Example: A Pairwise Interaction

Consider a hypothetical concentration endpoint \(Y\), measured in mg/L, with independent dimensionless inputs \(X_1,X_2\sim U(-1,1)\):

$$ Y=10+X_1+X_2+2X_1X_2. $$

Both inputs have mean zero and variance \(1/3\). The three nonconstant terms are mutually orthogonal: the product term has zero covariance with each single-input term because the inputs are independent and centered. Its variance is $$ \operatorname{Var}(2X_1X_2)=4\operatorname{E}[X_1^2]\operatorname{E}[X_2^2] =4\left(\frac{1}{3}\right)\left(\frac{1}{3}\right)=\frac{4}{9}. $$ The two main-effect variances are each \(1/3\), so the total output variance is $$ V=\frac{1}{3}+\frac{1}{3}+\frac{4}{9}=\frac{10}{9}. $$

Conditioning on \(X_1\) gives \(\operatorname{E}[Y\mid X_1]=10+X_1\), because \(\operatorname{E}[X_2]=0\). Hence \(S_1=(1/3)/(10/9)=0.300\). By symmetry, \(S_2=0.300\). For \(X_1\), the components involving it are its main effect and the interaction, so $$ S_{T_1}=\frac{1/3+4/9}{10/9}=\frac{7}{10}=0.700. $$ Likewise, \(S_{T_2}=0.700\). The 0.400 gap for each input reflects the pairwise interaction variance. The first-order indices sum to 0.600, leaving 0.400 of total variance in the interaction. The total-effect indices sum to 1.400 because the same interaction is included in both indices.

Worked Example: A Nonlinear Main Effect

Suppose a hypothetical PBPK endpoint \(Y\), in mg/L, follows $$ Y=5+X_1^2+X_2, $$ where \(X_1\) and \(X_2\) are independent and uniformly distributed on \([-1,1]\). This model has a nonlinear effect of \(X_1\), but no interaction. For a uniform variable on this interval, \(\operatorname{E}[X_1^2]=1/3\) and \(\operatorname{E}[X_1^4]=1/5\). Thus, $$ \operatorname{Var}(X_1^2)=\frac{1}{5}-\left(\frac{1}{3}\right)^2=\frac{4}{45}, \qquad \operatorname{Var}(X_2)=\frac{1}{3}=\frac{15}{45}. $$ The total variance is \(19/45\) in squared output units. The conditional means are \(\operatorname{E}[Y\mid X_1]=5+X_1^2\) and \(\operatorname{E}[Y\mid X_2]=5+1/3+X_2\). Therefore, $$ S_1=S_{T_1}=\frac{4/45}{19/45}=\frac{4}{19}\approx0.211, \qquad S_2=S_{T_2}=\frac{15}{19}\approx0.789. $$ The indices for \(X_1\) are nonzero even though the derivative of \(X_1^2\) with respect to \(X_1\) is zero at \(X_1=0\). The global index reflects variation over the full specified input range, not only behavior at one nominal value. Since the model is additive in its transformed effects, each total-effect index equals its first-order index.

Assumptions, Diagnostics, and Alternatives

Sobol indices answer a well-defined question only when the output, inputs, and uncertainty distributions are specified. The standard interpretation above assumes mutually independent inputs, a scalar output with finite nonzero variance, and a model that can be evaluated over the input space. A deterministic PBPK model naturally fits the output definition. If an output also contains stochastic simulation noise, that noise can contribute to the variance denominator and alter the indices; it should not be mistaken for parameter-driven variation.

Check that the selected input ranges or distributions represent the uncertainty question of interest, that the output is the intended endpoint and time point, and that the model produces finite values throughout the relevant input space. A very small output variance makes ratios unstable or practically uninformative. With numerical estimates, sampling error or model-evaluation error can also produce estimates outside theoretical bounds or inconsistent rankings; such results warrant investigation rather than literal interpretation. These checks do not replace the more detailed design and convergence checks needed for a Sobol simulation.

The independence assumption is especially consequential. With correlated parameters, the usual independent-input variance components do not have the same clean interpretation, because information about one input can be carried by another. The indices should not be applied mechanically to correlated PBPK parameters as if the standard decomposition still held. Likewise, changing parameter distributions can change the indices even when the PBPK equations are unchanged.

ApproachWhat it summarizesImportant distinction
Local sensitivityOutput response near a specified parameter valueCan miss behavior elsewhere in the input range
Morris screeningChanges from selected input perturbations across the domainUseful for screening, but does not directly allocate output variance
First-order Sobol indexVariance due to an input’s main effectDoes not include that input’s interaction contributions
Total-effect Sobol indexVariance due to the input’s main effect and all interactions containing itDoes not identify which other input forms each interaction

Common Pitfall

Wrong reasoning: “The first-order indices sum to 0.6, so the remaining 0.4 of variance must be explained by measurement error.” Correction: Under the independent-input variance decomposition, the remaining variance can be due to interactions, as in the pairwise example. The sum of first-order indices alone does not distinguish interactions from other sources omitted from the model or analysis. Also, a total-effect index is not an additional variance share that can be added to the first-order indices: it already includes the input’s main effect.

Key takeaway: First-order indices quantify individual main effects; total-effect indices include every interaction involving the input. Interpret both relative to the specified PBPK output, model, and independent input distributions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

  1. For independent inputs, what does a positive difference \(S_{T_i}-S_i\) indicate?
  2. If first-order indices sum to 0.75, can you conclude that the remaining 0.25 is measurement error? Explain.
  3. A parameter has \(S_i=0.10\) and \(S_{T_i}=0.65\). Interpret both values in context of output variance.
  4. What happens to the standard interpretation of these indices if the input parameters are correlated?

Solutions. (1) It indicates variance from one or more interaction components involving \(X_i\); it does not identify the interaction partner by itself. (2) No. Under the independent-input decomposition, interactions can account for variance not represented by first-order terms; other omitted variation may also matter. (3) The parameter’s main effect accounts for 10% of output variance, while its main effect plus all interactions involving it account for 65%. The difference, 55 percentage points, is interaction contribution involving that parameter. (4) The standard independent-input decomposition no longer applies in the same straightforward way, so the usual indices cannot be interpreted mechanically as clean independent variance shares.