Diagnosing Monte Carlo Error in Sobol Estimates

Why Diagnose Monte Carlo Error?

A Sobol index is a property of a model output and a specified distribution of its inputs. In practice, it is usually estimated from a finite set of model evaluations, so the reported value is also affected by Monte Carlo error: the variation that would arise if the sampling design were drawn again at the same size. A point estimate alone does not show whether that variation is small enough for the intended interpretation.

This matters when a sensitivity analysis is used to prioritize experiments, model refinement, or computational effort. If two estimated indices are close, a small change in the sampled design may reverse their ranking. Conversely, an estimate can be stable and still be small in a practically relevant sense. These are different questions: Is the estimate precise? and Is the effect large enough to matter?

The diagnostics here concern Monte Carlo error in Sobol estimates under a specified input distribution and model. They do not measure uncertainty about whether that distribution is scientifically appropriate, whether the model is structurally adequate, or whether the model parameters are known. Those sources of uncertainty require separate consideration.

Indices and the Sampling Structure

Definition: For independent inputs $X_1,\ldots,X_d$, model output $Y=f(X_1,\ldots,X_d)$, and $0<\operatorname{Var}(Y)<\infty$, the first-order Sobol index for input $X_i$ is $$ S_i=\frac{\operatorname{Var}\bigl(\operatorname{E}[Y\mid X_i]\bigr)}{\operatorname{Var}(Y)}. $$ The total-effect index is $$ S_{T_i}=\frac{\operatorname{E}\bigl[\operatorname{Var}(Y\mid X_{-i})\bigr]}{\operatorname{Var}(Y)}, $$ where $X_{-i}$ denotes all inputs except $X_i$. The first-order index measures the contribution attributable to $X_i$ alone; the total-effect index includes contributions involving $X_i$ in interactions.

The definitions assume independent inputs, as in the standard Sobol variance decomposition. For an estimator based on Saltelli sampling, each base row supplies an input vector from each of two independently sampled matrices, conventionally called $A$ and $B$, together with hybrid vectors that replace one input at a time. All evaluations associated with one base row form a row block. The row blocks are the independent Monte Carlo sampling units; the evaluations within a block are deliberately related.

One common estimator pair uses the first-order estimate $ \widehat S_i=\widehat{\operatorname{E}[Y_B(Y_{A_B^{(i)}}-Y_A)]}/\widehat{\operatorname{Var}(Y)} $ and the Jansen total-effect estimate $ \widehat S_{T_i}=\widehat{\operatorname{E}[(Y_A-Y_{A_B^{(i)}})^2]}/(2\widehat{\operatorname{Var}(Y)}) $. Here $A_B^{(i)}$ is the hybrid vector with input $i$ taken from $B$ and the other inputs from $A$. The variance estimate and the numerator averages must be calculated according to a consistent convention. Finite-sample first-order estimates can fall below zero or above one; this is a warning about estimator variability, not evidence that the exact index is negative or exceeds one.

Three Complementary Diagnostics

No single diagnostic answers every question. Independent replicate designs reveal how much the estimate changes when the whole design is redrawn. A row bootstrap estimates sampling variability by resampling the observed design blocks. A sample-size convergence plot shows how estimates and their variability change as more base rows are used.

Key distinction: A Monte Carlo interval describes variability of an estimator under repeated sampling, conditional on the model, input distributions, estimator, and sample size. It is not an interval for uncertainty about the true model or a probability statement that the index lies in a particular interval.

Independent Replicates and Ranking Uncertainty

Run $R$ independent Saltelli designs of the same size $N$, calculate the index separately in each, and denote the estimates by $\widehat S_i^{(1)},\ldots,\widehat S_i^{(R)}$. Their average is $ \overline S_i=R^{-1}\sum_{r=1}^R\widehat S_i^{(r)} $. If the replicate estimates have sample standard deviation $s_{\mathrm{rep}}$, the estimated standard error of their average is $s_{\mathrm{rep}}/\sqrt R$. This estimates the Monte Carlo uncertainty in the average of $R$ replicate estimates. The replicate spread itself, $s_{\mathrm{rep}}$, describes typical variation among estimates from a single design of size $N$.

Replicates should use independent random designs. Reusing the same base draws with only a different ordering does not provide independent replicates. Also, a small number of replicates gives a noisy estimate of the replicate spread; report $R$, $N$, and whether the reported estimate is one design or an average across designs.

