Gaussian Process

Gaussian Process in Health Economic Modelling

A Gaussian process is a way of putting a probability model on an unknown function itself, rather than on a handful of parameters. In health economics the unknown function is usually the link between a decision model's uncertain inputs and an output such as net monetary benefit, which can be slow to compute for each set of inputs. Once a few runs have been observed, the process gives a best guess and an uncertainty for the output at inputs that were never run, so a probabilistic sensitivity analysis or value of information calculation becomes feasible on a slow model. Its predictions follow from a mean function and a covariance function, and the covariance function carries most of the assumptions.

A distribution over functions built from finite pieces

Rasmussen and Williams define a Gaussian process as a collection of random variables, any finite number of which have a joint Gaussian distribution. Each random variable is the value of the function at one input, so a decision model evaluated at ten parameter sets corresponds to ten jointly normal outputs. The pieces must agree, a consistency or marginalisation property, so looking at a larger set of inputs never changes the distribution already assigned to a smaller set.

$$m(x) = E[f(x)], \qquad k(x, x') = E\big[(f(x) - m(x))(f(x') - m(x'))\big], \qquad f(x) \sim \mathcal{GP}\big(m(x), k(x, x')\big)$$

where $f(x)$ is the unknown function at input $x$, $m(x)$ is the mean function, $k(x, x')$ is the covariance function between the outputs at inputs $x$ and $x'$, and $\mathcal{GP}$ denotes a Gaussian process.

For any finite set of inputs the process is an ordinary multivariate normal distribution, with $m(x)$ in the mean vector and $k(x_i, x_j)$ in the variance-covariance matrix. The covariance function must be positive semidefinite so that this matrix is valid for every set of inputs. Oakley and O'Hagan treat a complex computer model as a deterministic function $y = \eta(x)$ whose response to its inputs is not transparent, with uncertainty about $\eta(\cdot)$ represented by a Gaussian process prior.

The covariance function sets smoothness and scale

The covariance function says how much the output at one input reveals about the output at another. The squared exponential function is the running example in Rasmussen and Williams.

$$k_{SE}(x, x') = \sigma_f^2 \exp\left(-\frac{(x - x')^2}{2\ell^2}\right)$$

where $\sigma_f^2$ is the signal variance, the prior variance of the output at any single input, and $\ell$ is the characteristic length-scale. Rasmussen and Williams describe the length-scale as roughly the distance one has to move in input space before the function value can change significantly.

This kernel is infinitely differentiable, so the functions it generates are very smooth. Rasmussen and Williams report Stein's argument that such strong smoothness is unrealistic for many physical processes and his recommendation of the Matérn class, in which a parameter $\nu$ controls how rough the function may be. Decision models have many inputs, and a separate length-scale for each is the usual extension, which Rasmussen and Williams call automatic relevance determination because a very large length-scale marks an input as close to irrelevant. Strong, Oakley and Brennan use the same idea for EVPPI, with a correlation length $\delta_j$ for each parameter of interest describing how smoothly conditional net benefit changes with it.

The mean function is the second choice. Rasmussen and Williams usually take it as zero, while Oakley and O'Hagan use regressor functions $h(x)$ for beliefs about the general way the output responds, setting $h(x)^T = (1, x^T)$ in one example, a constant plus a linear term in each input. The mean matters most away from the observed runs, where predictions drift back towards it.

Hyperparameters and the marginal likelihood

The signal variance, the length-scales and any noise variance are hyperparameters, and their values affect every prediction. Rasmussen and Williams describe setting them by maximising the marginal likelihood, the probability of the observed outputs with the function values integrated out, and show that it trades fit against complexity automatically, so it does not simply favour the most flexible model. Oakley and O'Hagan note that integrating over the hyperparameters by Markov chain Monte Carlo is highly intensive and that experience suggests plugging in estimates is adequate. Plug-in estimates leave uncertainty about the hyperparameters out of the reported prediction uncertainty.

Why a smooth prior saves model runs

The case for a Gaussian process rests on smoothness. Oakley and O'Hagan state the requirement as a belief that knowing $\eta(x)$ gives some idea of $\eta(x')$ for nearby $x'$, information that plain Monte Carlo simulation ignores after each run. Their design points are chosen to give good information about the function rather than drawn at random, starting in their examples from a Latin hypercube sampling design.

