rsm: Response-surface regression

In rsm: Response-Surface Analysis

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, ...)

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

loftest (object)

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

Arguments

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

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:

order	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.
b	The first-order response-surface coefficients.
B	The matrix of second-order response-surface coefficients, if present.
labels	Labels for the response-surface terms. These make the summary much more readable.
coding	Coding formulas, if provided in the codings argument or if the data argument passed to lm is a coded.data object.

Summary and print methods

The print method for rsm objects just shows the call and the regression coefficints.

The summarymethod 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 analysis and stationary point

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.

Other functions

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.

emmeans support

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. doi: 10.18637/jss.v032.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)  # 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)))

## End(Not run)

Example output

Call:
rsm(formula = Yield ~ FO(x1, x2), data = CR, subset = 1:7)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 82.81429    0.54719 151.3456 1.143e-08 ***
x1           0.87500    0.72386   1.2088    0.2933    
x2           0.62500    0.72386   0.8634    0.4366    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.3555,	Adjusted R-squared:  0.0333 
F-statistic: 1.103 on 2 and 4 DF,  p-value: 0.4153

Analysis of Variance Table

Response: Yield
            Df Sum Sq Mean Sq F value  Pr(>F)
FO(x1, x2)   2 4.6250  2.3125  1.1033 0.41534
Residuals    4 8.3836  2.0959                
Lack of fit  2 8.2969  4.1485 95.7335 0.01034
Pure error   2 0.0867  0.0433                

Direction of steepest ascent (at radius 1):
       x1        x2 
0.8137335 0.5812382 

Corresponding increment in original units:
    Time     Temp 
4.068667 2.906191 


Call:
rsm(formula = Yield ~ Block + SO(x1, x2), data = CR)

             Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 84.095427   0.079631 1056.067 < 2.2e-16 ***
BlockB2     -4.457530   0.087226  -51.103 2.877e-10 ***
x1           0.932541   0.057699   16.162 8.444e-07 ***
x2           0.577712   0.057699   10.012 2.122e-05 ***
x1:x2        0.125000   0.081592    1.532    0.1694    
x1^2        -1.308555   0.060064  -21.786 1.083e-07 ***
x2^2        -0.933442   0.060064  -15.541 1.104e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9981,	Adjusted R-squared:  0.9964 
F-statistic: 607.2 on 6 and 7 DF,  p-value: 3.811e-09

Analysis of Variance Table

Response: Yield
            Df Sum Sq Mean Sq   F value    Pr(>F)
Block        1 69.531  69.531 2611.0950 2.879e-10
FO(x1, x2)   2  9.626   4.813  180.7341 9.450e-07
TWI(x1, x2)  1  0.063   0.063    2.3470    0.1694
PQ(x1, x2)   2 17.791   8.896  334.0539 1.135e-07
Residuals    7  0.186   0.027                    
Lack of fit  3  0.053   0.018    0.5307    0.6851
Pure error   4  0.133   0.033                    

Stationary point of response surface:
       x1        x2 
0.3722954 0.3343802 

Stationary point in original units:
     Time      Temp 
 86.86148 176.67190 

Eigenanalysis:
eigen() decomposition
$values
[1] -0.9233027 -1.3186949

$vectors
         [,1]       [,2]
x1 -0.1601375 -0.9870947
x2 -0.9870947  0.1601375


$xs
        A         B 
-5.176505 -2.706733 

$eigen
eigen() decomposition
$values
[1]  -0.509419 -12.706370

$vectors
        [,1]       [,2]
A -0.8396245 -0.5431673
B -0.5431673  0.8396245


$xs
        A         B 
-5.176505 -2.706733 

$eigen
eigen() decomposition
$values
[1]  -0.509419 -12.706370

$vectors
        [,1]       [,2]
A -0.8396245 -0.5431673
B -0.5431673  0.8396245


 x1 x2 emmean    SE df lower.CL upper.CL
 -1 -2   75.0 0.298  7     74.3     75.7
  0 -2   77.0 0.240  7     76.4     77.5
  1 -2   76.4 0.298  7     75.6     77.1
 -1  2   76.8 0.298  7     76.1     77.5
  0  2   79.3 0.240  7     78.7     79.9
  1  2   79.2 0.298  7     78.5     79.9

Results are averaged over the levels of: Block 
Confidence level used: 0.95 
 Time Temp emmean    SE df lower.CL upper.CL
   80  165   75.0 0.298  7     74.3     75.7
   85  165   77.0 0.240  7     76.4     77.5
   90  165   76.4 0.298  7     75.6     77.1
   80  185   76.8 0.298  7     76.1     77.5
   85  185   79.3 0.240  7     78.7     79.9
   90  185   79.2 0.298  7     78.5     79.9

Results are averaged over the levels of: Block 
Confidence level used: 0.95 

rsm documentation built on Oct. 7, 2021, 1:08 a.m.