Bayesian Optimisation

Bayesian Optimisation for Calibrating Expensive Disease Models

Many health economic models contain natural history parameters that no study measures directly, so analysts adjust them until the model reproduces observed incidence, prevalence or mortality. When one microsimulation run takes minutes or hours, that search becomes the slowest step in building the model. Bayesian optimisation fits a Gaussian process to the runs already made and uses an acquisition function to decide where the model should be run next. Applied to model calibration, it can reach a good fit in few runs, as a small illustrative calibration and uses in cancer and malaria models show, and it differs in what it returns from other search and calibration methods.

Why calibrating simulation models calls for a sparing search

Calibration compares model outputs with empirical data and keeps the parameter values that fit well. Vanni and colleagues divide it into seven steps, of which the fourth is the search strategy and the sixth the stopping rule. The search strategy matters most when runs are expensive, because the number of runs then sets the cost.

Published calibrations have used manual adjustment, random search, grid search, the Nelder-Mead method, simulated annealing and genetic algorithms. In a Markov model of human papillomavirus and cervical disease, Taylor and colleagues found mean deviations from the targets of 7% for Nelder-Mead, 10% for manual calibration and 39% for random search. Kong and colleagues fitted 28 parameters of the Lung Cancer Policy Model to 11 targets, where simulated annealing beat a genetic algorithm. None of these builds a model of the fit surface from all earlier runs, and grid search uses none of them. A grid with 11 values for each of five parameters needs $N = 11^5 = 161{,}051$ runs, out of reach for a microsimulation that takes minutes per run.

The calibration problem as an optimisation

Bayesian optimisation needs one objective to minimise. The usual form folds the distances between model outputs and targets into a weighted goodness-of-fit score. The search then runs over a bounded region.

$$\theta^{*} = \arg\min_{\theta \in \Theta} L(\theta), \qquad L(\theta) = \sum_{j=1}^{J} w_j \, d\big(g_j(\theta), T_j\big)$$

where $\theta$ is the vector of parameters being calibrated, $\Theta$ is the region searched, $L(\theta)$ is the total loss, $J$ is the number of targets, $T_j$ is the observed value of target $j$, $g_j(\theta)$ is the model's prediction of it, $d$ is a distance such as an absolute or squared difference, and $w_j$ is the weight on target $j$. Each evaluation of $L(\theta)$ is a full model run with no formula or derivatives, the black-box, derivative-free setting for which Frazier describes the method as designed.

The surrogate, a Gaussian process over the fit surface

The first component is a statistical model of the unknown loss. A Gaussian process treats the losses at any set of points as jointly normal, with a kernel that makes nearby parameter values have similar losses. A common choice is the squared exponential kernel, shown below.

$$k(\theta, \theta') = \sigma_f^{2} \exp\left( -\frac{\lVert \theta - \theta' \rVert^{2}}{2\ell^{2}} \right)$$

