Monte Carlo Integration

Monte Carlo Integration in Health Economic Models

The expected cost per patient, the expected QALYs and the expected net benefit that a health economic model reports are integrals over the joint distribution of its uncertain inputs. Monte Carlo integration evaluates such an integral by drawing input sets at random from that distribution, running the model for each and taking the average. Its error shrinks with the number of draws at a rate that does not depend on how many parameters the model has, so it copes with many uncertain inputs where grids fail. The same averaging sits behind probabilistic sensitivity analysis, patient-level simulation and value of information analysis.

Expected costs and QALYs written as integrals

A model output depends on a vector of uncertain inputs, each given a distribution that represents the evidence. The expected value of the output is its probability-weighted average over every possible input combination. For continuous inputs it is written as follows.

$$I = \mathbb{E}_{\theta}\left[g(\theta)\right] = \int g(\theta)\, p(\theta)\, d\theta$$

where $I$ is the quantity sought, $\theta$ is the vector of uncertain inputs, $g(\theta)$ is the model output for one input set (such as discounted cost or QALYs per patient) and $p(\theta)$ is the joint probability density of the inputs.

An integral that is not already an expectation can be turned into one. Over a finite interval, Owen's textbook writes the integral as the interval length times the mean of the integrand at a uniformly distributed point.

$$\int_a^b f(t)\, dt = (b - a)\, \mathbb{E}\left[f(U)\right], \qquad U \sim \text{Uniform}(a, b)$$

where $f(t)$ is the integrand, $a$ and $b$ are the limits and $U$ is a random point spread evenly over the interval. The area under a survival function and discounted cost accumulated over a time horizon fit this form as readily as an expectation over parameters.

The estimator, its error and its standard error

The estimate replaces the integral with an average over independent draws. Owen calls this simple Monte Carlo and shows that it is unbiased, so its expected value equals the target whatever the number of draws.

$$\hat{I}_n = \frac{1}{n} \sum_{i=1}^{n} g(\theta_i), \qquad \operatorname{Var}\left(\hat{I}_n\right) = \frac{\sigma^2}{n}$$

where $\hat{I}_n$ is the estimate from $n$ draws, $\theta_i$ is the $i$-th input set drawn independently from $p(\theta)$ and $\sigma^2$ is the variance of the output $g(\theta)$ across draws. The law of large numbers guarantees that the average settles on $I$ as $n$ grows, provided the expectation exists. When $\sigma^2$ is finite, the root mean squared error is $\sigma / \sqrt{n}$, and the central limit theorem gives an approximate 95% interval of $\hat{I}_n \pm 1.96\, s / \sqrt{n}$, where $s$ is the sample standard deviation of the outputs.

Owen spells out the exchange rate. One more decimal digit of accuracy needs a tenth of the error and so a hundredfold increase in computation, which makes simple Monte Carlo poorly suited to problems that need high precision. The draws themselves supply a rough error estimate, the Monte Carlo error checked when judging whether enough draws have been run.

Why the number of parameters does not enter the error

The formula $\sigma / \sqrt{n}$ contains no term for the dimension of $\theta$. Owen stresses that the dimension can be 2 or 1,000 and the root mean squared error is still $\sigma / \sqrt{n}$, whereas a straightforward extension of Simpson's rule to $d$ dimensions has error of order $n^{-4/d}$, which he calls useless for large $d$. In one dimension the comparison runs the other way, since Simpson's rule achieves error of order $n^{-4}$ for an integrand with a continuous fourth derivative.

Metropolis and Ulam gave the classic illustration in their 1949 paper on the method. Finding the volume of a region in a twenty-dimensional unit cube by dividing each coordinate into ten parts would mean examining ten to the power twenty lattice points, which they called obviously impossible, while about 10,000 random points would give an estimate valid within a few per cent with high probability. A probabilistic model with twenty uncertain inputs poses the same problem, and Owen shows that in worst-case terms any quadrature rule in high dimensions suffers this curse of dimensionality. Numerical integration by quadrature rules such as Gaussian quadrature therefore suits survival areas and other smooth low-dimensional integrals, while Monte Carlo suits expectations over many parameters.

Probabilistic sensitivity analysis as an integral over parameters

When a model is linear in its inputs, one run at the mean of every parameter returns the expected output. Otherwise that run is in general biased, and the expected output must be found by integrating over the parameter distribution. Doubilet and colleagues described such simulation for medical decision trees in 1985, assigning every probability and utility a random value from its distribution in each simulation and recording the mean and standard deviation of each strategy's expected utility.

