Response Surface Modeling (RSM)

• Response Surface Modeling (RSM) is a statistical framework that approximates system responses by designing experiments and fitting regression models.
• It employs first- and second-order models, Gaussian processes, and classical designs like CCD and Box-Behnken to guide local optimization.
• RSM is applied across fields—from physical experiments to computer simulations and causal inference—offering robust strategies for prediction and uncertainty quantification.

Response Surface Modeling (RSM) denotes a family of statistical frameworks in which a response is represented as a function of controllable inputs and then used for prediction, interpretation, or optimization. In the classical literature, the response surface is a first- or second-order approximation to an unknown response over a region of interest, fitted by regression to data from designed experiments and then interrogated for directions of improvement or stationary points (Diaz-Garcia et al., 2012). In computer-experiment and machine-learning settings, the same idea extends to surrogate or metamodel construction with Gaussian processes and related emulators when direct evaluations are expensive (Toutiaee et al., 2021). The acronym is not uniform across disciplines: in high-contrast imaging, “RSM” denotes a “regime switching model” used for exoplanet detection rather than response surfaces (Dahlqvist et al., 2021).

1. Classical response-surface formulation

A standard RSM setup assumes a true but unknown response function over controllable factors $x=(x_1,\dots,x_n)'$, approximated locally by a quadratic polynomial. One canonical form is

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$

with matrix representation

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$

least-squares estimator

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

and fitted surface

$\hat y(x)=z'(x)\hat\beta$

(Diaz-Garcia et al., 2012). A closely related generic statement used in surrogate-based work is

$y=f'(x)\beta+\epsilon,$

where $f(x)$ is a vector of basis functions and $\epsilon$ is mean-zero random error (Afrin et al., 2018).

The local character of the approximation is central. When the region is roughly linear, a first-order model is used,

$$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \cdots + \beta_kx_k + e,$$

whereas curvature motivates a second-order specification,

$$y = \beta_0 + \sum_i \beta_i x_i + \sum_i \beta_{ii}x_i^2 + \sum_{i<j}\beta_{ij}x_ix_j + e$$

(Gamero-Salinas et al., 2024). In manufacturing optimization, the same principle appears in the two-factor second-order model

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$0

with linear, interaction, and quadratic terms explicitly separated (Srinivasan et al., 2020).

This formulation is broader than physical experimentation alone. In one line of work, the “true” response is the output of a trained black-box model rather than a physical process; the response surface is then a metamodel for the learned predictor itself (Toutiaee et al., 2021). A plausible implication is that RSM is best understood as a representation-and-optimization paradigm, not as a single fixed regression family.

2. Experimental design and sequential search

RSM is inseparable from design of experiments. The cited literature repeatedly treats optimization as a staged process: broad exploration, model fitting, movement toward improvement, and local refinement. One building-performance study follows an explicit sequence of fractional-factorial screening, first-order response-surface fitting, repeated movement along steepest ascent, and then a central composite design (CCD) once curvature appears (Gamero-Salinas et al., 2024). An actuarial neural-network study uses an analogous progression: a full factorial $$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$1 first-order screening design, a steepest-descent stage, and then a second-order CCD for local refinement of the hyperparameter region (Ariuntugs et al., 2024).

Classical named designs remain prominent. The manufacturing study notes that CCD and Box-Behnken designs are widely used, and implements a simplified two-factor CCD with factorial points and center points while omitting axial points for practical reasons (Srinivasan et al., 2020). In simulation-based tissue modeling, Latin Hypercube Sampling (LHS) is used as a stratified Monte Carlo design to generate well-distributed samples for fitting a quadratic surrogate, and is reported to give smaller estimation error than simple random sampling as sample size increases (Afrin et al., 2018).

The design problem also admits direct optimal-design formulations. For large linear and full quadratic response-surface models, the D-optimal design problem is written as

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$2

where each candidate run corresponds to a row $$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$3 of the design matrix (Ponte et al., 2023). In the linear case,

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$4

and in the full quadratic case,

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$5

(Ponte et al., 2023).

A major contemporary extension concerns restricted randomization. Response surface designs are often described under complete randomization, but many practical experiments are multi-stratum because some factors are hard to set or because units are crossed or nested. A recent general solution constructs designs stratum by stratum, with criterion matrices $$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$6 adapted to completely randomized, blocked, or row $$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$7 column strata, and uses compound criteria to balance parameter estimation, pure error, lack-of-fit, and treatment degrees of freedom efficiency (Trinca et al., 2024). This suggests that experimental structure is part of the response-surface problem rather than a secondary implementation detail.

3. Optimization, uncertainty, and robustness of the estimated optimum

