Fit a linear model with a response-surface component, and produce appropriate analyses and summaries.

`rsm (formula, data, ...)`# S3 method for rsm
summary(object, adjust = rev(p.adjust.methods), ...)
# S3 method for summary.rsm
print(x, ...)

# S3 method for rsm
codings(object)

loftest (object)

canonical (object, threshold = 0.1*max.eigen)
xs (object, ...)

formula

Formula to pass to `lm`

.
The model must include at least one `FO()`

, `SO()`

, `TWI()`

, or `PQ()`

term to
define the response-surface portion of the model.

data

`data`

argument to pass to `lm`

.

…

In `rsm`

, arguments that are passed to `lm`

,
`summary.lm`

, or `canonical`

, as appropriate.
In `summary`

, and `print`

, additional arguments
are passed to their generic methods.

object

An object of class `rsm`

adjust

Adjustment to apply to the P values in the coefficient matrix, chosen from among the available `p.adjust`

methods in the stats package. The default is `"none"`

.

threshold

Threshold for canonical analysis -- see "Canonical analysis" below.

x

An object produced by `summary`

`rsm`

returns an `rsm`

object, which is a `lm`

object with
additional members as follows:

The order of the model: 1 for first-order, 1.5 for first-order plus interactions, or 2 for a model that contains square terms.

The first-order response-surface coefficients.

The matrix of second-order response-surface coefficients, if present.

Labels for the response-surface terms. These make the summary much more readable.

Coding formulas, if provided in the `codings`

argument or
if the `data`

argument passed to `lm`

is a `coded.data`

object.

The `print`

method for `rsm`

objects just shows the call and the regression
coefficints.

The `summary`

method for `rsm`

objects returns an object of class
`summary.rsm`

, which is an extension of the `summary.lm`

class with these additional list elements:

- sa
Unit-length vector of the path of steepest ascent (first-order models only).

- canonical
Canonical analysis (second-order models only) from

`canonical`

- lof
ANOVA table including lack-of-fit test.

- coding
Coding formulas in parent

`rsm`

object.
Its print method shows the regression summary,
followed by an ANOVA and lack-of-fit test.
For first-order models, it shows the direction of
steepest ascent (see steepest), and for second-order models, it shows the canonical analysis of the
response surface.

`canonical`

returns a list with elements `xs`

, the stationary point, and `eigen`

, the eigenanalysis of the matrix **B** of second-order coefficients. Any eigenvalues less than `threshold`

are taken to be zero, and a message is displayed.
If this happens, the stationary point is determined using only the surviving eigenvectors,
and stationary ridges or valleys are assumed to exist in their
corresponding canonical directions. The default threshold is one tenth
of the maximum eigenvalue, internally named `max.eigen`

.
Setting a small `threshold`

may move the stationary point much farther from the origin.

When uncoded data are used, the canonical analysis and stationary point are not
very meaningful and those results should probably be ignored.
See `vignette("rsm")` for more details.

The function `xs`

returns just the stationary point.

`loftest`

returns an `anova`

object that tests the fitted model against a model
that interpolates the means of the response-surface-variable combinations.

`codings`

returns a `list`

of coding formulas if the model was fitted to
`coded.data`

, or `NULL`

otherwise.

Support is provided for the emmeans package: its `emmeans`

and related functions work with special provisions for models fitted to coded data. The optional `mode`

argument can have values of `"asis"`

(the default), `"coded"`

, or `"decoded"`

. The first two are equivalent and simply return LS means based on the original model formula and the variables therein (raw or coded), without any conversion. When coded data were used and the user specifies `mode = "decoded"`

, the user must specify results in terms of the decoded variables rather than the coded ones. See the illustration in the Examples section.

In `rsm`

, the model formula must contain at least an `FO`

term; optionally, you can add
one or more `TWI()`

terms and/or a `PQ()`

term. All variables that appear
in `TWI`

or `PQ`

*must* be included in `FO`

.
For convenience, specifying `SO()`

is the same as including `FO()`

, `TWI()`

, and `PQ()`

,
and is the safe, preferred way of specifying a full second-order model.

The variables in `FO`

comprise the variables to consider in response-surface methods. They need not all appear in `TWI`

and `PQ`

terms; and more than one `TWI`

term is allowed. For example, the following two model formulas are equivalent:

resp ~ Oper + FO(x1,x2,x3,x4) + TWI(x1,x2,x3) + TWI(x2,x3,x4) + PQ(x1,x3,x4) resp ~ Oper + FO(x1,x2,x3,x4) + TWI(formula = ~x1*x2*x3 + x2*x3*x4) + PQ(x1,x3,x4)

The first version, however, creates duplicate `x2:x3`

terms -- which `rsm`

can handle but there may be warning messages if it is subsequently used for predictions or plotted in `contour.lm`

.

In `summary.rsm`

, any `…`

arguments are passed to `summary.lm`

, except for `threshold`

, which is passed to `canonical`

.

Lenth RV (2009) ``Response-Surface Methods in R, Using rsm'',
*Journal of Statistical Software*, 32(7), 1--17.
https://www.jstatsoft.org/v32/i07/.

`FO`

, `SO`

,
`lm`

, `summary`

, `coded.data`

# NOT RUN { library(rsm) CR <- coded.data (ChemReact, x1~(Time-85)/5, x2~(Temp-175)/5) ### 1st-order model, using only the first block CR.rs1 <- rsm (Yield ~ FO(x1,x2), data=CR, subset=1:7) summary(CR.rs1) ### 2nd-order model, using both blocks CR.rs2 <- rsm (Yield ~ Block + SO(x1,x2), data=CR) summary(CR.rs2) ### Example of a rising-ridge situation from Montgomery et al, Table 6.2 RRex <- ccd(Response ~ A + B, n0 = c(0, 3), alpha = "face", randomize = FALSE, oneblock = TRUE) RRex$Response <- c(52.3, 5.3, 46.7, 44.2, 58.5, 33.5, 32.8, 49.2, 49.3, 50.2, 51.6) RRex.rsm <- rsm(Response ~ SO(A,B), data = RRex) canonical(RRex.rsm) # rising ridge is detected canonical(RRex.rsm, threshold = 0) # xs is far outside of the experimental region # } # NOT RUN { # Illustration of emmeans support emmeans::emmeans(CR.rs2, ~ x1 * x2, mode = "coded", at = list(x1 = c(-1, 0, 1), x2 = c(-2, 2))) # The following will yield the same results, but based on the decoded data emmeans::emmeans(CR.rs2, ~ Time * Temp, mode = "decoded", at = list(Time = c(80, 85, 90), Temp = c(165, 185))) # }