Response surface methodology

What a response surface is, and when to use it

Response surface methodology (RSM) is "a collection of mathematical and statistical techniques" for modeling a response that depends on several variables, where "the objective is to optimize this response" (Montgomery, 2020, p. 408). The expected response is a function of the factors, η = f(x₁, x₂), and its graph is the "response surface"; it is displayed in a contour plot – "lines of constant response" in the plane of two factors (pp. 408-409).

The difference from an ordinary analysis of variance: there the factors are categorical, and the question is whether the means differ. Here the factors are quantitative and continuous – lesson length, number of exercises, dose, temperature – and the question is where in their range the response is best, and in which direction it is worth moving. The true f is unknown, so it is approximated by a polynomial: first order (Equation 11.1) or second order (11.2; p. 409). This is only a local approximation – "it is unlikely that a polynomial model will be a reasonable approximation over the entire space", but "in a relatively small region they usually work quite well" (p. 409).

Montgomery's examples come from chemical engineering – reaction time and temperature against yield. The same logic fits any study in which the researcher sets the levels of continuous variables: the length and frequency of an intervention, dose, number of items, exposure time. The sample data file on the page – fictitious data on a practice lesson – is built this way.

Working in stages

"RSM is a sequential procedure" (Montgomery, 2020, p. 410). When you are far from the optimum, the surface has almost no curvature and a first-order model is enough; you climb with it toward the optimum. Near the optimum curvature appears, and a second-order model is fitted to locate it – like "climbing a hill" (p. 410). The calculator is built for both stages: a linear model with the path of steepest ascent, and a second-order model with a stationary point.

In the first stage a 2^k experiment (or a fraction of it) with a few center points is enough. The center points make three things possible: estimating the error, checking interactions, and checking curvature (p. 412). The curvature check has one degree of freedom: if the surface is planar, the mean of the factorial points ȳ_F equals the mean of the center points ȳ_C; the difference estimates Σβ_jj, and SS = n_Fn_C(ȳ_F − ȳ_C)²/(n_F + n_C) is tested against the pure error (p. 413). The page computes it when the data are exactly a 2^k with center points, and breaks down the residuals of the linear model into interactions, pure quadratic and pure error – as in Tables 11.2 and 11.5 (pp. 413, 415).

The path of steepest ascent

"The method of steepest ascent is a procedure for moving sequentially in the direction of the maximum increase in the response"; for minimization – steepest descent (Montgomery, 2020, p. 411). The direction is perpendicular to the contour lines of the linear model, and "the steps along the path are proportional to the regression coefficients" (p. 411). The algorithm (p. 415): choose a step size Δx_j in one factor – "usually the factor we know the most about, or the one with the largest absolute regression coefficient"; in every other factor Δx_i = b_i/(b_j/Δx_j); and convert to natural units. In Example 11.1 the step is 5 minutes in time and 2 degrees in temperature (p. 416).

"The actual step size is determined by the experimenter based on process knowledge" (p. 411). Experiments are run along the path "until no further increase in response is observed", and there a new linear model is fitted; arriving in the vicinity of the optimum "is usually indicated by lack of fit of a first-order model" (p. 411). In Example 11.1 the response increased up to the tenth step and decreased after it (pp. 413-414), and the new model around the tenth step showed significant curvature – "the first-order model is not an adequate approximation" (p. 414). On the page – choose a base factor, a step size and a number of steps, and the table gives the path in coded and natural units, with a prediction only within the range of the data.

The second-order model: stationary point and canonical analysis

In matrix notation, ŷ = b₀ + x'b + x'Bx (Equation 11.5): b – the vector of the linear coefficients, and B – a symmetric matrix with the quadratic coefficients on the diagonal and half the interaction coefficients off it (Montgomery, 2020, pp. 417-418). The stationary point – where all the derivatives are zero – is x_s = −½B⁻¹b, and the response there is ŷ_s = b₀ + ½x_s'b (Equations 11.7-11.8, p. 418). It "could represent a point of maximum response, a point of minimum response, or a saddle point" (p. 416).

To find out which it is, you move to the canonical form: shift the origin to the stationary point and rotate the axes to the principal axes of the surface, so that ŷ = ŷ_s + λ₁w₁² + … + λ_kw_k² (Equation 11.9), where the λ are the eigenvalues of B (p. 418). All negative – a maximum; all positive – a minimum; mixed signs – a saddle point; and "the surface is steepest in the w_i direction for which |λ_i| is the greatest" (p. 418). The new axes are w = M'(x − x_s), where M – the normalized eigenvectors (Equation 11.10, p. 423). In Montgomery's Example 11.2: x_s = (0.389, 0.306), that is, 87 minutes and 176.5 degrees, ŷ_s = 80.21, and both eigenvalues negative (−0.963, −1.414) – a maximum (p. 420); the page reproduces all these numbers.

