Response Surface Methodology

Response Surface Methodology in Health Economic Modelling

Response surface methodology, usually shortened to RSM, treats a slow decision model as an experiment that is run at a planned set of input combinations, and then fits a low-order polynomial to the outputs. The polynomial then stands in for the model, for example to find the screening interval with the highest expected net benefit or to run analyses that need too many evaluations. In health economics the method is one of several metamodelling techniques, and one of the more interpretable ones, because each term describes a slope, a curvature or an interaction between two inputs. Its weakness is that a quadratic is a local approximation, so the answer is only as good as the region the design covered.

The second-order model and coded inputs

The workhorse of the method is a full quadratic in the inputs. Degeling and colleagues, writing for health economists, note that the method is linear in the regression parameters but does not assume a linear input-output relationship.

$$y = \beta_0 + \sum_{i=1}^{k} \beta_i x_i + \sum_{i=1}^{k} \beta_{ii} x_i^2 + \sum_{i<j} \beta_{ij} x_i x_j + \varepsilon$$

where $y$ is the model output, such as incremental net monetary benefit, $x_1, \ldots, x_k$ are the $k$ coded inputs, $\beta_0$ is the intercept, $\beta_i$ are the first-order (slope) coefficients, $\beta_{ii}$ are the pure quadratic (curvature) coefficients, $\beta_{ij}$ are the two-factor interaction coefficients, and $\varepsilon$ is the error left by the approximation and by simulation noise.

Inputs are coded so that the low and high levels of each factor become minus 1 and plus 1. The NIST/SEMATECH e-Handbook of Statistical Methods recommends this because the estimates then do not depend on units of measurement.

$$x_i = \frac{X_i - (X_{\text{low}} + X_{\text{high}})/2}{(X_{\text{high}} - X_{\text{low}})/2}$$

where $X_i$ is input $i$ in its natural units, such as years between screens, and $X_{\text{low}}$ and $X_{\text{high}}$ are the low and high factorial levels of that input, which star points can exceed.

The coefficients are estimated by least squares, which Barton's tutorial on simulation metamodels writes in matrix form.

$$\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}$$

where $\hat{\boldsymbol{\beta}}$ is the vector of estimated coefficients, $\mathbf{X}$ is the design matrix with a column of ones, a column for each input and further columns for the squares and cross-products, and $\mathbf{y}$ is the vector of model outputs at the design points.

Designs that make a quadratic estimable

A quadratic needs at least three levels of each input, and a full grid at three levels grows quickly. The NIST handbook notes that a four-factor quadratic has 15 coefficients while the full three-level factorial needs 81 runs. Response surface designs estimate the quadratic with far fewer.

The central composite design starts from a two-level factorial design with centre points and adds two star (axial) points per input, which allow curvature to be estimated. In the rotatable version the star points sit at a coded distance $\alpha = (2^k)^{1/4}$ from the centre, 1.414 for two factors. The face-centred version sets $\alpha = 1$, needs only three levels and stays inside the original ranges, at the cost of poor precision for the pure quadratic coefficients. The Box-Behnken design also uses three levels and needs 15 runs for three factors against 20, an advantage the handbook says disappears at four or more.

Guidance for health economic metamodelling favours space-filling designs. Degeling and colleagues judge Latin hypercube designs preferable when few runs are affordable, and suggest starting with ten runs per input. A central composite design targets a local quadratic near a suspected optimum, while a space-filling design covers the whole range a later analysis will use. Both apply experimental design to a computer model, where the analyst controls every input.

The two-phase search for an optimum

The NIST handbook sets out a two-phase sequence it attributes to Box and Wilson. Far from the optimum, a first-order model fitted to a small factorial gives a direction of improvement, and the analyst keeps stepping along it until the output stops improving.

$$x_i^{*} = \rho \, \frac{b_i}{\sqrt{\sum_{j=1}^{k} b_j^2}}$$

where $x_i^{*}$ is the coded setting of input $i$ on the path of steepest ascent, $\rho$ is the chosen distance from the current centre, and $b_i$ is the fitted first-order coefficient for input $i$. When centre runs show significant curvature or lack of fit, axial runs are added and the full quadratic is fitted.

The fitted quadratic is then written in matrix form, and its stationary point is found by setting the partial derivatives to zero.

