Response Surface Methodology in R: rsm Package, CCD & Box-Behnken

Response Surface Methodology (RSM) is a suite of techniques for finding the factor settings that maximize or minimize a response when you don't yet know the underlying equation. You run a planned experiment, fit a second-order polynomial to the results, then walk the surface to the peak.

This guide uses R's rsm package to build central composite and Box-Behnken designs, fit surface models, and extract the optimum in under twenty lines of code, with every block runnable right on the page.

How does response surface methodology work in R?

Factorial experiments answer which factors matter. RSM answers the next question, where exactly should we set them to hit the peak. The rsm package wraps design generation, second-order fitting, canonical analysis, and contour plots into one tight workflow. The block below fits a response surface to the classic ChemReact dataset and pulls out the optimum coordinates in four lines.

Read through the #> outputs. The canonical$xs row is what RSM is really about, the factor settings that maximize yield.

The stationary point (x1 = 0.35, x2 = 0.98) is where the fitted surface flattens out. In natural units that's roughly Time = 86.7 minutes and Temp = 179.9 degrees. One rsm() call, two factor settings, problem solved, at least for this model. The rest of the guide walks through how to get here from scratch on a new process, and how to decide whether the surface you fit is a peak, a bowl, a saddle, or a ridge.

Figure 1: The five-stage RSM workflow, from design to optimum.

When should you use CCD versus Box-Behnken?

Central Composite Designs (CCD) and Box-Behnken Designs (BBD) are the two workhorses of RSM. They both support second-order fits, but they place their runs very differently. Your process constraints, not statistical theory, pick the winner.

The quick decision table:

Let's get a concrete feel for the sizes. The block below builds one of each for three factors and reports the row count.

The CCD needs 20 runs (8 cube + 6 axial + 6 center), the BBD needs 15 (12 edge midpoints + 3 center). The BBD saves five experiments and avoids the four corner combinations entirely. In catalysis, that matters, a reactor run at maximum temperature and maximum pressure and maximum feed rate can destroy equipment.

Figure 2: CCD vs Box-Behnken: level count and where each design places its runs.

How do you build a central composite design with ccd()?

A CCD has three parts. The cube block is a 2^k factorial (or fractional factorial for large k) with factors at ±1. The star block adds 2k points along each axis at ±alpha, with alpha usually set so the design is rotatable. Center points sit at the origin and appear in both blocks for pure-error estimation.

ccd() takes four key arguments. basis is the number of factors (or a formula for fractional designs). n0 = c(cube_centers, star_centers) controls center points per block. alpha defaults to "rotatable", which picks the axial distance that makes prediction variance depend only on distance from the origin. coding is a list of formulas mapping coded units (x1, x2, ...) to your real-world variables.

The first eight rows are the 2^3 cube block, all ±1 combinations. Block 2 holds the axial (star) points at roughly ±1.68 on each axis plus centers. With alpha = "rotatable" the axial distance is (2^k)^(1/4) = 1.682 for three factors.

In practice, you'll want coded variables tied to your real units. Pass a list of formulas in coding. The object that comes back remembers the coding, so decode.data() or code2val() will flip between natural and coded units at any time.

decode.data() shows the settings the operator will actually dial in: temperature 140 or 160, pressure 45 or 55, feed rate 3 or 5. The star block's axial points will push past these boundaries when alpha > 1, which brings us to the warning below.

How do you build a Box-Behnken design with bbd()?

A Box-Behnken design takes a sharper shortcut. It skips the corners entirely and places runs only at the midpoints of each edge of the design cube, plus center points. Each factor sits at only three levels, -1, 0, +1. That's the full menu.

The tradeoff is clean. You lose the cube corners (and the combinations of extreme settings they encode), but you gain two things, you never run a combination of simultaneous extremes, and the design is orthogonal for the linear and interaction terms.

Every row varies two factors together, holding the third at its center. You never see the (Temp=160, Pres=55, Feed=5) combination that a CCD cube corner would require. For high-pressure, high-temperature chemistry, that's exactly the insurance you want.

How do you fit a second-order model with rsm() and SO()?

rsm() is essentially lm() with formula shortcuts. The shortcuts build the model-matrix columns you would otherwise have to write by hand:

• FO(x1, x2) adds first-order (linear) terms
• TWI(x1, x2) adds two-way interactions
• PQ(x1, x2) adds pure quadratics
• SO(x1, x2) = FO + TWI + PQ, the full second-order model

The second-order form is the workhorse of RSM. In math:

