Control Variates

Control Variates in Health Economic Simulation Models

Control variates make a Monte Carlo estimate, such as mean discounted cost or QALYs from a patient-level simulation, more precise by removing the part of its random error that a second output reveals. That second output is computed from the same random draws and has an expected value already known from the model inputs, so the gap between its simulated mean and its known mean shows how far the run has drifted by chance. Shechter and colleagues noted that variance reduction methods are common in industrial engineering and operations research but seldom used in medical models, even though most Monte Carlo cohort simulations suit them well. Because the method works on the outputs after the run and leaves the sampling unchanged, it can be combined with other variance reduction methods, although an incorrect known mean biases the adjusted result.

How a known mean corrects a noisy estimate

Owen's textbook describes the idea as using a second problem, close to the one being solved, whose exact answer is known; in practice "known" can mean known with an error that is negligible next to the Monte Carlo error. If a simulated cohort happens to contain unusually long survivors, an output that tracks survival, such as undiscounted life years, will sit above its known mean, and the target output will usually be inflated too. Subtracting a scaled version of that gap pulls the target back towards its true value.

$$\hat{\mu}_{b} = \bar{Y} - b \left( \bar{X} - \mu_X \right), \qquad b^{*} = \frac{\operatorname{Cov}(Y, X)}{\operatorname{Var}(X)}$$

where $\hat{\mu}_{b}$ is the adjusted estimate of the expected model output, $\bar{Y}$ is the mean of the target output over $n$ independent replications, $\bar{X}$ is the mean of the control over the same replications, $\mu_X$ is the known expected value of the control, $b$ is a fixed coefficient and $b^{*}$ is the value of $b$ that minimises variance. Because the expected value of $\bar{X}$ equals $\mu_X$, the adjusted estimate is unbiased for any fixed $b$. Setting $b = 0$ gives the plain Monte Carlo mean and $b = 1$ gives the simple difference estimator.

$$\operatorname{Var}(\hat{\mu}_{b^{*}}) = \frac{\sigma_Y^{2}}{n} \left( 1 - \rho^{2} \right)$$

where $\sigma_Y^{2}$ is the variance of a single target output, $n$ is the number of replications and $\rho$ is the correlation between the target and the control. With the optimal coefficient the variance is never worse than the plain estimate. Owen notes that any correlated control helps, including one that is negatively correlated, since the coefficient takes the sign of the covariance.

Estimating the coefficient and the standard error

The optimal coefficient is rarely known, so it is estimated as the least squares slope of the target on the control from the same run. Owen shows that this makes the estimator slightly biased, but the bias shrinks in proportion to $n^{-1}$ while the random error shrinks in proportion to $n^{-1/2}$, so it is usually ignored. When an unbiased estimate is required, the coefficient can come from an independent pilot sample and then be applied to the main run.

The calculation fits ordinary linear regression software. The target is regressed on the control with its known mean subtracted, and the fitted intercept is the adjusted estimate, with the intercept's standard error as its standard error. The centring step matters, because only then does the intercept equal the expected value of the target. Several controls can enter the same regression, and the asymptotic variance is then multiplied by the factor $(1 - R^{2})$, where $R^{2}$ is the proportion of variance in the target explained by the controls. Rubinstein and Marcus quantified the loss in variance reduction caused by estimating the optimal coefficients from the run, and their numerical results showed good reductions when the number of controls was relatively small.

When the extra work pays

A control is worth using only if the precision gained outweighs the extra computing time. Owen gives the condition in terms of the average cost of one model evaluation and the extra cost of computing the control alongside it.

$$\left( 1 - \rho^{2} \right) \left( c_f + c_h \right) < c_f$$