$$\hat{y}(\mathbf{x}) = b_0 + \mathbf{b}'\mathbf{x} + \mathbf{x}'\mathbf{B}\mathbf{x}, \qquad \mathbf{x}^{*} = -\tfrac{1}{2}\mathbf{B}^{-1}\mathbf{b}$$

where $\hat{y}$ is the predicted output, $b_0$ is the fitted intercept, $\mathbf{b}$ is the vector of first-order coefficients, $\mathbf{B}$ is the symmetric matrix with the pure quadratic coefficients on the diagonal and half of each interaction coefficient off the diagonal, and $\mathbf{x}^{*}$ is the stationary point in coded units. Substituting back gives $\hat{y}(\mathbf{x}^{*}) = b_0 + \tfrac{1}{2}\mathbf{b}'\mathbf{x}^{*}$.

Canonical analysis classifies $\mathbf{x}^{*}$. If all eigenvalues of $\mathbf{B}$ are negative the point is a maximum, if all are positive a minimum, and mixed signs mean a saddle. An eigenvalue near zero means a ridge of near-optimal settings, which in a policy model gives freedom to pick a practical option. If the optimum falls outside the experimental region, the handbook points to ridge analysis, a form of constrained optimisation, to keep the search inside it.

Response surfaces with decision models

In health technology assessment the experiment is usually a run of a slow simulation model, such as a discrete event simulation or microsimulation. Degeling and colleagues list value of information analysis, model calibration, optimisation and probabilistic analysis among earlier health economic uses of metamodels, and note that Tappenden and colleagues named the method among five techniques for value of information work.

Modelling choices shape the surface. The design ranges should match the ranges the analysis will explore, so in the example of Degeling and colleagues a screening interval feasible between 1 and 10 years spans that range in the training runs. For optimisation they prefer outputs from a probabilistic analysis at each design point, because those are the values expected in reality given current information. One metamodel can approximate net health benefit at a given willingness to pay, while costs and QALYs need one each.

Jalal and colleagues show the interpretive side. Regressing probabilistic sensitivity analysis outputs on standardised inputs makes the intercept the base-case outcome and lets the coefficients rank parameters, and an interaction term relates to two-way sensitivity analysis. They add that higher-order terms complicate reading the coefficients.

Worked example optimising a screening interval and starting age

All numbers in this example are illustrative. A cancer screening microsimulation reports mean incremental net monetary benefit (INMB) per person invited, against no screening, from a probabilistic analysis at each design point. The value of a QALY is GBP 25,000, the lower of the two values at which the NICE manual (PMG36, section 4.2.16) asks for net health benefits to be presented. The inputs are the screening interval, 1 to 5 years, and the starting age, 50 to 70 years.

1. Code the inputs. With centres of 3 and 60 years and half-ranges of 2 and 10 years, $x_1 = (X_1 - 3)/2$ and $x_2 = (X_2 - 60)/10$.

2. Run a face-centred central composite design. Nine runs (four factorial, four axial, one centre) keep every setting inside the feasible ranges.

3. Estimate the slopes and the interaction. The first-order and interaction columns are orthogonal to every other column in this design, so each of these coefficients is a simple contrast of runs at coded values of plus and minus 1. For the interval, $b_1 = (-43 + 51 - 23 + 47 - 42 + 58)/6 = 48/6 = 8$. For the starting age, $b_2 = (-43 - 51 + 23 + 47 - 59 + 53)/6 = -30/6 = -5$. For the interaction, $b_{12} = (43 - 51 - 23 + 47)/4 = 16/4 = 4$.

4. Fit the full quadratic. Least squares on all nine runs gives $\hat{y} = 60 + 8x_1 - 5x_2 - 12x_1^2 - 6x_2^2 + 4x_1x_2$. The residual sum of squares is 48 on 3 degrees of freedom, a residual standard deviation of GBP 4 per person, and $R^2 = 0.95$.

5. Locate the stationary point. Here $\mathbf{b} = (8, -5)'$ and $\mathbf{B}$ has diagonal elements of minus 12 and minus 6 and off-diagonal elements of 2. Its determinant is 72 minus 4, or 68, and $\mathbf{x}^{*} = -\tfrac{1}{2}\mathbf{B}^{-1}\mathbf{b} = (19/68, -22/68) = (0.279, -0.324)$.

6. Classify the point. The eigenvalues of $\mathbf{B}$ solve $\lambda^2 + 18\lambda + 68 = 0$, giving $-9 \pm \sqrt{13}$, or minus 5.39 and minus 12.61. Both are negative, so the point is a maximum.

