Response-surface regression

Description

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

Usage

rsm (formula, data, ...)

## S3 method for class 'rsm'
summary(object, adjust = rev(p.adjust.methods), ...)
## S3 method for class 'summary.rsm'
print(x, ...)

loftest (object)

canonical (object, threshold = 1e-04)
xs (object, ...)

## S3 method for class 'rsm'
codings(object)

Arguments

Details

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.

Value

rsm returns an rsm object, which is a lm object with additional members as follows:

summary is the summary method for rsm objects. It returns an object of class summary.rsm, which is an extension of the summary.lm class with these additional list elements:

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, and for second-order models, it shows the canonical analysis of the response surface.

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

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, thus modeling stationary ridges or valleys in their corresponding canonical directions. Setting a larger threshold may improve the numerical conditioning and bring the stationary point much closer to the design center, thus avoiding as much extrapolation. See vignette("rsm") for more details.

xs returns just the stationary point.

codings returns a list of coding formulas if the model was fitted to coded.data, or NULL otherwise.

Note

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.

Author(s)

Russell V. Lenth

References

Lenth RV (2009) “Response-Surface Methods in R, Using rsm”, Journal of Statistical Software, 32(7), 1–17. http://www.jstatsoft.org/v32/i07/.

See Also

FO, SO, lm, summary, coded.data

Examples

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)
canonical(RRex.rsm, threshold = 1)  # xs is MUCH closer to the experiment

## 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:
emmeans::emmeans(CR.rs2, ~ Time * Temp, mode = "decoded", 
        at = list(Time = c(80, 85, 90), Temp = c(165, 185)))

## End(Not run)