Ridges. When one of the eigenvalues is close to zero, the response is "very insensitive" to movement along its axis (p. 424). If the stationary point is inside the region of the experiment – this is a "stationary ridge": a line of nearly equal points, and flexibility to choose among them by cost or convenience. If it is far outside the region – a rising (or falling) ridge: "we cannot draw inferences about the true surface or about the stationary point", but it is worth exploring further in the direction of the flat axis (pp. 424-425). The page flags an eigenvalue smaller than a tenth of the largest one in absolute value – the threshold of canonical in the rsm package (site convention) – and reports when the point is outside the range of the data.

Multiple responses. When there are several responses, a model is fitted for each and a combination that satisfies all of them is sought – by overlaying contour plots, by constrained optimization, or with a desirability function that combines all the responses into a single measure (pp. 425-429). The page analyzes one response per run; for several responses – run it for each one and compare the maps.

Which design?

For a first-order model – a 2^k or a fraction in which the main effects are not aliased with each other; center points do not change the linear coefficients or the orthogonality, and add an estimate of the error (Montgomery, 2020, p. 430). On the "Design of experiments" page – 2^k and fractions.

Central composite design (CCD) – "the most popular class" for a second-order model (p. 430): a 2^k (or a resolution V fraction) with n_F runs, 2k axial points at ±α on each axis, and n_C center points (p. 431). It lends itself to sequential construction: first run the 2^k with the center points, and if there is lack of fit – add the axial points (p. 431). Two decisions: α and the number of center points.

• Rotatable – α = n_F^1/4 (1.414 for k = 2, 1.682 for k = 3): the prediction variance is equal at all points at the same distance from the center – "the variance of predicted response is constant on spheres" (pp. 431-432).
• Spherical – α = √k: all the points on one sphere; for a spherical region this is the good choice in terms of prediction variance, and "exact rotatability is not necessary to have a good design" (p. 432). Three to five center points are recommended (p. 432).
• Face-centered – α = 1: only three levels, for a cuboidal region, when it is impossible to go beyond the range; it is not rotatable, and "two or three" center points are enough (p. 433).
• Orthogonal blocking – when the factorial and the axial parts run at different times, they are two blocks. The design "blocks orthogonally" when the block effect does not affect the coefficients (p. 437); in a CCD this happens when α = [n_F(n_A + n_CA)/(2(n_F + n_CF))]^1/2 (Equation 11.18, p. 439). It is possible to "sacrifice a little rotatability or orthogonal blocking without a great loss of information" (p. 439). Only center points that ran in the same block are replicates (p. 440).

Box-Behnken – a three-level design, "formed by combining 2^k factorials with incomplete block designs": every pair of factors at ±1 and the rest at the center, plus center points. It is "very efficient in terms of the number of required runs", "rotatable or nearly rotatable", all its points lie on a sphere of radius √2, and it has no corners – an advantage when the extreme combination of all the factors is expensive or impossible (p. 432; Table 11.8 for k = 3, p. 433). On the page – k = 3 to 5, with the same construction.

In the design, the page shows the scaled prediction variance N·x'(X'X)⁻¹x (Equation 11.16, p. 436) at the center, at a factorial point and at an axial point – in a rotatable design the last two are equal when they are at the same distance. Optimal designs (D, I), mixture designs and designs for computer experiments (pp. 440-472) – are not on the page.

Is the model adequate?

Lack of fit versus pure error. The residuals of the model break down into two: pure error – the variance between replicate runs at the same point, which does not depend on the model; and lack of fit – what remains. If the lack of fit is large relative to the pure error, the model is missing something (Montgomery, 2020, Tables 11.2, 11.5 and 11.7). This is why center points matter: without replicates there is no pure error and no such test.

Model comparison. Table 1 in the output compares the three models – linear, with interactions, and second order – like "Model Summary Statistics" and "Lack of Fit Tests" in Montgomery (Table 11.7, p. 421): the standard deviation of the residuals, R², adjusted R², predicted R² and PRESS, and the lack of fit. In Example 11.2 the linear model and the one with interactions fail (negative predicted R², significant lack of fit), and the second-order model fits (predicted R² .918, nonsignificant lack of fit, p = .290). Predicted R² relies on PRESS – the sum of the squared residuals when each run is predicted without itself.