where $k(\theta, \theta')$ is the prior covariance between the losses at $\theta$ and $\theta'$, $\sigma_f^{2}$ is the prior variance of the loss, $\lVert \theta - \theta' \rVert$ is the distance between the two parameter vectors, and $\ell$ is the length-scale, which sets how quickly the loss is expected to change. After $n$ runs, Gaussian process regression gives a normal posterior for the loss at any untried point, as follows.

$$\mu_n(\theta) = \mu_0 + \mathbf{k}_n(\theta)^{\top} \left[ K_n + \sigma_{\varepsilon}^{2} I \right]^{-1} (\mathbf{y} - \mu_0 \mathbf{1}), \qquad \sigma_n^{2}(\theta) = k(\theta, \theta) - \mathbf{k}_n(\theta)^{\top} \left[ K_n + \sigma_{\varepsilon}^{2} I \right]^{-1} \mathbf{k}_n(\theta)$$

where $\mu_n(\theta)$ and $\sigma_n^{2}(\theta)$ are the posterior mean and variance of the loss at $\theta$, $\mu_0$ is a constant prior mean, $\mathbf{k}_n(\theta)$ holds the covariances between $\theta$ and the evaluated points, $K_n$ is the covariance matrix of the evaluated points, $\sigma_{\varepsilon}^{2}$ is the simulation noise variance, $I$ is the identity matrix, $\mathbf{y}$ is the vector of observed losses and $\mathbf{1}$ is a vector of ones. The noise term is zero for a deterministic cohort model but not for a microsimulation with a finite number of simulated people.

The acquisition function and the choice of the next model run

The second component turns the posterior into a decision. Expected improvement, which Frazier traces to Močkus and to its popularisation by Jones and colleagues, values a candidate by how far it is expected to beat the best loss so far. For a loss that is minimised it takes the following form.

$$EI_n(\theta) = \big(L_n^{*} - \mu_n(\theta)\big)\,\Phi(z) + \sigma_n(\theta)\,\phi(z), \qquad z = \frac{L_n^{*} - \mu_n(\theta)}{\sigma_n(\theta)}$$

where $EI_n(\theta)$ is the expected improvement at $\theta$ after $n$ runs, $L_n^{*}$ is the lowest loss observed so far, $\sigma_n(\theta)$ is the posterior standard deviation, and $\Phi$ and $\phi$ are the standard normal distribution and density functions. The first term rewards points predicted to fit well and the second rewards points where the surrogate is uncertain. With a noisy microsimulation the best observed loss is itself uncertain, so implementations typically use the lowest posterior mean at the evaluated points in its place. The knowledge gradient and entropy search are alternatives, and Frazier reports that the knowledge gradient can substantially outperform expected improvement when evaluations are noisy. The procedure starts from a small space-filling design, such as random points or a Latin hypercube, then repeats three steps until the run budget is spent. Each round fits the surrogate, maximises the acquisition function (cheap, since it uses only the surrogate), and runs the model at the chosen point.

Worked example fitting one transition probability to a ten-year incidence target

All numbers here are illustrative. A cohort model has one unknown annual probability $p$ of moving from a healthy to a disease state, with no deaths, so it predicts ten-year cumulative incidence $I_{10} = 1 - (1 - p)^{10}$. The target is 30% and the loss is the absolute difference in percentage points. The exact answer, $p = 1 - 0.7^{0.1} \approx 0.0350$, is known, so the search can be checked. The search runs from $p = 0.01$ to $p = 0.06$, and the surrogate has $\mu_0 = 20$, $\sigma_f^{2} = 100$, $\ell = 0.02$ and no noise, all fixed for clarity.

1. Initial runs. At $p = 0.01$ the model predicts $I_{10} = 1 - 0.99^{10} = 0.0956$, or 9.56%, a loss of 20.44 points. At $p = 0.06$ it predicts $I_{10} = 1 - 0.94^{10} = 0.4614$, or 46.14%, a loss of 16.14 points, the best so far.

2. Surrogate at the candidate 0.04. Writing each covariance as 100 times a correlation, the two runs have correlation $\exp(-0.05^2 / 0.0008) = 0.0439$, so the determinant of the two-run correlation matrix is $D = 1 - 0.0439^2 = 0.99807$. The candidate's correlations with the two runs are $\exp(-1.125) = 0.3247$ and $\exp(-0.5) = 0.6065$. The weights on the two observed losses are as follows.

$$a_1 = \frac{0.3247 - 0.0439 \times 0.6065}{0.99807} = \frac{0.29807}{0.99807} = 0.2986, \qquad a_2 = \frac{0.6065 - 0.0439 \times 0.3247}{0.99807} = \frac{0.59225}{0.99807} = 0.5934$$

where $a_1$ and $a_2$ apply to the deviations of the first and second losses from the prior mean. The posterior mean and variance then follow.

$$\mu_2 = 20 + 0.2986 \times 0.44 + 0.5934 \times (-3.86) = 20 + 0.1314 - 2.2905 = 17.84$$

where 0.44 and $-3.86$ are the two losses minus the prior mean of 20, and $\mu_2$ is the predicted loss at $p = 0.04$.

$$\sigma_2^{2} = 100 \times (1 - (0.3247 \times 0.2986 + 0.6065 \times 0.5934)) = 100 \times (1 - 0.45686) = 54.31$$

where $\sigma_2^{2}$ is the posterior variance, giving a standard deviation of $\sqrt{54.31} \approx 7.37$.

3. Expected improvement. The standardised gap and the expected improvement are as follows.

$$z = \frac{16.14 - 17.84}{7.37} = -0.23, \qquad EI = -1.70 \times 0.4090 + 7.37 \times 0.3885 = -0.6953 + 2.8632 = 2.17$$

where 0.4090 and 0.3885 are $\Phi(-0.23)$ and $\phi(-0.23)$. The predicted loss is worse than the best so far, but the uncertainty makes the expected improvement positive.

4. Choosing the next run. Repeating the calculation for every candidate from 0.020 to 0.050, 0.005 apart, gives the values in the table.

The candidates at 0.015 and 0.055 score lower still. The highest value is at 0.040, where the model predicts $I_{10} = 1 - 0.96^{10} = 0.3352$, a loss of 3.52 points.

5. Second iteration. Refitted to three runs, the surrogate gives the highest expected improvement at $p = 0.035$ (0.46, against 0.29 at 0.030). The model predicts $I_{10} = 1 - 0.965^{10} = 0.2997$, a loss of 0.03 points, at a value of $p$ within 0.0001 of the exact answer. Four runs found what a grid at the same spacing would need 11 runs to find, and a full grid grows multiplicatively with every added parameter. The search itself cannot tell that it has reached the exact answer, so in practice a stopping rule or run budget ends it.

Design choices that change the result

Several settings decide where the search ends. Gómez-Guillén and colleagues searched 50% either side of literature-based starting values, so a true value outside that range could not be found. The kernel encodes how smooth the loss is assumed to be; the absolute-difference loss above has a sharp minimum, and one more iteration of the example pushed the posterior mean below zero, a loss that cannot occur.

Frazier describes the usual setting as 20 or fewer inputs and a budget limited to a few hundred evaluations. The weights $w_j$ decide which targets dominate when they conflict.

Evidence from cancer and malaria models

Gómez-Guillén and colleagues took a lung cancer Markov microsimulation from a published cost-effectiveness analysis and compared Bayesian optimisation, using expected improvement, with Nelder-Mead, simulated annealing and particle swarm optimisation. For all 99 parameters at once, standard Bayesian optimisation was stopped at 1,000 evaluations without a satisfactory fit, although two variants built for many parameters did better. For one age group with 11 parameters it was the fastest method once a run took longer than 0.2 seconds, and calibrating age group by age group cut the projected break-even run time for the full model from about 5 minutes to 0.24 seconds.

Reiker and colleagues calibrated 23 parameters of OpenMalaria, an individual-based malaria model, against 11 epidemiological objectives; malaria models of this kind have been used for economic analyses of new interventions. Bayesian optimisation with Gaussian process or machine learning emulators quickly beat the fit of the genetic algorithm used before.

What a calibrated point can and cannot support in an evaluation

Bayesian optimisation returns a single best-fitting parameter set. Menzies and colleagues note that such a set can be a useful first step, but that estimates of uncertainty around modelled results need a distribution of parameter values. A probabilistic sensitivity analysis therefore needs Bayesian calibration or a similar method, which combines priors, model structure and a likelihood into a posterior.

A single optimum can also hide nonidentifiability. Alarid-Escudero and colleagues showed that two equally good-fitting parameter sets for a cancer survival model gave life expectancy gains for a hypothetical treatment of 0.67 and 0.31 years. An optimiser that stops at one of them gives no warning that the other exists, so calibration reports should check identifiability and whether near-optimal sets change the decision.

Bayesian optimisation beside grid search, Nelder-Mead and approximate Bayesian computation

The methods differ in what they return and in how they spend runs. Approximate Bayesian computation answers a different question, since it targets a posterior rather than a best point.

Shewmaker and colleagues set out rejection and Markov chain Monte Carlo versions of approximate Bayesian computation for microsimulation.

Where Bayesian optimisation is misread or misapplied

The word Bayesian can mislead. The Bayesian updating applies to the surrogate's beliefs about the loss surface. The model parameters receive no prior or posterior, so the output cannot be read as a Bayesian analysis of them. It also differs from the Bayesian calibration of Kennedy and O'Hagan, whose predictions allow for remaining uncertainty over the fitted parameters and for model inadequacy, and from Bayesian adaptive designs, which update a clinical trial as data accrue.

The method pays only when runs are slow. For a fast model, fitting the surrogate can cost more than many extra runs, as the lung cancer comparison showed. Its result is the best point found within a given budget, region and kernel, with no guarantee that it is the global optimum, so reports should state these with the acquisition function and initial design.

Sources

• Alarid-Escudero F, MacLehose RF, Peralta Y, Kuntz KM, Enns EA. Nonidentifiability in model calibration and implications for medical decision making. Medical Decision Making. 2018;38(7):810-821. doi:10.1177/0272989X18792283
• Frazier PI. A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.02811 [stat.ML], version 1, 8 July 2018. https://arxiv.org/abs/1807.02811
• Gómez-Guillén D, Díaz M, Arcos JL, Cerquides J. Bayesian optimization with additive kernels for a stepwise calibration of simulation models for cost-effectiveness analysis. International Journal of Computational Intelligence Systems. 2024;17:249. doi:10.1007/s44196-024-00646-x
• Kennedy MC, O'Hagan A. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2001;63(3):425-464. doi:10.1111/1467-9868.00294
• Kong CY, McMahon PM, Gazelle GS. Calibration of disease simulation model using an engineering approach. Value in Health. 2009;12(4):521-529. doi:10.1111/j.1524-4733.2008.00484.x
• Menzies NA, Soeteman DI, Pandya A, Kim JJ. Bayesian methods for calibrating health policy models: a tutorial. PharmacoEconomics. 2017;35(6):613-624. doi:10.1007/s40273-017-0494-4
• Rasmussen CE, Williams CKI. Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press; 2006. Chapter 2.
• Reiker T, Golumbeanu M, Shattock A, Burgert L, Smith TA, Filippi S, Cameron E, Penny MA. Emulator-based Bayesian optimization for efficient multi-objective calibration of an individual-based model of malaria. Nature Communications. 2021;12:7212. doi:10.1038/s41467-021-27486-z
• Shewmaker P, Chrysanthopoulou SA, Iskandar R, Lake D, Jutkowitz E. Microsimulation model calibration with approximate Bayesian computation in R: a tutorial. Medical Decision Making. 2022;42(5):557-570. doi:10.1177/0272989X221085569
• Taylor DC, Pawar V, Kruzikas D, Gilmore KE, Pandya A, Iskandar R, Weinstein MC. Methods of model calibration: observations from a mathematical model of cervical cancer. PharmacoEconomics. 2010;28(11):995-1000. doi:10.2165/11538660-000000000-00000
• Vanni T, Karnon J, Madan J, White RG, Edmunds WJ, Foss AM, Legood R. Calibrating models in economic evaluation: a seven-step approach. PharmacoEconomics. 2011;29(1):35-49. doi:10.2165/11584600-000000000-00000