The NICE manual (PMG36, section 4.7.12) states that the committee's preferred cost-effectiveness estimate should be derived from a probabilistic analysis when possible unless the model is linear. Section 4.7.11 asks that enough model simulations be used to minimise the effect of Monte Carlo error, and suggests reviewing the variance around net benefits or ICERs as the number of simulations increases. A Monte Carlo simulation of a cohort model studies the whole distribution of outputs, and its mean is a Monte Carlo integral whose integrand is the whole model, with one full run per draw. The probability that a strategy is cost-effective is also such an integral, the average of an indicator.

Patient-level models nest a second integral inside the first. In a microsimulation, each run averages over simulated patients for one parameter set, and the probabilistic analysis averages those runs over parameter sets. The ISPOR-SMDM task force report on parameter estimation (Briggs and colleagues, 2012) notes that assessing parameter uncertainty in such models requires the stochastic uncertainty, sometimes called Monte Carlo error, to be eliminated. O'Hagan, Stevenson and Madan noted that each run may simulate many thousands of patients, so only a relatively small number of runs is possible, and gave analysis-of-variance methods with formulae for optimal sample sizes.

Worked example with restricted mean survival over ten years

The numbers in this example are illustrative. A model needs the restricted mean survival time over a ten-year horizon for patients whose survival follows an exponential distribution with a hazard of 0.2 per year. The target has a closed form, so each Monte Carlo step can be checked against it.

$$\text{RMST}(\tau) = \int_0^{\tau} e^{-\lambda t}\, dt = \frac{1 - e^{-\lambda \tau}}{\lambda}$$

where $\tau = 10$ years is the horizon, $\lambda = 0.2$ per year is the hazard and $e^{-\lambda t}$ is the probability of surviving to time $t$. The exact value is $(1 - e^{-2}) / 0.2 = 0.864665 / 0.2 = 4.3233$ life-years.

1. Write the integral as an expectation. With $U$ uniform on 0 to 10 years, the restricted mean equals $\mathbb{E}[10\, e^{-0.2U}]$, so each draw at time $t$ contributes $(b - a)\, f(t) = 10\, e^{-0.2 t}$, with $f(t) = e^{-0.2 t}$ as the integrand.

2. Take five draws. Suppose the random times are 1.2, 4.7, 8.1, 3.3 and 6.5 years. The survival probabilities are 0.7866, 0.3906, 0.1979, 0.5169 and 0.2725, and the contributions are 7.8663, 3.9063, 1.9790, 5.1685 and 2.7253.

3. Average. The contributions sum to 21.6454, so the estimate is 21.6454 / 5 = 4.3291 life-years, within 0.006 of the exact value.

4. Estimate the standard error. The sample standard deviation of the five contributions is 2.3179, so the standard error $s / \sqrt{n}$ is 2.3179 / 2.2361 = 1.0366 life-years. The closeness in step 3 is luck, since five draws cannot be relied upon to land so near.

5. Scale up. The mean squared contribution is $\tfrac{10 (1 - e^{-4})}{0.4} = 24.542$, so the true variance of one contribution is 24.542 minus 4.3233 squared, about 5.851, and its standard deviation about 2.419. With 10,000 draws the standard error is 2.419 / 100 = 0.0242 life-years, and a standard error of 0.01 would need about 5.851 / 0.0001, or roughly 58,500 draws.

6. Change the sampling density. A patient-level simulation estimates the same integral by drawing survival times $T$ from the exponential distribution itself and averaging the time lived within the horizon, $\min(T, 10)$. Its mean square is $\tfrac{2 (1 - 3e^{-2})}{0.04} = 29.700$, so its variance is about 11.009, about 1.88 times the uniform version, and it needs roughly 110,100 draws for a standard error of 0.01. Both estimators are unbiased, and the choice of sampling density changes only their precision.

7. Compare a deterministic rule. The trapezoidal rule uses six survival values, 1.0000, 0.6703, 0.4493, 0.3012, 0.2019 and 0.1353, and with step $h = 2$ years gives $h \times (0.5 + 1.6227 + 0.0677) = 4.3808$. Its error of 0.0575 is fixed and far smaller than the five-draw standard error, so for a smooth one-dimensional integral quadrature wins. The advantage reverses once the expectation runs over many uncertain parameters at once.