where $c_f$ is the average cost of generating the inputs and evaluating the model once, $c_h$ is the extra average cost of the control and $\rho$ is the correlation defined above. Rearranged, the absolute correlation must exceed the square root of $c_h / (c_f + c_h)$. A control that costs as much as the model itself therefore needs an absolute correlation above about 0.71, and with $J$ controls each costing as much as the model, $R^{2}$ must exceed $J/(J+1)$. Controls that the model computes anyway, such as the sampled inputs or a survival time drawn on the way to a discounted outcome, cost almost nothing, so even a modest correlation helps.

Worked example with discounted life years in a discrete event simulation

The figures are illustrative. A discrete event simulation samples each patient's time to death from an exponential distribution with a hazard of 0.1 a year, by inverse transform sampling (using $U$ in place of one minus $U$, which has the same distribution), and the target output is life years discounted continuously at 0.035 a year. The control is the undiscounted survival time, whose expected value is known exactly from the input distribution as one divided by the hazard, or 10 years. The true expected discounted life years can also be found by integration as 1 / (0.1 + 0.035) = 7.407, which allows each estimate to be checked.

$$T = \frac{-\ln(U)}{h}, \qquad Y = \frac{1 - e^{-rT}}{r}$$

where $T$ is survival time in years, $U$ is a uniform random number, $h$ is the hazard of 0.1 a year, $Y$ is discounted life years and $r$ is the continuous discount rate of 0.035 a year.

1. Simulate five patients. Values are rounded to three decimal places.

2. Take the plain means. Mean survival is 38.837 / 5 = 7.767 years, well below the known 10, so this small run has drawn short survivors by chance. Mean discounted life years is 31.288 / 5 = 6.258, which understates the true 7.407 for the same reason.

3. Find the optimal coefficient and adjust. For this model the covariance of $Y$ and $T$ is $(h + r)^{-2} = 54.87$ and the variance of $T$ is $h^{-2} = 100$.

$$b^{*} = \frac{54.87}{100} = 0.5487, \qquad \hat{\mu}_{b^{*}} = 6.258 - 0.5487 \times (7.767 - 10) = 7.483$$

where $b^{*}$ is the optimal coefficient and $\hat{\mu}_{b^{*}}$ is the adjusted estimate of mean discounted life years. The error falls from 1.149 years to 0.076 years.

4. Compare a coefficient estimated from the sample. The least squares slope from the five patients is 0.717, which gives 6.258 + 0.717 × 2.233 = 7.859. Five points are far too few to estimate a coefficient, which is why it is taken from a run of thousands of patients or a separate pilot.

5. Compare precision at 1,000 patients. The variance of discounted life years for one patient is 32.28, so the plain estimate from 1,000 patients has a standard error of the square root of 0.03228, or 0.180 years. The correlation between $Y$ and $T$ is 0.966.

$$1 - \rho^{2} = 1 - 0.966^{2} = 0.067$$

where $\rho$ is the correlation between discounted and undiscounted survival. The standard error is multiplied by the square root of 0.067, which is 0.259, and falls to 0.180 × 0.259 = 0.047 years. Without the control, about 1 / 0.067 = 15 times as many patients would be needed for the same precision, and the control costs nothing extra because $T$ is computed anyway. A check of 4,000 repeated runs of 1,000 patients gave variances of 0.0317 for the plain mean and 0.0022 with the control, against 0.0323 and 0.0022 from the formulas. A cost accrued at a constant annual rate is a fixed multiple of $Y$, so discounted lifetime cost would gain in the same proportion.

Choosing a control in health economic models

A good control is cheap, strongly correlated with the target and has a mean that is known for the same population and input distribution as the simulation. A mean estimated from the same run cannot serve, because its error is the error the method is trying to remove.

In a PSA each iteration samples parameters from distributions whose means are known, such as a beta distribution for a probability, so the sampled parameters themselves can act as controls. Regressing net benefit on the centred parameter values and reading off the intercept corrects for chance imbalance in the parameter draws. Owen lists low order polynomials in the inputs as suitable inexpensive controls, and the gain depends on how well those terms predict the output.