RSM is often presented as producing a best operating condition, but the estimated optimum is itself a function of estimated regression coefficients and therefore has sampling variability. For quadratic RSM over a spherical region

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$8

the fitted-surface optimization problem is

$$y(x)=\beta_0+\sum_{i=1}^n \beta_i x_i+\sum_{i=1}^n \beta_{ii}x_i^2+\sum_{i<j}\beta_{ij}x_ix_j,$$9

and the optimizer is characterized by Kuhn–Tucker conditions for the Lagrangian

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$0

(Diaz-Garcia et al., 2012). In the unconstrained interior case, where $$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$1, the stationary point has the explicit form

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$2

(Diaz-Garcia et al., 2012).

The sensitivity of this optimizer to perturbations in $$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$3 yields an asymptotic distribution for the estimated optimum. If

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$4

then

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$5

(Diaz-Garcia et al., 2012). The practical implication stated in that work is that one can construct confidence regions or interval estimates for optimal operating conditions rather than relying only on a point estimate (Diaz-Garcia et al., 2012).

Robustness to model misspecification is a second major issue. One methodological framework assumes a primary polynomial model

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$6

where $$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$7 contains potential contamination terms, and then optimizes a compound criterion that simultaneously supports primary-model inference, lack-of-fit detection, and minimization of the variance and bias effects induced by omitted terms (Egorova et al., 2022). The same work adopts a model-independent pure-error approach to variance estimation and extends the resulting compound criteria to blocked and multistratum experiments via stratum-by-stratum construction and point-exchange search (Egorova et al., 2022). A common misconception is therefore that RSM design can be optimized purely for a single fitted polynomial; the cited literature treats design adequacy under contamination as a central objective rather than a peripheral diagnostic.

4. Surrogate, Gaussian-process, and multi-fidelity response surfaces

In computer experiments and black-box emulation, the response surface is frequently modeled by a Gaussian process rather than a low-order polynomial. One formulation writes

$$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$8

