Designing and Checking Sobol Simulations

From Sobol Definitions to a Reproducible Simulation

Sobol indices describe how uncertainty in model inputs contributes to the variance of an output under specified input distributions. Their definitions are population quantities; a PBPK analysis obtains estimates by evaluating the model at a finite collection of sampled parameter sets. Designing that collection well is essential: a formula cannot recover reliable indices from a biased, dependent, or inadequately sized sample.

Consider a deterministic PBPK model with \(d\) uncertain parameters and a scalar output \(Y=f(X_1,\ldots,X_d)\), such as the predicted area under the concentration–time curve in a stated dose scenario. The simulation design in this tutorial assumes the input parameters are mutually independent, their probability distributions have been specified, and the output has finite, nonzero variance. The model, dose, scenario, endpoint, units, and solver settings are held fixed throughout the design. The indices therefore describe this output under these distributions and conditions; they are not universal properties of the model.

To construct a row, first generate independent draws from convenient base distributions, often uniform draws on \([0,1]\), and transform each column using the inverse cumulative distribution function of its specified parameter distribution. For example, a positive parameter may use a lognormal distribution if that is the declared uncertainty model. The transformation must preserve the specified marginal distribution. Generating independent columns is appropriate only under the independence assumption; it does not establish that independence is scientifically justified.

Constructing the Matrices

The two base matrices must be sampled independently from one another, and rows within each matrix must follow the intended joint distribution. For ordinary Monte Carlo, rows are independent random draws. A space-filling or randomized low-discrepancy design can be used instead, but its randomization and analysis must be specified; a low-discrepancy sequence is not a substitute for defining the input distributions. In either case, the same row index must be kept together across \(A\), \(B\), and every hybrid when forming estimators.

In general, the base design with \(A\), \(B\), and one hybrid per parameter requires \(N(d+2)\) model evaluations. If the model is expensive, this cost should be planned before running it. Reusing a single matrix or changing the row pairing after the runs breaks the intended estimator structure. Record the random seed or design-generation settings, parameter transformations, row order, model version, solver tolerances, and output extraction procedure.

Estimating First-Order and Total-Effect Indices

Let \(Y_A\), \(Y_B\), and \(Y_{A_B^{(i)}}\) denote the vectors of \(N\) model outputs from the corresponding matrices. The first-order index \(S_i\) measures the fraction of output variance attributable to parameter \(i\) alone, while the total-effect index \(S_{T_i}\) includes its contribution through interactions with other inputs. The following estimators use the matrix convention just defined. They are finite-sample estimates, not exact identities for a particular finite design.

Here \(\widehat{V}\) is an estimate of the output variance, obtained consistently for all parameters from the simulated outputs. One simple choice is the empirical variance of the \(2N\) outputs pooled from \(A\) and \(B\), with divisor \(2N\) when treating the simulated collection as the variance estimate. Other established estimator conventions use different variance normalizations or index formulas; do not combine a numerator from one convention with a denominator from another. The convention and variance estimate should be reported.

The first-order numerator compares the model output under \(B\) with the output in which only parameter \(i\) is switched to its \(B\) value while the other parameters come from \(A\). Under independent inputs, its expectation is the variance of the conditional mean of \(Y\) given \(X_i\). The total-effect numerator instead measures the squared output difference after all inputs except \(i\) are held at the matched \(B\) values. Its expectation, after division by \(2\), is the variance contribution not explained by holding \(X_i\) fixed. These identities explain why input independence and matched rows matter.

In the population, the usual independent-input Sobol indices lie between zero and one, and the total-effect index is at least as large as the first-order index. Finite-sample estimates need not obey either constraint. A negative first-order estimate or a total-effect estimate above one can occur from Monte Carlo error, especially at small \(N\) or when output variance is poorly estimated. Do not silently truncate such estimates to zero or one: truncation hides a diagnostic signal and changes the estimator.

Checking Monte Carlo Error and Convergence

A single run at a single sample size does not demonstrate convergence. For ordinary Monte Carlo, a practical check is to calculate estimates at increasing \(N\) and compare them with estimates from independent designs. A batch or replicate approach gives a direct view of between-design variation: generate multiple independent designs, calculate the indices for each, and inspect the spread and stability of the estimates. This spread is a diagnostic for simulation error, not a guarantee that the design has explored every relevant model feature.

Useful diagnostics include the index estimate versus \(N\), replicate-to-replicate spread, the output variance estimate, and the first-order and total-effect estimates together. Investigate estimates that drift substantially with increasing \(N\), fluctuate widely across independent designs, or repeatedly violate expected population relationships. If a parameter’s total-effect estimate is close to zero, its squared-difference numerator can be particularly sensitive to numerical noise. A flat convergence plot can also be misleading if every design shares the same coding error, wrong distribution, or inadequate parameter range.

For a deterministic PBPK model, repeat evaluations at selected parameter sets when numerical solver error may be consequential. If the same inputs produce materially different outputs, resolve solver tolerances, event handling, or stochastic simulation noise before interpreting design-to-design variability as Monte Carlo error. Failed or non-finite model evaluations must be investigated. Dropping those rows without a prespecified treatment can change the sampled input distribution and bias the indices; report the failure pattern and address its cause or justify an explicit model-domain rule.

Alternatives, Limitations, and Common Mistakes

A direct variance-decomposition calculation is conceptually simple but usually needs many model evaluations as the number of parameters grows. The A–B hybrid design reuses runs to estimate effects efficiently and supports first-order and total-effect estimates from a common sample. A randomized low-discrepancy design may improve coverage for some problems, while ordinary independent Monte Carlo has a straightforward replicate-based error assessment. Neither design is universally superior: assess actual stability, preserve the required matrix structure, and document the method. The phrase “Sobol sequence” refers to one kind of low-discrepancy sampling sequence; it is distinct from the Sobol sensitivity indices.

These estimators rely on a correctly specified joint input distribution, in particular the independence assumption used in the derivation. They also require finite, nonzero output variance. If the output is constant over the sampled domain, the denominator is zero and variance-based indices are undefined. Results can be misleading if input ranges or probability distributions do not represent the uncertainty of interest, if the model has unhandled failures, or if the sample is too small for a rare but influential region. A scalar endpoint’s indices also need not describe another endpoint or a different dose scenario.

• Mistake: treating one set of indices as converged because the calculation completed. Completion means the model ran; it does not establish adequate precision. Compare increasing sample sizes and independent designs.
• Mistake: interpreting a finite-sample estimate outside zero to one as a true negative contribution. Population indices under the stated assumptions are nonnegative and bounded, but noisy estimates can fall outside those limits. Report the diagnostic and assess stability rather than silently clipping.
• Mistake: permuting hybrid rows or using correlated columns while applying independent-input formulas. These changes invalidate the estimator’s pairing or assumptions. Preserve row alignment and use a method appropriate to the actual joint distribution.
• Mistake: treating failed runs as harmless missing values. Failure can depend on parameter values and therefore alter the sampled domain. Diagnose failures and document any justified treatment.