In a patient-level model, the Decision Support Unit's Technical Support Document 15 (TSD 15) describes the sampling error between simulated patients as first-order uncertainty, which can be decreased by simulating more patients at the price of longer run times. Baseline characteristics sampled from a known distribution, such as age at entry, offer a control for this patient-level noise. Owen shows that post-stratification, which reweights a sample to known stratum shares, is a special case of control variates with stratum indicators as the controls. A simplified patient-level version of the model, such as one without history-dependent features, run on the same random numbers and with an exact expected value from a cohort calculation, is another option; Owen's finance example of an average-price option controlled by a closely related option with a closed-form value follows this pattern.

Interpretation in probabilistic and incremental analysis

Economic evaluation turns on differences between strategies, so the control is best applied to the incremental output, such as incremental net benefit, with a coefficient estimated for that difference. When strategies share common random numbers, a baseline characteristic enters both arms in the same way, and it reduces the variance of the difference only to the extent that the treatment effect varies with that characteristic. A control that sharpens each strategy's mean can therefore do little for the increment.

In a PSA the adjusted estimate corrects the mean across iterations. It does not narrow the spread of the iterations, which describes parameter uncertainty and must be kept for the cost-effectiveness acceptability curve and value of information analysis. Li and colleagues used metamodels for a different aim, fitting the PSA outcome as a function of the parameters plus an error term so that stochastic noise in each iteration could be removed while the parameter-driven variation was kept. The probability that a strategy is cost-effective is the mean of a binary indicator, and Owen notes that large reductions are hard to obtain for binary outcomes because few controls predict a binary value strongly. O'Hagan and colleagues observed that the many patients needed per run have led to a common view that PSA is impractical for patient-level models, which is where cheap controls offer the most.

Limitations and common errors

An incorrect known mean biases the result, since an error of $\delta$ in $\mu_X$ shifts the adjusted estimate by $b\delta$. The standard error must come from the regression residuals, not from the raw spread of the target, which describes the unadjusted estimator. Owen notes that a control variate estimate can fall outside natural limits, for example a negative value for a quantity that must be positive, and a probability can likewise be pushed above one. A coefficient fitted under one parameter set or strategy may not hold under another. The control must be computed from the same draws without altering the model, since forcing a correlation by changing patient pathways changes the quantity being estimated.

Control variates belong to the wider family of variance reduction methods used within Monte Carlo simulation. Antithetic variates, importance sampling and Latin hypercube sampling change how the random numbers are drawn, while control variates leave the sampling unchanged and work on the outputs after the run, so they can be combined with those methods. None of them reduces parameter or structural uncertainty, since a more precise estimate of what the model implies says nothing about whether the model and its evidence are right.

Sources

• Davis S, Stevenson M, Tappenden P, Wailoo A. NICE DSU Technical Support Document 15: Cost-effectiveness modelling using patient-level simulation. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; April 2014. Section 5.2. https://www.sheffield.ac.uk/nice-dsu/tsds/patient-level-simulation
• Li Z, Knowlton GS, Wheatley MM, Jenness SM, Enns EA. Meta-modeling as a variance-reduction technique for stochastic model-based cost-effectiveness analyses. Medical Decision Making. 2025;45(8):976-986. https://doi.org/10.1177/0272989X251352210
• 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. 2013. Chapter 8, Variance reduction, sections 8.4 and 8.9 and end notes. https://artowen.su.domains/mc/Ch-var-basic.pdf
• Rubinstein RY, Marcus R. Efficiency of multivariate control variates in Monte Carlo simulation. Operations Research. 1985;33(3):661-677. https://doi.org/10.1287/opre.33.3.661
• Shechter SM, Schaefer AJ, Braithwaite RS, Roberts MS. Increasing the efficiency of Monte Carlo cohort simulations with variance reduction techniques. Medical Decision Making. 2006;26(5):550-553. https://doi.org/10.1177/0272989X06290489