Worked Example: Two Inputs with an Uncertain Ranking

Suppose five independent hypothetical designs, each with $N=2{,}000$ base rows, give first-order estimates for inputs A of 0.41, 0.43, 0.39, 0.42, and 0.40. For input B they give 0.40, 0.41, 0.42, 0.39, and 0.41. Treat these as separate replicate sets for the two indices.

For A, the mean estimate is $ (0.41+0.43+0.39+0.42+0.40)/5=0.410 $. The squared deviations from 0.410 sum to 0.00100, so the replicate sample variance is $0.00100/(5-1)=0.000250$, and the sample standard deviation is $s_A=\sqrt{0.000250}=0.0158$. The standard error of the mean across five replicates is $0.0158/\sqrt5=0.00707$.

For B, the mean is $0.406$. Its squared deviations sum to 0.000520, giving $s_B=\sqrt{0.000520/4}=0.0114$ and standard error $0.0114/\sqrt5=0.00510$. The estimated difference is $0.410-0.406=0.004$. Under the simplifying assumption that the replicate sets are independent, its standard error is approximately $ \sqrt{0.00707^2+0.00510^2}=0.00872 $.

The estimated difference is smaller than its own Monte Carlo standard error. These results do not establish a reliable ranking of A above B. They also do not establish that the two effects are practically equivalent: that requires a prespecified practical margin and an uncertainty assessment aimed at that question. This calculation is a diagnostic, not a formal claim that the index difference is normally distributed.

Bootstrap Intervals That Preserve the Design

A nonparametric row bootstrap resamples the $N$ row blocks with replacement. For every resampled row, keep together the outputs for $A$, $B$, and all relevant hybrid evaluations. Recompute the variance estimate and Sobol index from each bootstrap sample. Repeating this procedure $B$ times gives bootstrap estimates; their chosen quantiles, such as the 2.5th and 97.5th percentiles, form a percentile 95% bootstrap interval.

Resampling individual model evaluations separately would break the Saltelli design relationships and is generally inappropriate. The bootstrap also has limits: it approximates the sampling behavior represented in the observed rows, does not remove bias in the estimator, and can be unreliable when there are very few rows, extreme outputs, or highly unstable estimates. A percentile interval crossing zero for a first-order estimate does not mean the exact index is negative; it signals that the finite-sample estimate does not resolve a positive first-order contribution well.

Worked Example: A Small First-Order Estimate with a Larger Total Effect

In a hypothetical QSP model, input C has estimated first-order index $\widehat S_C=0.018$. A bootstrap with $B=2{,}000$ resampled row blocks gives a 95% percentile interval from 0.000 to 0.061. The estimated total-effect index is $\widehat S_{T_C}=0.24$, with bootstrap interval 0.19 to 0.30. Suppose the analysis plan specified 0.05 as the threshold for a practically small first-order contribution.

The first-order point estimate is below the threshold because $0.018<0.05$, but the interval includes values above it. Thus the estimate does not give a precise conclusion that C's first-order contribution is below the threshold. More importantly, the total-effect estimate is much larger: its interval is entirely above 0.05. C's effect may therefore arise substantially through interactions with other inputs.

The defensible conclusion is that C's estimated first-order contribution is small, but its overall contribution is not shown to be negligible; the total-effect result indicates a meaningful contribution under the specified input distribution. A small first-order index alone cannot justify excluding C from the model or sensitivity analysis.

Sample-Size Convergence

A convergence plot displays an index estimate against the number $N$ of base rows. A useful version shows estimates from several independent designs at each $N$, or shows an estimate with an uncertainty band. A single cumulative sequence can be useful, but its successive points are dependent because they reuse earlier rows; apparent smoothness in such a plot is not independent confirmation of convergence.

Under independent sampling and suitable finite-moment conditions, Monte Carlo error typically decreases at a rate proportional to $N^{-1/2}$. This is a guide, not a guarantee that each estimate will improve monotonically. A stable-looking curve can still be biased by the estimator or affected by model-solver error. Compare estimates and uncertainty across increasing sample sizes, and check that the model evaluations themselves are numerically reliable.

Worked Example: Reading a Sample-Size Convergence Plot

For a hypothetical first-order index, five independent designs at each sample size produce the following means and between-design standard deviations: at $N=500$, mean 0.320 and SD 0.060; at $N=1{,}000$, mean 0.300 and SD 0.040; at $N=2{,}000$, mean 0.290 and SD 0.025; and at $N=4{,}000$, mean 0.295 and SD 0.018. A second input's corresponding estimates are near 0.28 at the two largest sample sizes.