where $$y=X\beta+\varepsilon,\qquad \varepsilon\sim N_N(0,\sigma^2 I_N),$$9 is a trend function and $\hat\beta=(X'X)^{-1}X'y,$0 is a stationary Gaussian stochastic process with covariance

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

or, in Gaussian product form,

$\hat\beta=(X'X)^{-1}X'y,$2

(Toutiaee et al., 2021). In that framework, the correlation parameters $\hat\beta=(X'X)^{-1}X'y,$3 are estimated by maximum likelihood and interpreted as variable-importance indicators because they correspond to the correlation of individual variables with the target response (Toutiaee et al., 2021).

Sequential design can be made response-aware rather than merely space-filling. For expensive simulators modeled by a Gaussian process

$\hat\beta=(X'X)^{-1}X'y,$4

two contour-based strategies have been proposed for improving global prediction: Multiple Contour Estimation, which estimates several fixed contours simultaneously, and Sequential Contour Estimation, which adaptively chooses the next contour level as the predicted response at the point of maximum predictive variance (Yang et al., 2019). The empirical conclusion reported there is that contour-guided designs can outperform maximin Latin hypercube and several other sequential criteria, especially in maximum prediction error (Yang et al., 2019).

A related line of work uses prior surrogates from related problems. In aerodynamic computer experiments, database functions $\hat\beta=(X'X)^{-1}X'y,$5 are aligned by admissible affine transformations, compressed with POD or PCA into a low-dimensional basis, and then incorporated as a generic surrogate model inside hierarchical Kriging,

$\hat\beta=(X'X)^{-1}X'y,$6

so that sparse high-fidelity samples need only correct the residual mismatch (Rosenbaum et al., 2012). This suggests a reusable structural prior for families of related response surfaces.

Multi-fidelity RSM makes the same point in a different form. For two fidelities, the hierarchical surrogate

$\hat\beta=(X'X)^{-1}X'y,$7

is trained on nested designs $\hat\beta=(X'X)^{-1}X'y,$8, and the cited study proposes a heuristic budget-allocation rule derived from an empirical error grid over $\hat\beta=(X'X)^{-1}X'y,$9 combinations rather than from fidelity correlation alone (Rijn et al., 2021). In structural model updating, a multi-response Gaussian process is used to emulate the error response surface between finite-element predictions and measurements, and adaptive sampling iteratively concentrates new simulations in low-error regions of parameter space (Zhou et al., 2020). Across these formulations, “response surface” denotes an emulator of an expensive response, not necessarily a polynomial regression surface.

5. Extensions to complex inputs and structured factor spaces

Classical low-dimensional coordinates are not mandatory. For black-box optimization with high- or infinite-dimensional inputs $\hat y(x)=z'(x)\hat\beta$0, where $\hat y(x)=z'(x)\hat\beta$1 is a Hilbert space, RSM has been extended by combining dimension reduction with classical multivariate design of experiments. First- and second-order models take the forms

$\hat y(x)=z'(x)\hat\beta$2

and

$\hat y(x)=z'(x)\hat\beta$3

while design points are generated by lifting ordinary $\hat y(x)=z'(x)\hat\beta$4-dimensional designs into $\hat y(x)=z'(x)\hat\beta$5 through a chosen basis (Roche, 2015). A central theoretical observation in that work is that orthogonality, rotatability, and alphabetic optimality transfer directly because the model matrix in reduced coordinates is the same as in ordinary multivariate regression (Roche, 2015).

Order-of-addition experiments supply another nonstandard design space. There the response-surface formulation is built from component position numbers $\hat y(x)=z'(x)\hat\beta$6, standardized as

$\hat y(x)=z'(x)\hat\beta$7

The proposed second-order response-surface model for order effects is

$\hat y(x)=z'(x)\hat\beta$8

with $\hat y(x)=z'(x)\hat\beta$9 free parameters, placing it between pairwise-ordering and component-position models in complexity (Piepho et al., 2021). The cited examples show that this RS model can be competitive or decisively best, and that model averaging is advisable under model uncertainty (Piepho et al., 2021).

In nanotoxicology, response surfaces may be dose–duration maps rather than operating-condition polynomials. A Bayesian multivariate framework represents each nanoparticle–outcome surface as

$y=f'(x)\beta+\epsilon,$0

with linear B-splines, random interior knots, and hierarchical borrowing across outcomes (Patel et al., 2013). In that model, the first dose and time knots have direct risk-assessment meanings as maximal safe dose and maximal safe exposure time in the absence of interaction (Patel et al., 2013). This broadens the scope of RSM from optimization surfaces to probabilistic risk surfaces.

6. Applications, variant meanings, and disciplinary scope

The cited literature uses RSM across markedly different domains. In laser–tissue interaction modeling, a quadratic response surface trained on LHS-generated samples acts as the surrogate inside a surrogate-based optimization loop for maximum temperature and maximum thermal damage; the reported maximum relative error between the generalized DPL model and the surrogate model is $y=f'(x)\beta+\epsilon,$1 (Afrin et al., 2018). In diesel engine nozzle hydro-abrasive grinding, RSM is used as a local process-optimization tool with a two-factor CCD, and the study reports that abrasive liquid concentration had a positive effect on absolute error, with the response minimum at the low–low factor corner in the studied region (Srinivasan et al., 2020). In tropical housing optimization, RSM is coupled with Derringer–Suich desirability functions to minimize Indoor Overheating Hours and maximize Useful Daylight Illuminance, using screening, steepest ascent, CCD-based quadratic fitting, and $y=f'(x)\beta+\epsilon,$2 bootstrap replications for confidence intervals around the optimum (Gamero-Salinas et al., 2024). In actuarial machine learning, factorial design and RSM are used to optimize Combined Actuarial Neural Networks, reducing runs from $y=f'(x)\beta+\epsilon,$3 to $y=f'(x)\beta+\epsilon,$4 after dropping statistically insignificant hyperparameters, with negligible loss in out-of-sample Poisson deviance performance (Ariuntugs et al., 2024).

The term also appears in causal inference, where it does not refer to experimental design but to outcome regression. There the conditional response surface is

$y=f'(x)\beta+\epsilon,$5

and the average treatment effect is estimated by

$y=f'(x)\beta+\epsilon,$6

That study concludes that RSM performs well only when the outcome model is correctly specified, whereas augmented inverse probability weighting is more robust because of its doubly robust property (Das et al., 20 May 2026). This is a substantively different usage from design-based process optimization, but it preserves the core idea of fitting a surface for conditional response.

Finally, acronym ambiguity is real. In high-contrast exoplanet imaging, “RSM” refers not to response surfaces at all but to a regime switching model that uses residual cubes from PSF-subtraction methods to generate a planetary-regime probability map from ADI data (Dahlqvist et al., 2021). That contrast clarifies a broader point: across fields, “response surface” language may denote polynomial approximation, Gaussian-process emulation, hierarchical risk surfaces, or conditional outcome regression, while the acronym RSM can also be repurposed for unrelated models. The common statistical theme, when the term does mean Response Surface Modeling, is the construction of an explicit surface linking inputs to responses so that prediction, optimization, uncertainty quantification, or interpretability can be carried out in a structured way.