Diagnostics. As in any regression – residuals, leverage and influence (Table 11.7, p. 422). On the page – the collapsed diagnostics table at the end of the output.

How to read the output?

In the design: Table 1 – the design and its properties: the parts, α, the number of runs and coefficients, the degrees of freedom for lack of fit and pure error, the blocking, and the prediction variance. Table 2 – the run sheet in random order, with the coded and natural values. Figure 1 – the design points. The Excel download includes a codebook sheet with the coding, and the page reads it back in the analysis.

In the analysis: before the run – the coding table: for each factor a center and a half-range (x = (ξ − center)/half-range), from the codebook sheet or from the values. Table 1 – the model comparison. Table 2 – the analysis of variance of the selected model: each term with a partial sum of squares, and residuals broken down into lack of fit and pure error. Table 3 – the coefficients in coded units, with confidence intervals and VIF, and below it the equation in original units. In a second-order model – Table 4, the stationary point, and Table 5, the canonical analysis; in a linear model – the path of steepest ascent. Figure – a contour map for each pair of factors, with the others at the center or at the stationary point. After it – prediction at a point with a confidence interval and a prediction interval, a diagnostics table, a results paragraph and the "Interpreting the result" box.

Assumptions and warnings

The approximation is local. The polynomial fits the region that was examined, not beyond it (Montgomery, 2020, p. 409). Any prediction outside the design points is extrapolation; the page does not compute predictions on the path outside the range, and warns when the stationary point is outside it. The analysis of the stationary point assumes that it is "within the region of exploration" (p. 418).

Randomization. As in any experiment, the run order is random – the run sheet the page builds. The assumptions about the error are those of ordinary regression: independent, equal variance, a distribution close to normal – checked in the residuals.

Coding. The coefficients, the stationary point and the canonical analysis – in coded units, as in Montgomery. The equation in original units is equivalent, but its coefficients depend on the units and should not be compared with each other (p. 422). The eigenvalues depend on the coding; a coding in which ±1 are the factorial levels preserves their meaning.

What this page does not compute. Optimal designs (D, I, G), mixture designs, computer experiments, EVOP (Montgomery, 2020, Sections 11.4.4-11.7), a desirability function for several responses, 3^k experiments (chapter 9) and robust designs (chapter 12). On the "Design of experiments" page – categorical factors, 2^k and fractions; and on the "Random, nested and split-plot factors" page – designs with a restriction on randomization.

And how do you do it in SPSS?

SPSS Statistics has no procedure for response surfaces. Compute the coded variables, their products and their squares with Transform › Compute Variable, and run Analyze › General Linear Model › Univariate with all of them as Covariates, with Type III sums of squares – these are the partial SS of Table 2 – and with Lack of fit in Options. The syntax the page generates does exactly this (COMPUTE and then UNIANOVA … WITH … /PRINT=PARAMETER LOF). SPSS does not compute the stationary point, the canonical analysis and the maps – for them use R or the page. Without blocks, the syntax also adds REGRESSION with the same terms – the same coefficients. The menu paths are general information about the software, and the syntax was not run in SPSS. PSPP, the free software that reads SPSS syntax, has no UNIANOVA; we ran the REGRESSION part in it on the sample data file, and the coefficients are identical to the page.

And how do you do it in R?

The rsm package: coded.data for the coding, rsm with FO, TWI, PQ or SO for the model, and summary gives the coefficients and the analysis of variance by groups of terms with lack of fit and pure error; canonical – the stationary point and the eigenvalues; contour – the maps; ccd and bbd – the designs. canonical zeroes eigenvalues smaller than a tenth of the largest before computing the point, so on a surface with a ridge its point differs from Montgomery's −½B⁻¹b used on the page. The ascent path of steepest is measured in uniform distance and not by a step in one factor. The script was run as is on the English sample data file, and the values in the comments are identical to the page. The site rule for R: a UTF-8 file, every comment line starts with an English word, and English variable names.

# Response surfaces – Professor Stat. This file is saved in UTF-8; every comment line starts with an English word
library(readxl)   # Reading Excel
library(rsm)      # Response surfaces (install.packages("rsm"))

# Loading: the standard template – row = run; the sheet "Data"
d <- as.data.frame(read_excel("C:/Rdata/DATA_rsm_example_eng.xlsx", sheet = "Data"))

# Coding: x = (ξ − center)/half-range – the factorial levels become −1 and +1
cd <- coded.data(d, x1 ~ (lesson_length - 45) / 10, x2 ~ (practice_items - 10) / 3, x3 ~ (review_weeks - 3) / 1)

