Antithetic Variates

Antithetic Variates for Lower Monte Carlo Error in Health Economic Models

Health economic models that rely on random sampling, such as probabilistic sensitivity analysis and patient-level simulations, give answers that carry Monte Carlo error, so a rerun with fresh random numbers moves the expected costs, QALYs and net benefit a little. Antithetic variates reduce that error by running the model on each set of random numbers and again on their mirror images, so that a high result from one run tends to be offset by a low result from its partner. The cancellation is guaranteed to reduce variance when the output moves in one direction with each input, and the method applies both to parameter draws in a PSA and to event times in microsimulation and discrete event simulation.

Why mirrored draws cancel simulation noise

Most simulation models start from uniform random numbers between 0 and 1 and transform them into the quantities the model needs, such as a transition probability or a time to an event. If a uniform draw $U$ produces an unusually long survival time, its complement $\tilde{U} = 1-U$ produces an unusually short one, and the average of the two sits closer to the true mean than either alone. For a normally distributed input the mirror image of a draw is its reflection about the mean. The method was introduced by Hammersley and Morton in 1956.

With $n$ model evaluations arranged as $n/2$ independent pairs, the antithetic estimate of an expected output and its variance are as follows.

$$\hat{\mu}_{\text{anti}} = \frac{1}{n} \sum_{i=1}^{n/2} \left[ f(X_i) + f(\tilde{X}_i) \right], \qquad \operatorname{Var}(\hat{\mu}_{\text{anti}}) = \frac{\sigma^2}{n} (1 + \rho)$$

where $\hat{\mu}_{\text{anti}}$ is the estimate of the expected model output, $f$ is the model viewed as a function of its random inputs, $X_i$ is the $i$th set of sampled inputs, $\tilde{X}_i$ is its mirror image, $\sigma^2$ is the variance of a single model output and $\rho$ is the correlation between $f(X_i)$ and $f(\tilde{X}_i)$. Ordinary sampling with $n$ independent evaluations has variance $\sigma^2/n$, so the factor $(1+\rho)$ measures the gain. Because $\rho$ lies between $-1$ and $+1$, the factor lies between 0 and 2; in the best case one pair gives the exact answer, and in the worst case the variance doubles.

A useful way to see when the correlation is negative is to split the output into a part that changes sign when the inputs are mirrored (the odd part) and a part that does not (the even part). Pairing removes the odd part entirely but doubles the contribution of the even part. Antithetic sampling therefore works best when the output is close to linear in its inputs and worst when it is symmetric, so that both members of a pair give the same answer.

Monotone outputs and the guarantee of a gain

The main theoretical result is that if the model output is monotone in every random input, the correlation within each pair is negative, so antithetic sampling cannot increase the variance. The direction can differ between inputs, because mirroring a single input turns a falling relationship into a rising one without changing the sampling. In many decision models incremental net benefit rises with the treatment effect and falls with the price of the new treatment, but the guarantee needs the output to be monotone in every random input, and a single input with a turning point removes it.

Monotonicity guarantees some gain, not a large one. The correlation can be only slightly below zero for a monotone output, and a non-monotone output can still have a strongly negative correlation. The practical test is a pilot run of the actual model that compares the variance of pair means with the variance of independent runs at the same number of model evaluations, for the quantity that matters to the decision.

Antithetic pairs in probabilistic sensitivity analysis

In a PSA of a cohort model, each iteration draws a value for every uncertain parameter and records the expected costs, QALYs and net monetary benefit for each strategy. Antithetic sampling replaces half of the independent iterations with mirrored partners. If a beta distribution for a response rate is sampled at its 90th percentile in one iteration, the partner samples it at its 10th percentile, and every other parameter is mirrored in the same way.