$$y = b_0 + \sum_{i=1}^{k} b_i x_i + \sum_{i=1}^{k} b_{ii} x_i^2 + \sum_{i

Where $b_0$ is the intercept, $b_i$ are linear slopes, $b_{ii}$ are curvatures, and $b_{ij}$ are interactions. The curvature terms are what let the model bend and therefore have an interior optimum rather than a monotonic slope.

The typical workflow is to start with a first-order fit, check whether the lack-of-fit F-test blows up (signalling curvature), then upgrade to second-order. Let's walk that path on the ChemReact data. We'll use the smaller first-stage dataset, ChemReact1, for the first-order fit.

The lack-of-fit p-value of 0.042 is a red flag. The linear model is missing structure, almost certainly curvature, so the response is already near the peak and a first-order approximation no longer captures the shape. Time to upgrade to a second-order fit using the full 10-run CCD from ChemReact.

The PQ (pure-quadratic) row jumps out with F = 13.97, p = 0.016, confirming real curvature. Lack-of-fit is no longer significant at the 5% level, meaning the second-order model captures the signal. The FO row is still highly significant (linear trend is still present), and TWI is quiet (no strong interaction between Time and Temperature on this surface).

How do you read the canonical analysis to find the optimum?

Once you have a second-order fit, summary() computes a canonical analysis. It finds the stationary point (the spot where the gradient of the fitted surface is zero) and computes the eigenvalues of the pure-quadratic matrix. The signs of those eigenvalues tell you the shape of the surface.

Pull the pieces straight out of the summary object.

Both eigenvalues are negative, so the stationary point is a genuine maximum. If one had been positive we'd be on a saddle, if one had been near zero we'd be on a ridge and would need follow-up experiments along the ridge direction to nail down the optimum.

Coded units are great for fitting, but the operator running the reactor needs natural units. code2val() converts from coded space back to Time/Temperature, and predict() evaluates the fitted yield at the optimum.

Run the process at Time = 86.7 min, Temp = 179.9°C and the model expects a yield near 84.6. That's the RSM deliverable, one number for each factor, validated by a second-order model fit.

How do you visualize the response surface?

Pictures close the loop. contour() draws 2D slices of the fitted surface, one for each pair of factors. persp() draws the 3D perspective version. Both accept the fitted model directly and use an at = argument to hold the other factors fixed while one pair varies.

The contour plot shades the yield landscape across the (x1, x2) plane and draws iso-yield contours. The brightest region (top-right quadrant in this fit) is where yield peaks, and the contour lines curve smoothly around the stationary point, the visual signature of a genuine maximum.

For a more intuitive shape, pair the contour with a 3D perspective view.

The rising mound near (x1 ≈ 0.35, x2 ≈ 1.0) is the same peak the canonical analysis pointed to. Seeing it as a 3D surface makes it easy to explain to non-statisticians who glazed over at the eigenvalues.

Practice Exercises

Exercise 1: Recover a known peak with a 3-factor CCD

Simulate a response surface with a known peak at (x1, x2, x3) = (0.5, -0.5, 0) and recover it using rsm(). Use the formula y = 50 - (x1 - 0.5)^2 - (x2 + 0.5)^2 - x3^2 + rnorm(n, 0, 0.5) to generate the response. Fit a second-order model and check that canonical$xs is within 0.1 of the true peak on each axis.

Exercise 2: Classify the surface from eigenvalues

You've fit a 4-factor BBD and found canonical$xs = c(x1 = 0.2, x2 = -0.3, x3 = 0.8, x4 = 0.1) with eigenvalues c(-1.2, -0.9, 0.3, -0.05). Without running more code, answer: (a) What kind of surface is this? (b) Which factor directions are worth exploring further? (c) Would you recommend running the process at xs?

Complete Example

End-to-end on a synthetic chemical reaction. Build a 2-factor CCD, simulate yields with a known peak at Time = 85, Temp = 175, fit a second-order model, extract the optimum in natural units.

The pipeline recovers the true peak (Time = 85, Temp = 175) within 0.05 on each axis, all from 11 simulated runs. In a real experiment the same five lines would turn n = 11 yield measurements into concrete operating settings to try on the plant.

Summary

• RSM = plan a design, run it, fit a second-order surface, locate the stationary point, verify.
• CCD uses 5 levels per factor, corners plus axial extensions, rotatable by default.
• Box-Behnken uses 3 levels, edge midpoints only, no extreme-corner runs.
• rsm(y ~ Block + SO(x1, x2, ...)) fits the full second-order model in one line.
• summary(model)$canonical$xs gives the stationary point in coded units, code2val() converts back.
• Eigenvalues classify the surface: all negative = max, all positive = min, mixed = saddle, near-zero = ridge.
• contour(model, ~ x1 + x2, image = TRUE) and persp(model, ~ x1 + x2) turn the fit into pictures.

References

1. Lenth, R.V. (2009). Response-Surface Methods in R, Using rsm. Journal of Statistical Software 32(7), 1-17. Link
2. Myers, R.H., Montgomery, D.C., & Anderson-Cook, C.M. (2016). Response Surface Methodology: Process and Product Optimization Using Designed Experiments (4th ed.). Wiley.
3. Box, G.E.P. & Draper, N.R. (2007). Response Surfaces, Mixtures, and Ridge Analyses (2nd ed.). Wiley.
4. Box, G.E.P. & Behnken, D.W. (1960). Some New Three Level Designs for the Study of Quantitative Variables. Technometrics 2(4), 455-475.
5. CRAN rsm package reference manual. Link
6. CRAN rsm vignette (JSS article). Link
7. Penn State STAT 503 Lesson 11: Response Surface Methods and Designs. Link

Continue Learning

• Factorial Experiments in R, the 2^k designs that feed the first-order screening stage of RSM.
• One-Way ANOVA in R, the F-test and lack-of-fit machinery rsm() uses under the hood.
• Linear Regression in R, the fitting engine that rsm() wraps for second-order models.