Latin Hypercube Sampling

Latin Hypercube Sampling in Probabilistic Sensitivity Analysis

In a health economic model, Latin hypercube sampling is a rule for choosing the parameter sets on which a probabilistic sensitivity analysis is run, so that every uncertain input is drawn evenly across its whole distribution instead of wherever independent draws happen to fall. The model is then run once per parameter set as in ordinary Monte Carlo simulation, and expected costs, QALYs and net benefit are averages over the runs. The gain is precision, which matters most when each run of a patient-level or dynamic model is slow. The design also changes how simulation error is estimated and how correlated inputs are handled.

Why the choice of draws matters

The NICE manual (PMG36, section 4.7.10) asks that distributions characterise uncertainty about mean parameter values and represent the available evidence. Section 4.7.11 adds that enough simulations should be used to minimise the effect of Monte Carlo error, and section 4.7.12 prefers a probabilistic estimate when possible unless the model is linear. Every such model needs a way of turning input distributions into a finite list of parameter sets.

The default is simple random sampling, in which each run draws every input independently, and with few runs the draws cluster by chance. With five runs and an input cut into five equally likely intervals, the probability that independent draws land once in every interval is 5!/5^5, or 120/3,125 = 0.0384, so about 96% of such samples miss at least one interval. Latin hypercube sampling removes this chance element for every input at once, and is one of the variance reduction techniques.

Building the sample from equal-probability strata

McKay, Beckman and Conover introduced the method in 1979 for computer codes too slow to run many times. They divided the range of each input into n strata of equal marginal probability 1/n, sampled once from each stratum and matched the inputs at random. They noted its advantage when the output is dominated by a few inputs, since each of those is fully stratified whichever turns out to be important.

$$x_{ij} = F_j^{-1}\left(\frac{\pi_j(i-1) + U_{ij}}{n}\right), \qquad i = 1, \ldots, n, \quad j = 1, \ldots, d$$

where $x_{ij}$ is the value of input $j$ in run $i$, $F_j^{-1}$ is the inverse cumulative distribution function of that input, $\pi_j(i-1)$ is the stratum, a member of $\{0, 1, \ldots, n-1\}$, given to run $i$ by a random permutation drawn independently for each input, $U_{ij}$ is an independent uniform draw on $[0, 1)$ placing the value within its stratum, $n$ is the number of runs and $d$ is the number of inputs.

Owen shows that each point keeps the correct uniform distribution, so the mean model output is unbiased. Passing each coordinate through the quantile function carries this over to any marginal, such as a beta distribution for a probability or a gamma distribution for a cost. Owen also notes that n points balance strata in all d inputs at once, that the method works even with more inputs than runs, and that a centred version using stratum midpoints has a small bias.

How much precision the method buys

McKay and colleagues proved that if the output is monotonic in each input, and the summary is a monotonic function of the output, the variance is no larger than under random sampling. Owen reports Stein's result that for large samples the additive part of the model drops out of the variance, and his own bound that the method cannot be much worse than random sampling.

$$\operatorname{Var}\left(\hat{\mu}_{\mathrm{SRS}}\right) = \frac{\sigma^2}{n}, \qquad \operatorname{Var}\left(\hat{\mu}_{\mathrm{LHS}}\right) \approx \frac{\sigma_e^2}{n}, \qquad \operatorname{Var}\left(\hat{\mu}_{\mathrm{LHS}}\right) \le \frac{\sigma^2}{n-1}$$

where $\hat{\mu}_{\mathrm{SRS}}$ and $\hat{\mu}_{\mathrm{LHS}}$ are mean outputs from a simple random sample and a Latin hypercube sample of $n$ runs, $\sigma^2$ is the output variance, $\sigma_e^2$ is the variance left after removing the best additive approximation (a sum of one-input functions), and the bound holds for $n \ge 2$.

Where a decision model's output is largely a sum of single-input effects, the method helps most. Interactions, such as a response probability multiplied by the QALY gain per responder, are the part stratification cannot remove. That non-additive share is one minus the sum of the first-order Sobol indices used in global sensitivity analysis. Helton and Davis add that if a model is cheap to run, a large sample can be used and the choice of design has little effect on cost or results.

Worked example with a three-input net benefit model