The mirroring is applied to the uniform numbers before they are transformed, which requires the inverse distribution function method, under which a parameter equal to $F^{-1}(U)$, where $F$ is its distribution function, has its partner at $F^{-1}(\tilde{U})$. Generators that use rejection steps do not map each uniform number to one parameter value in a monotone way, so feeding them the complement does not produce a mirrored draw. Correlated parameters sampled through a Cholesky decomposition can be mirrored by reversing the sign of the underlying standard normal draws, which reflects the whole parameter vector about its mean and keeps the correlation structure. The monotone guarantee then applies to the independent normal draws, not to the parameters, because when parameters are negatively correlated, an output that is monotone in each parameter need not be monotone in each underlying draw.

Each member of a pair is drawn from the correct distribution, so the full set of iterations still estimates mean outcomes, the probability that a strategy is cost-effective and the cost-effectiveness acceptability curve without bias. The iterations are not independent, however, so the Monte Carlo standard error of any of these estimates has to be computed from the $n/2$ pair averages, treating each pair as one observation.

Antithetic pairs in microsimulation and discrete event models

Patient-level models add a second layer of randomness. Each simulated patient's history depends on random numbers for event times and event types, so even with fixed parameter values the average cost and QALY per patient carry noise that falls only as more patients are simulated. The Decision Support Unit's Technical Support Document 15 notes that a PSA of a patient-level model usually needs two nested loops, an inner loop over simulated patients and an outer loop over parameter values. The cost of this nesting has led to a common view that PSA is not practical for patient-level models, which prompted methods based on analysis of variance that set the numbers of patients and parameter sets to simulate.

Antithetic pairing can be applied to the inner loop by simulating each patient twice, once with the original random numbers and once with their complements. The pairing only holds if both runs use each random number for the same purpose. In a discrete event model the number of random numbers consumed depends on which events occur, so a single shared stream soon falls out of step between the two runs; a separate stream for each patient, and for each stochastic process within a patient, keeps the draws aligned, the same discipline that common random numbers require.

A 2010 study that coded a type 2 diabetes simulation based on the UKPDS 68 outcomes equations reported that antithetic variates typically cut the number of simulation replications, and the run time, needed for predicted QALYs over a 40-year horizon by 53%.

Worked example with exponential survival times in a discrete event simulation

The following figures are illustrative. A discrete event simulation samples each patient's time to death from an exponential distribution with a hazard of 0.1 per year, using the inverse distribution function, and the quantity of interest is mean undiscounted life expectancy. The true answer is 10 years, which makes the precision of each method easy to check.

$$T = \frac{-\ln(U)}{h}, \qquad \tilde{T} = \frac{-\ln(1-U)}{h}$$

where $T$ is the sampled survival time in years, $\tilde{T}$ is the survival time of the antithetic partner, $U$ is a uniform random number and $h$ is the hazard of 0.1 per year.

1. Sample two pairs. With $U = 0.20$ the patient survives 16.094 years and the partner, drawn at 0.80, survives 2.231 years; the pair average is 9.163 years. With $U = 0.90$ the two survival times are 1.054 and 23.026 years, averaging 12.040 years. Each long survival time is matched with a short one, so the pair averages sit closer to 10 years than any of the four single draws.

2. Find the correlation within a pair. A single exponential survival time with hazard 0.1 has variance $h^{-2}$, or 100 years squared. The correlation between $-\ln(U)$ and $-\ln(1-U)$ can be found by integration and equals $(1 - \pi^2/6)$, or $-0.6449$ to four decimal places. Survival time is monotone in $U$, so a negative value was guaranteed.

3. Compare precision at 1,000 simulated patients. Ordinary sampling of 1,000 independent patients gives an estimate of mean survival with variance $\sigma^2/n = 100/1000 = 0.1$, a standard error of 0.316 years. Antithetic sampling with 500 pairs uses the same 1,000 model evaluations, as follows.

$$\operatorname{Var}(\hat{\mu}_{\text{anti}}) = \frac{100}{1000} \times (1 - 0.6449) = 0.0355$$

where 100 is the variance of a single survival time, 1,000 is the number of model evaluations and $-0.6449$ is the correlation within each pair. The standard error falls to 0.188 years.