Nested integrals in value of information

The EVPI is a difference between two integrals over the same parameter distribution, and both are estimated from the same probabilistic draws. The first averages a maximum and the second maximises an average, so their gap is the expected loss from choosing before the uncertainty is resolved.

$$\text{EVPI} = \mathbb{E}_{\theta}\left[\max_j NB_j(\theta)\right] - \max_j \mathbb{E}_{\theta}\left[NB_j(\theta)\right]$$

where $NB_j(\theta)$ is the net benefit of strategy $j$ for parameter set $\theta$.

The partial value for a subset of parameters, the EVPPI, places one expectation inside a maximum inside another expectation. Brennan, Kharroubi, O'Hagan and Chilcott showed that its two nested expectations must be evaluated separately because of the maximum between them, and set out a two-level algorithm with an outer loop over the parameters of interest and an inner loop over the rest. They also showed that maxima of Monte Carlo estimates of expectations are biased upward, so small samples give biased EVPPI estimates. Monte Carlo error therefore enters value of information as bias as well as noise.

Sampling choices that change precision but not the target

Because $\sigma^2$ sets the error, Owen notes that halving it while leaving the target unchanged gains as much as doubling $n$. Variance reduction methods such as common random numbers across strategies and control variates work on this principle. Importance sampling draws from a different density and reweights, and Owen warns that it sometimes brings enormous variance reductions and in other circumstances delivers an estimate with infinite variance.

Other methods change how points are placed. Quasi-Monte Carlo uses evenly spread deterministic points such as a Sobol sequence, and Owen notes that its error is hard to estimate unless the points are randomised so that independent replicates can be compared. Latin hypercube sampling stratifies the range of each input. Markov chain Monte Carlo, including the Metropolis-Hastings algorithm, applies the same averaging when independent draws cannot be made, sampling instead from a Markov chain whose limiting distribution is the one wanted. Its draws are correlated, so the $s / \sqrt{n}$ formula, which assumes independence, does not give its standard error.

Limits and common misreadings

The commonest errors come from treating the Monte Carlo average as more than an estimate of one stated integral. Each is avoided by naming the integral a reported number estimates.

1. Taking precision for accuracy. More draws shrink Monte Carlo error only. Owen points out that model errors can be larger than Monte Carlo errors, and no number of draws corrects a wrong distribution or structure.
2. Averaging ratios. Owen notes that the expectation of a ratio may not exist when the denominator has positive density at zero, and that ratios are ordinarily estimated as a ratio of sample averages. An ICER is therefore formed from mean incremental cost and mean incremental QALYs, not by averaging per-draw ratios whose QALY gain can be near zero.
3. Assuming the expectation exists. Owen shows that when the mean is infinite the running average keeps jumping and the usual interval fails.

Sources

• Brennan A, Kharroubi S, O'Hagan A, Chilcott J. Calculating partial expected value of perfect information via Monte Carlo sampling algorithms. Medical Decision Making. 2007;27(4):448-470. https://doi.org/10.1177/0272989X07302555
• Briggs AH, Weinstein MC, Fenwick EAL, Karnon J, Sculpher MJ, Paltiel AD. Model parameter estimation and uncertainty analysis: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force Working Group-6. Medical Decision Making. 2012;32(5):722-732. https://doi.org/10.1177/0272989X12458348
• Doubilet P, Begg CB, Weinstein MC, Braun P, McNeil BJ. Probabilistic sensitivity analysis using Monte Carlo simulation: a practical approach. Medical Decision Making. 1985;5(2):157-177. https://doi.org/10.1177/0272989X8500500205
• Metropolis N, Ulam S. The Monte Carlo method. Journal of the American Statistical Association. 1949;44(247):335-341. https://doi.org/10.1080/01621459.1949.10483310
• 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.7.11 and 4.7.12. https://www.nice.org.uk/process/pmg36/chapter/economic-evaluation-2
• O'Hagan A, Stevenson M, Madan J. Monte Carlo probabilistic sensitivity analysis for patient level simulation models: efficient estimation of mean and variance using ANOVA. Health Economics. 2007;16(10):1009-1023. https://doi.org/10.1002/hec.1199
• Owen AB. Monte Carlo Theory, Methods and Examples. Stanford University; 2013. Chapters 1, 2 and 7. https://artowen.su.domains/mc/