The between-design spread shrinks as $N$ increases, while the mean for the first index moves from 0.320 to values near 0.29. The $N^{-1/2}$ guide predicts that increasing $N$ from 500 to 4,000, an eightfold increase, would reduce a standard deviation of 0.060 to approximately $0.060\sqrt{500/4{,}000}=0.0212$. The observed 0.018 is of a similar order, consistent with substantial Monte Carlo variability reduction.

At $N=4{,}000$, the difference between the two inputs' estimates is only about $0.295-0.28=0.015$, while the first index's between-design SD is 0.018. The plot therefore does not support a confident ordering of these inputs. More sampling may help, but the decision to spend that computation should depend on whether the ranking matters and on the required precision, not on a desire for a visually smooth curve.

Choosing and Interpreting Diagnostics

Approach	What it helps assess	Main caution
Independent replicate designs	Variation across wholly redrawn designs; stability of rankings	Requires multiple designs and enough replicates to estimate spread usefully
Row bootstrap	Sampling uncertainty conditional on the observed design, via resampled row blocks	Cannot correct estimator bias or guarantee reliable intervals in difficult samples
Sample-size convergence plot	Whether estimates and variability settle as the number of base rows grows	A smooth curve alone does not prove accuracy; cumulative points are dependent

These methods are complementary rather than interchangeable. Replicates are especially transparent for checking rankings across complete independent designs. A bootstrap can provide an interval from one sufficiently large design, but its assumptions and resampling unit must be respected. Convergence plots expose changes with sample size, but do not by themselves provide a calibrated confidence level.

For stochastic simulators, distinguish Monte Carlo variation from randomness in the model output itself. Decide whether the target is sensitivity of a deterministic expected output, sensitivity of a particular stochastic realization, or another explicitly defined quantity. Use the same target consistently across designs; otherwise differences between estimates can mix design error with simulation noise.

Common Pitfalls and Practical Decisions

Calling a close ranking certain. A larger point estimate is not automatically a meaningfully larger effect. Examine the uncertainty in the difference or compare rankings across independent replicates.
Calling an uncertain estimate negligible. An interval that includes small values may also include values large enough to matter. State the prespecified practical threshold and compare the interval with it.
Equating a small first-order index with a negligible input. First-order indices omit interaction contributions. Check the total-effect index or other appropriate interaction diagnostics before concluding that an input has little overall effect.
Bootstrapping evaluations independently. Resample complete Saltelli row blocks so that each design's paired evaluations remain together.
Confusing Monte Carlo precision with scientific certainty. A narrow interval can coexist with a misspecified input distribution, model structural error, or solver error. The interval addresses none of those automatically.

A practical report should specify the estimator, $N$, the number of independent replicates or bootstrap resamples, the uncertainty summary, and the input distribution. It should distinguish a stable index estimate from a practically important effect, and a small first-order contribution from a small total effect. If the intended decision depends on ordering or a threshold, state that decision rule before interpreting the diagnostics.

Key takeaway: Monte Carlo diagnostics quantify how much Sobol estimates may change because of finite sampling. Use replicate variability, design-preserving bootstrap intervals, and sample-size convergence together when decisions warrant it. Interpret first-order and total-effect indices separately: low first-order contribution does not rule out important interactions.

Check Your Understanding

Answer each question, then compare with the concise solutions below.

What is the independent resampling unit in a bootstrap of Saltelli estimates, and why?
Five independent replicate estimates have sample standard deviation 0.020. What is the estimated standard error of their mean?
An input has a first-order estimate of 0.02 and a total-effect estimate of 0.31. Can the input be called negligible overall? Explain.
Why does a smooth estimate from a single cumulative sample-size sequence not independently demonstrate convergence?

Solutions. (1) Resample complete row blocks, including the paired $A$, $B$, and hybrid evaluations, because the evaluations within a Saltelli row are structurally related. (2) The standard error is $0.020/\sqrt5=0.00894$. (3) No. The first-order estimate is small, but the much larger total-effect estimate is consistent with an important contribution involving interactions; inspect its uncertainty as well. (4) Successive cumulative estimates reuse earlier rows and are dependent, so a smooth trajectory is not a sequence of independent confirmations.