With the hyperparameters fixed, conditioning the prior on the runs gives a posterior that is again a Gaussian process, with a variance that shrinks near the runs, and the Gaussian process regression page gives the matrix formulas. Oakley and O'Hagan report that the standard deviations of their estimates were generally much smaller, often by orders of magnitude, than those from Monte Carlo methods with the same number of runs. They also state that the approach is not appropriate if the model is known to respond erratically or discontinuously to its inputs.

Two uses in economic evaluation

Gaussian processes appear in health economic work in two different ways. The difference changes what the fitted function means and what it is fitted to.

The first is emulation of the whole model as a function of all its uncertain inputs. Stevenson, Oakley and Chilcott applied it to an individual patient simulation of treatments for established osteoporosis, with minimal loss of accuracy, and each further run, which had taken 150 minutes, became virtually instantaneous, so that probabilistic sensitivity analysis became possible. Rojnik and Naveršnik used Gaussian process and multiple linear regression metamodels to compute partial EVPI for six groups of parameters in a breast cancer screening model. The Gaussian process performed better, and with it the computing time was 47 days instead of 44 years.

The second use, set out by Strong, Oakley and Brennan, works on an existing probabilistic sample and models only the conditional expected net benefit as a function of the parameters of interest. The sampled net benefits are treated as noisy data about that function, and the fitted values give the inner expectation needed for EVPPI without a nested Monte Carlo loop. The authors call this subtly different from the earlier applications, which modelled net benefit as a function of all the inputs. They recommend the Gaussian process over a generalised additive model when the parameters of interest number more than 6 or so and are expected to interact.

Worked example with a squared exponential kernel (illustrative)

All numbers here are illustrative. An analyst places a Gaussian process prior on incremental net monetary benefit (INMB), in GBP thousands, as a function of one input $x$, the hazard ratio for progression on the new treatment. The prior mean is zero, $\sigma_f^2 = 4$ (a prior standard deviation of GBP 2,000) and $\ell = 0.1$. The inputs of interest are $x_1 = 0.70$, $x_2 = 0.75$ and $x_3 = 0.90$.

1. Prior covariances. For $x_1$ and $x_2$ the exponent is $(x_1 - x_2)^2 / (2\ell^2) = 0.0025 / 0.02 = 0.125$, and for $x_1$ and $x_3$ it is $(x_1 - x_3)^2 / (2\ell^2) = 0.04 / 0.02 = 2$.

$$k(x_1, x_2) = 4 \times e^{-0.125} = 4 \times 0.8825 = 3.530, \qquad k(x_1, x_3) = 4 \times e^{-2} = 4 \times 0.1353 = 0.541$$

where each $k$ is a prior covariance between INMB values, in squared GBP thousands. Dividing by the prior variance of 4 gives correlations of 0.88 for the nearby pair and 0.14 for the distant pair.

2. Condition on one model run. The model is run once at $x_1 = 0.70$ and returns an INMB of 3, or GBP 3,000, with no noise. For a single noise-free run the conditional normal formulas are as follows.

$$\mu(x_*) = m(x_*) + \frac{k(x_*, x_1)}{k(x_1, x_1)}\big(y_1 - m(x_1)\big), \qquad v(x_*) = k(x_*, x_*) - \frac{k(x_*, x_1)^2}{k(x_1, x_1)}$$

where $x_*$ is the prediction input, $y_1 = 3$ is the observed output, $m(\cdot) = 0$ is the prior mean, $k(x_1, x_1) = k(x_*, x_*) = 4$ is the prior variance, $\mu(x_*)$ is the predictive mean and $v(x_*)$ the predictive variance.

3. Prediction at the nearby input. At $x_2$ the mean is $\mu(x_2) = 0.8825 \times 3 = 2.65$ and the variance is $v(x_2) = 4 - 12.461 / 4 = 4 - 3.115 = 0.885$, a standard deviation of 0.94. With $\Phi$ the standard normal distribution function, the probability that INMB is positive is $\Phi(2.65 / 0.94) = \Phi(2.82) = 0.998$.

4. Prediction at the distant input. At $x_3$ the mean is $\mu(x_3) = 0.1353 \times 3 = 0.41$ and the variance is $v(x_3) = 4 - 0.293 / 4 = 4 - 0.073 = 3.927$, a standard deviation of 1.98, so the prediction has fallen most of the way back to the prior, and $\Phi(0.41 / 1.98) = \Phi(0.21) = 0.58$.