4. Interpret. Ordinary sampling would need about 2,816 independent patients, 1,000 divided by 0.3551, to match the precision that antithetic sampling achieves with 1,000 evaluations, so for this output the pairing is worth roughly 2.8 times as many patients. The 64% cut in variance here is close to a best case, with one smooth monotone input; the 53% fewer replications reported for the diabetes model implies a variance factor of roughly 0.47 per evaluation, a smaller but still useful gain in a full patient-level model.

Where antithetic sampling gives little or makes things worse

The gain depends on the shape of the output, and some common health economic quantities weaken or reverse it. Two cases are common in health economic models.

The first case is a rare binary event. In a microsimulation an event with annual probability $p$ is typically simulated as occurring when $U < p$. For $p$ below 0.5 the two members of a pair cannot both have the event, which gives the following correlation.

$$\rho = -\frac{p}{1-p}$$

where $\rho$ is the correlation between the event indicators of the two members of a pair and $p$ is the event probability. For an illustrative annual probability of 0.02, $\rho$ is about $-0.020$ and the variance factor is about 0.980, a gain of only 2%. Models driven mainly by rare events gain little.

The second case is an output that is symmetric in its inputs. The expected value of perfect information depends on the expected maximum net benefit across strategies. In an illustrative extreme, suppose one input raises the net benefit of strategy A by exactly the amount it lowers that of strategy B, the two strategies have equal net benefit at the centre of that input's symmetric distribution, and that input is the only random one. The maximum is then V-shaped about the centre, both members of a pair give the same maximum and $\rho$ equals 1, so the antithetic estimate has twice the variance of ordinary sampling at the same cost. A near-tie between strategies is where the expected value of perfect information is largest and matters most to a decision, so the risk is real, and outputs built from maxima, absolute values or squared deviations deserve a pilot check.

Two reporting errors are common. Treating paired runs as independent when computing the standard error gives the wrong answer. It is too large when pairing helps, as in the worked example (0.316 rather than 0.188 years), and falsely reassuring when pairing hurts, so the standard error should be computed from the pair means. Averaging ICERs is wrong regardless of sampling method, because the mean of ratios is not the ratio of means; expected costs, expected QALYs and net benefit are the quantities to estimate.

Antithetic variates, common random numbers and other variance reduction methods

Antithetic variates and common random numbers both manipulate random numbers to control correlation, but they aim at opposite signs. Antithetic sampling induces negative correlation within the estimate of one expected output. Common random numbers give competing strategies the same random numbers so that their outputs are positively correlated, which reduces the variance of the difference between them because the covariance is subtracted. Technical Support Document 15 describes common random numbers as a popular variance reduction technique for comparing treatments in discrete event simulation, with the same patient characteristics sampled under intervention and control.

Since economic evaluation turns on incremental net benefit, common random numbers usually address the quantity of interest more directly than antithetic pairing. The two can be combined by running every strategy on the same original stream and on the same mirrored stream, so the strategies stay matched within each run while the pairs offset each other. Whether the combination helps the incremental estimate has to be checked, because incremental outputs can be less monotone than the outputs of each strategy.

Other variance reduction methods include control variates, Latin hypercube sampling and other forms of stratified sampling, and importance sampling. Antithetic pairing is among the simplest to add to an existing model, since it needs only control over the random numbers. None of these methods reduces parameter uncertainty or structural uncertainty, since they make a finite simulation a more precise estimate of what the model implies and leave the evidence behind the model unchanged.

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. Sections 3.6, 4.1 and 5.2. https://www.sheffield.ac.uk/nice-dsu/tsds/patient-level-simulation
• Hammersley JM, Morton KW. A new Monte Carlo technique: antithetic variates. Mathematical Proceedings of the Cambridge Philosophical Society. 1956;52(3):449-475. https://doi.org/10.1017/S0305004100031455
• McEwan P, Bergenheim K, Yuan Y, Tetlow AP, Gordon JP. Assessing the relationship between computational speed and precision: a case study comparing an interpreted versus compiled programming language using a stochastic simulation model in diabetes care. PharmacoEconomics. 2010;28(8):665-674. https://doi.org/10.2165/11535350-000000000-00000
• 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.2 and 8.6. https://artowen.su.domains/mc/Ch-var-basic.pdf