7. Convert and predict. In natural units the optimum is $X_1 = 3 + 2 \times 0.279 = 3.56$ years and $X_2 = 60 + 10 \times (-0.324) = 56.8$ years. The predicted value is $\hat{y}(\mathbf{x}^{*}) = 60 + 0.5 \times (8 \times 0.279 + (-5) \times (-0.324)) = 61.93$, about GBP 62 per person invited.

8. Compare practical options. At a starting age of 57, a 3-year interval predicts GBP 60.96 per person and a 4-year interval GBP 61.36. The gap of GBP 0.40, or GBP 40,000 per 100,000 invited, is a tenth of the residual standard deviation.

9. Interpretation. The best region is an interval of 3 to 4 years starting at about 57, and predicted differences there are small relative to the residual standard deviation. The surface cannot rank the two intervals, so both should be rerun in the original microsimulation with enough replications to separate them, and the fit repeated at GBP 35,000 per QALY.

What a fitted optimum can support in an evaluation

A response surface optimum is an estimate built on an estimate. The NIST handbook stresses that a repeated experiment would give a different surface and stationary point, and Degeling and colleagues recommend re-evaluating several of the best strategies in the original model to check their outcomes and ranking. They caution that a strategy found optimal may lie where the simulator's data give little support, and the result should then only start a discussion about the simulator's validity.

The output also carries a value judgement. A surface fitted to net benefit at GBP 25,000 per QALY can move at GBP 35,000, the other value PMG36 names. Separate surfaces for incremental costs and QALYs let net benefit be rebuilt at any value without new model runs. Capacity limits, such as the colonoscopies a service can provide, make the search a constrained problem, and multi-objective optimisation applies when health gain and cost are kept apart.

Response surfaces beside other metamodels and search methods

Response surface methodology is the polynomial member of a wider family. Degeling and colleagues describe Gaussian process regression as useful when few runs can be afforded and unique in giving prediction uncertainty directly, though slow as runs and inputs grow. Bayesian optimisation uses such a Gaussian process surrogate to choose each new model run, where classical response surface methodology fixes each block of runs in advance. Calibration fits model inputs to observed data and can use a metamodel to cut run time. A deterministic sensitivity analysis traces the output along one or two inputs at a time, while a response surface describes the joint effect of several inputs.

Limitations and common misreadings

Most problems arise when a local quadratic is asked to describe a whole model. Health economic simulators can contain sharp edges that polynomials handle badly.

1. Treating the quadratic as global. The NIST handbook calls fitted responses local approximations. Barton notes that high-order polynomials fit by cancelling large oscillations, which makes them unstable, and points to other bases such as splines.
2. Extrapolating beyond the design. Degeling and colleagues advise against predicting outside the fitted ranges.
3. Ignoring clinical decision rules. Degeling and colleagues note that rigid cut-offs in a simulator can be hard for metamodels to capture.
4. Fitting to noisy runs. Too few simulated patients per design point add noise, and Degeling and colleagues suggest variance reduction such as common random numbers.
5. Reading the surface as better than the model. Degeling and colleagues note that a metamodel can be at most as accurate as its simulator and will not correct errors in it, so model validation comes first.

Sources

• Barton RR. Simulation metamodels. In Medeiros DJ, Watson EF, Carson JS, Manivannan MS, eds. Proceedings of the 1998 Winter Simulation Conference. IEEE; 1998:167-174. https://doi.org/10.1109/WSC.1998.744912
• Degeling K, IJzerman MJ, Lavieri MS, Strong M, Koffijberg H. Introduction to metamodeling for reducing computational burden of advanced analyses with health economic models: a structured overview of metamodeling methods in a 6-step application process. Medical Decision Making. 2020;40(3):348-363. https://doi.org/10.1177/0272989X20912233
• Jalal H, Dowd B, Sainfort F, Kuntz KM. Linear regression metamodeling as a tool to summarize and present simulation model results. Medical Decision Making. 2013;33(7):880-890. https://doi.org/10.1177/0272989X13492014
• 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. Section 4.2.16. https://www.nice.org.uk/process/pmg36
• NIST/SEMATECH. e-Handbook of Statistical Methods. NIST Handbook 151. National Institute of Standards and Technology; 2012. Sections 5.3.3.6 (response surface designs) and 5.5.3 (process optimisation). https://doi.org/10.18434/M32189