5. Change the length-scale. With $\ell = 0.2$ the exponent for $x_1$ and $x_3$ is 0.04 / 0.08 = 0.5 and the covariance is $k(x_1, x_3) = 4 \times e^{-0.5} = 4 \times 0.6065 = 2.426$. The prediction at $x_3$ becomes $\mu(x_3) = 0.6065 \times 3 = 1.82$ with $v(x_3) = 4 - 5.885 / 4 = 4 - 1.471 = 2.529$, a standard deviation of 1.59, and $\Phi(1.82 / 1.59) = \Phi(1.14) = 0.87$.

6. Interpretation. One run says a great deal about hazard ratios near 0.70 and little about 0.90, and the length-scale alone moves the probability of positive INMB at 0.90 between 0.58 and 0.87. That probability describes uncertainty about what the slow model would return. It says nothing about how likely a hazard ratio of 0.90 is, which is the job of the input distributions in the probabilistic analysis. A real emulator uses many runs across all the uncertain inputs and estimates the length-scales from them.

Assumptions an analyst has to defend

The smoothness in the kernel does most of the work. Stopping rules or changes of treatment line can produce kinks or jumps that a smooth kernel rounds off, so predictions should be checked against held-out runs of the original model.

The noise assumption must match the model. A deterministic cohort model returns the same output for the same inputs and can be interpolated exactly, whereas an individual patient simulation carries Monte Carlo error, and Stevenson, Oakley and Chilcott note that large numbers of patients must be simulated to reduce this first-order uncertainty. In the EVPPI method of Strong, Oakley and Brennan the noise term is central, because variation due to parameters outside the group of interest is treated as noise around the conditional expectation.

The word Gaussian refers to the joint distribution of function values under the prior. It does not mean that net benefit across a probabilistic sample is normally distributed, and the model's outputs can be skewed or bounded.

Gaussian process regression and neighbouring methods

The Gaussian process is the prior model over functions, and Gaussian process regression is the procedure that conditions it on input and output pairs to predict. A Bayesian analysis of one parameter starts from a prior distribution over a number, whereas a Gaussian process is a prior over a whole function, and conditioning gives a posterior distribution over functions.

Outside value of information work, Oakley and O'Hagan used the fitted process to compute the main effects and interactions behind variance-based sensitivity analysis from a single set of runs, a strand of uncertainty quantification.

Limits and common misreadings

Most problems arise from forgetting that the process is a model of the model, with assumptions of its own. Five errors recur, and each can change a reported result.

1. Ignoring computational cost. Rasmussen and Williams give $n^3/6$ operations for the Cholesky factorisation with $n$ training points, and Strong, Oakley and Brennan note that the cubic growth limits a probabilistic sample to the order of tens of thousands.
2. Reading small predictive variance as small parameter uncertainty. Predictive variance measures how well the emulator knows the slow model's output, not how well the evidence pins down the inputs.
3. Trusting extrapolation. Outside the design, predictions return towards the prior mean, so the design should cover the input distributions used in the probabilistic analysis.
4. Treating hyperparameters as known. Plug-in length-scales and variances make predictive intervals too narrow when they are poorly determined.
5. Assuming one kernel suits every model. Different kernels can give materially different predictions between runs, so the choice should be reported with its validation.

Sources

• Oakley JE, O'Hagan A. Probabilistic sensitivity analysis of complex models: a Bayesian approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2004;66(3):751-769. https://doi.org/10.1111/j.1467-9868.2004.05304.x
• Rasmussen CE, Williams CKI. Gaussian Processes for Machine Learning. MIT Press; 2006. Chapters 2, 4 and 5. https://gaussianprocess.org/gpml/
• Rojnik K, Naveršnik K. Gaussian process metamodeling in Bayesian value of information analysis: a case of the complex health economic model for breast cancer screening. Value in Health. 2008;11(2):240-250. https://doi.org/10.1111/j.1524-4733.2007.00244.x
• Stevenson MD, Oakley J, Chilcott JB. Gaussian process modeling in conjunction with individual patient simulation modeling: a case study describing the calculation of cost-effectiveness ratios for the treatment of established osteoporosis. Medical Decision Making. 2004;24(1):89-100. https://doi.org/10.1177/0272989X03261561
• Strong M, Oakley JE, Brennan A. Estimating multiparameter partial expected value of perfect information from a probabilistic sensitivity analysis sample: a nonparametric regression approach. Medical Decision Making. 2014;34(3):311-326. https://doi.org/10.1177/0272989X13505910