All numbers in this example are illustrative. Incremental net monetary benefit is modelled as $\mathrm{INMB} = \lambda p g-C$, with $\lambda$ = GBP 25,000 per QALY, one of the two values PMG36 (section 4.2.16) uses for net health benefit. The inputs are uniform so every quantile can be checked by hand. They are the response probability $p$ between 0.20 and 0.60, the QALY gain per responder $g$ between 0.50 and 1.00, and the incremental cost $C$ between GBP 1,000 and GBP 5,000. The true expected INMB is 25,000 × 0.40 × 0.75 − 3,000 = GBP 4,500.

1. Cut each range into five strata. With $n = 5$, each probability scale is cut at 0.2, 0.4, 0.6 and 0.8, and a uniform input on $[a, b]$ has quantile function $a+(b-a)u$.

2. Permute the strata and place each value. Each input gets its own random ordering of strata, written $\pi$, and its own position $U$. For run 1, $p$ falls in stratum 2 with $U = 0.30$, so $u = (2 + 0.30)/5 = 0.46$ and $p = 0.20 + 0.40 × 0.46 = 0.384$.

3. Check the stratification. Each stratum column holds 0, 1, 2, 3 and 4 exactly once.

4. Run the model. For run 1, 25,000 × 0.384 × 0.96 = 9,216, and subtracting GBP 2,000 leaves GBP 7,216. The mean of the five runs is (7,216 + 168 + 4,588 + 5,400 + 4,266)/5 = 21,638/5 = GBP 4,327.60, against the true GBP 4,500.

5. Compare precision at 1,000 runs. A uniform input has variance equal to its squared range over 12, giving 0.16/12 = 0.013333 for $p$, 0.25/12 = 0.020833 for $g$ and 16,000,000/12 = 1,333,333 for $C$. With 625,000,000 as the square of 25,000, the output variance splits into main effects of 625,000,000 × 0.5625 × 0.16/12 = 4,687,500 for $p$ (0.5625 is the squared mean of $g$), 625,000,000 × 0.16 × 0.25/12 = 2,083,333 for $g$ (0.16 is the squared mean of $p$) and 1,333,333 for $C$, plus an interaction of 625,000,000 × 0.16/12 × 0.25/12 = 173,611. The total is 8,277,778. A random sample of 1,000 runs has standard error $\sqrt{8{,}277{,}778/1{,}000}$ = GBP 91.0, while for a Latin hypercube of the same size only the interaction remains to a good approximation, giving $\sqrt{173{,}611/1{,}000}$ = GBP 13.2.

6. Interpret. The ratio 8,277,778/173,611 = 47.7 means 1,000 Latin hypercube runs give about the precision of 47,700 independent runs for this model. The gain is large because only 2.1% of the output variance is non-additive, and it would shrink in a model dominated by interactions. Both designs estimate the same expected INMB, on which the decision rests.

Correlated inputs and restricted pairing

Random matching of strata leaves sample correlations near zero but not exactly zero, and many models need deliberate correlation. PMG36 (section 4.7.19) asks that evidence on correlation between parameters be reflected in the probabilistic analysis and the assumptions presented. Helton and Davis describe the restricted pairing method of Iman and Conover, which induces a target rank correlation structure while preserving the marginal distributions.

The procedure transforms a matrix of van der Waerden scores using the Cholesky decomposition of the target correlation matrix, then reorders the sampled values in each column to follow the ranks of the transformed scores. Only the pairing changes, so each input keeps one draw per stratum. Helton and Davis advise checking that the resulting Spearman rank correlations are close to the target and repeating the construction if not.

Estimating Monte Carlo error from replicated hypercubes

Draws within one hypercube are not independent, so the random-sample formula, the output standard deviation over the square root of n, tends to overstate the error for a nearly additive model. Owen notes that the variance is harder to estimate than under independent sampling and that independent replicates are the most direct route. Helton and Davis describe Iman's replicated sampling procedure, which regenerates the hypercube with different random seeds.

$$\bar{T} = \frac{1}{n_R} \sum_{r=1}^{n_R} T_r, \qquad \mathrm{SE}\left(\bar{T}\right) = \sqrt{\frac{\sum_{r=1}^{n_R} \left(T_r-\bar{T}\right)^2}{n_R \left(n_R-1\right)}}$$