# Second-order model: the rsm table – FO, TWI, PQ, residuals, lack of fit and pure error
fit <- rsm(retention ~ SO(x1, x2, x3), data = cd)
summary(fit)       # Coefficients: b0 71.74; x1 4.26, x2 −2.68, x3 1.69; x1:x2 1.10, x1:x3 0.23, x2:x3 −1.68; x1² −3.26, x2² −2.57, x3² −1.68; lack of fit: F(5, 5) = 0.89, p = .551

# Partial SS for each term (Table 2 on the page): explicit columns, then drop1 – each term after all the others
e <- transform(as.data.frame(cd), x12 = x1 * x2, x13 = x1 * x3, x23 = x2 * x3, x1q = x1^2, x2q = x2^2, x3q = x3^2)
drop1(lm(retention ~ x1 + x2 + x3 + x12 + x13 + x23 + x1q + x2q + x3q, data = e), test = "F")   # x1 F = 123.97 … x3q F = 21.07

# The stationary point and the eigenvalues (Table 4 on the page)
canonical(fit)     # xs = 0.568, −0.687, 0.883; in original units 50.7 minutes, 7.9 exercises, 3.9 weeks; eigenvalues −1.16, −2.64, −3.71 – a maximum

# Contour map for the first two factors, with the third at the stationary point
contour(fit, ~ x1 + x2, image = TRUE, at = canonical(fit)$xs)

# Model comparison (Table 1 on the page): adjusted R² and PRESS
for (f in list(retention ~ FO(x1, x2, x3), retention ~ FO(x1, x2, x3) + TWI(x1, x2, x3), retention ~ SO(x1, x2, x3))) {
  m <- rsm(f, data = cd)
  cat(sprintf("%-45s adjR2 = %.3f  PRESS = %.2f\n", deparse(f), summary(m)$adj.r.squared, sum((resid(m) / (1 - hatvalues(m)))^2)))
}   # Adjusted R² .482 / .432 / .944; PRESS 423.72 / 625.07 / 88.27

Reporting the results

The report should make it possible to reproduce the design and the analysis: the type of design (CCD, and its α; Box-Behnken), the factors and their ranges, the coding, the number of center points and blocks, and the number of runs; the model chosen and why (lack of fit, adjusted R² and predicted R²); a table of coefficients; and if there is one – the stationary point in original units, the predicted response there and its nature (the eigenvalues). APA (2020, section 6.43) asks for "sufficient information to allow the reader to corroborate the analyses". The format rules (section 6.36): R² and p without a zero before the decimal point; coefficients, means and eigenvalues – with one. The results paragraph in the output is built this way, in English; and before adopting the point – run confirmation runs at it.

Common mistakes

"One factor at a time". Each factor is moved separately and interactions and the true direction of improvement are missed – the path of steepest ascent changes all the factors together, in proportion to the coefficients (Montgomery, 2020, p. 415).

A second-order model from a design that does not support it. In a 2^k with center points all the quadratic terms are identical, and they cannot be estimated separately – only the overall curvature (p. 413). Add axial points.

Interpreting a stationary point outside the experimental region as if it were the optimum. Outside the region the model is not valid (p. 425).

Choosing the model by R² alone. R² always rises with more terms; adjusted R², predicted R² and the lack of fit say whether the addition is worthwhile (Table 11.7, p. 421).

Comparing coefficients in the natural equation. They depend on the units; compare in the coded one (p. 422).

Doing without center points. Without them there is no pure error, no lack-of-fit test, and no curvature check (pp. 412-413).

Rounding α. In Example 11.2 the table shows α = 1.414, and the output was computed with √2 – and the results differ in the third digit (a finding from validating the page). In the design use the exact α; if a level was rounded in practice, type in the value that was run.

References

American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000

Montgomery, D. C. (2020). Design and analysis of experiments (10th ed.). Wiley.

Box and Hunter (1957) and Box and Behnken (1960) – as cited in Montgomery (2020). According to the site rule, the reference list includes only sources that were read in the original; citations are by the printed page. Montgomery's chapter 11 was read in full from the page images (the text layer of the scan omits most of the equations). The computation on the page is verified against R 4.3.3 (lm, drop1, predict, hatvalues, rstudent, cooks.distance) and the rsm package (rsm, canonical, ccd, bbd), and against all the printed numbers in Examples 11.1 and 11.2 and in Tables 11.8 and 11.11.