where $T_r$ is the estimate from replicate $r$, $n_R$ is the number of replicates and $\bar{T}$ is the pooled estimate, with a confidence interval from the t distribution on $n_R-1$ degrees of freedom.

Continuing the illustrative example, ten replicates of 100 runs gave mean INMB values of 4,527, 4,432, 4,540, 4,467, 4,429, 4,517, 4,506, 4,514, 4,555 and 4,464. Their sum is 44,951, so the pooled mean is GBP 4,495.1. The squared deviations sum to 17,684.9, so $\mathrm{SE} = \sqrt{17{,}684.9/90}$ = GBP 14.0, and with $t = 2.262$ the 95% interval is 4,495.1 ± 31.7, or GBP 4,463.4 to GBP 4,526.8. The same check gives one way to judge whether enough simulations have been run, as PMG36 (section 4.7.11) asks.

Uses beyond probabilistic sensitivity analysis

The same designs choose the runs used to fit approximations of slow models. Degeling and colleagues, introducing metamodelling for health economic models, judge Latin hypercube designs preferable for training and testing data when only a limited number of runs is affordable. They describe maximin versions, which maximise the minimum distance between design points, and warn that holding out a random share of such a design for validation breaks its coverage, so a separate validation set should be simulated. The metamodels they review include Gaussian process regression, and they list model calibration and value of information analysis among earlier health economic uses of metamodels.

Helton and Davis describe analyses that use Latin hypercube sampling for uncertainty about fixed quantities and random sampling for stochastic variation, conditional on each sample element. The health economic analogue is a microsimulation with an outer loop of parameter sets and an inner loop of simulated patients. TreeAge Pro offers the method as an option for probabilistic analysis but samples distributions independently by default.

Where it sits among sampling methods

Neighbouring techniques target different parts of simulation error. The right one depends on the output needed.

• Grid stratification of the joint space. Owen notes that a regular grid gives only $n^{1/d}$ strata per input and becomes infeasible with many inputs, whereas a Latin hypercube stratifies each input separately.
• Sobol sequences. These deterministic points from a Sobol sequence fill the joint space evenly, and Degeling and colleagues note they may be more efficient than Latin hypercube designs in low-dimensional problems.
• Common random numbers. The common random numbers method aligns random streams across strategies to sharpen incremental results, a different target from spreading parameter draws.
• Stratified sampling for tails. Helton and Davis note that stratified sampling can suit the estimation of extreme quantiles better.

Limitations and common misreadings

The method reduces simulation noise only, and does nothing for parameter uncertainty itself, for distributions that misrepresent the evidence or for structural choices.

1. Expecting a large gain in every model. The gain is set by how additive the model is.
2. Using the independent-sample standard error. Replicated hypercubes give an honest error estimate.
3. Topping up a finished sample. A hypercube is built for a fixed n, so extra random runs lose the stratification, while a further replicate keeps it.
4. Reading even spread as more information. The decision uncertainty shown in a cost-effectiveness acceptability curve is the same under either design, and only the precision of its estimate changes.

Sources

• Degeling K, IJzerman MJ, Lavieri MS, Strong M, Koffijberg H. Introduction to metamodeling for reducing computational burden of advanced analyses with health economic models: a structured overview of metamodeling methods in a 6-step application process. Medical Decision Making. 2020;40(3):348-363. https://doi.org/10.1177/0272989X20912233
• Helton JC, Davis FJ. Latin Hypercube Sampling and the Propagation of Uncertainty in Analyses of Complex Systems. SAND2001-0417. Sandia National Laboratories; November 2002. https://doi.org/10.2172/806696
• McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 1979;21(2):239-245. Reprinted in Technometrics. 2000;42(1):55-61. https://doi.org/10.2307/1268522
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Sections 4.2.16, 4.7.10 to 4.7.12 and 4.7.19. https://www.nice.org.uk/process/pmg36
• Owen AB. Monte Carlo Theory, Methods and Examples. Draft chapters 2009 to 2013. Chapter 10, Advanced variance reduction, section 10.3 and end notes. https://artowen.su.domains/mc/Ch-var-adv.pdf
• TreeAge Software. TreeAge Pro user manual, section 32.14, Latin Hypercube Sampling. Accessed 26 October 2026. https://www.treeage.com/help/Content/66-Probabilistic-Sensitivty-Analysis-on-CE-Models/12-Latin-Hypercube.htm