Chapter 50 Response Surfaces Designs

50.1 Choosing Designs

The rsm packages has several functions that help you pick/generate central composite designs. Generally the assumption is that the design will be run in two blocks: (1) factorial points and some center points, and (2) axial points and some center points.

The ccd() function generates a CCD design. By default, it creates two blocks and chooses alpha for orthogonal blocking. Here we say not to randomize the design so that we can see the design more clearly.

> library(cfcdae)
Registered S3 method overwritten by 'DoE.base':
  method           from       
  factorize.factor conf.design
> library(rsm)
> ccd(3,randomize=FALSE)
   run.order std.order  x1.as.is  x2.as.is  x3.as.is Block
1          1         1 -1.000000 -1.000000 -1.000000     1
2          2         2  1.000000 -1.000000 -1.000000     1
3          3         3 -1.000000  1.000000 -1.000000     1
4          4         4  1.000000  1.000000 -1.000000     1
5          5         5 -1.000000 -1.000000  1.000000     1
6          6         6  1.000000 -1.000000  1.000000     1
7          7         7 -1.000000  1.000000  1.000000     1
8          8         8  1.000000  1.000000  1.000000     1
9          9         9  0.000000  0.000000  0.000000     1
10        10        10  0.000000  0.000000  0.000000     1
11        11        11  0.000000  0.000000  0.000000     1
12        12        12  0.000000  0.000000  0.000000     1
13         1         1 -1.825742  0.000000  0.000000     2
14         2         2  1.825742  0.000000  0.000000     2
15         3         3  0.000000 -1.825742  0.000000     2
16         4         4  0.000000  1.825742  0.000000     2
17         5         5  0.000000  0.000000 -1.825742     2
18         6         6  0.000000  0.000000  1.825742     2
19         7         7  0.000000  0.000000  0.000000     2
20         8         8  0.000000  0.000000  0.000000     2
21         9         9  0.000000  0.000000  0.000000     2
22        10        10  0.000000  0.000000  0.000000     2

Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is

You could tell it one block with a rotatable alpha.

> ccd(3,alpha="rotatable",oneblock=TRUE,randomize=FALSE)
   run.order std.order  x1.as.is  x2.as.is  x3.as.is
1          1         1 -1.000000 -1.000000 -1.000000
2          2         2  1.000000 -1.000000 -1.000000
3          3         3 -1.000000  1.000000 -1.000000
4          4         4  1.000000  1.000000 -1.000000
5          5         5 -1.000000 -1.000000  1.000000
6          6         6  1.000000 -1.000000  1.000000
7          7         7 -1.000000  1.000000  1.000000
8          8         8  1.000000  1.000000  1.000000
9          9         9  0.000000  0.000000  0.000000
10        10        10  0.000000  0.000000  0.000000
11        11        11  0.000000  0.000000  0.000000
12        12        12  0.000000  0.000000  0.000000
13         1         1 -1.681793  0.000000  0.000000
14         2         2  1.681793  0.000000  0.000000
15         3         3  0.000000 -1.681793  0.000000
16         4         4  0.000000  1.681793  0.000000
17         5         5  0.000000  0.000000 -1.681793
18         6         6  0.000000  0.000000  1.681793
19         7         7  0.000000  0.000000  0.000000
20         8         8  0.000000  0.000000  0.000000
21         9         9  0.000000  0.000000  0.000000
22        10        10  0.000000  0.000000  0.000000

Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is

You can specify a fractional factorial for the factorial part, in this case with the factorial part broken into two blocks.

> ccd(~A+B+C+D, generators=E~A*B*C*D, block = Day~A*B*C,randomize=FALSE)
    run.order std.order   A.as.is   B.as.is   C.as.is   D.as.is   E.as.is Day
1           1         1 -1.000000 -1.000000 -1.000000 -1.000000  1.000000   1
4           2         2  1.000000  1.000000 -1.000000 -1.000000  1.000000   1
6           3         3  1.000000 -1.000000  1.000000 -1.000000  1.000000   1
7           4         4 -1.000000  1.000000  1.000000 -1.000000  1.000000   1
9           5         5 -1.000000 -1.000000 -1.000000  1.000000 -1.000000   1
12          6         6  1.000000  1.000000 -1.000000  1.000000 -1.000000   1
14          7         7  1.000000 -1.000000  1.000000  1.000000 -1.000000   1
15          8         8 -1.000000  1.000000  1.000000  1.000000 -1.000000   1
11          9         9  0.000000  0.000000  0.000000  0.000000  0.000000   1
2          10        10  0.000000  0.000000  0.000000  0.000000  0.000000   1
3          11        11  0.000000  0.000000  0.000000  0.000000  0.000000   1
41         12        12  0.000000  0.000000  0.000000  0.000000  0.000000   1
13          1         1  1.000000 -1.000000 -1.000000 -1.000000 -1.000000   2
42          2         2 -1.000000  1.000000 -1.000000 -1.000000 -1.000000   2
61          3         3 -1.000000 -1.000000  1.000000 -1.000000 -1.000000   2
71          4         4  1.000000  1.000000  1.000000 -1.000000 -1.000000   2
91          5         5  1.000000 -1.000000 -1.000000  1.000000  1.000000   2
121         6         6 -1.000000  1.000000 -1.000000  1.000000  1.000000   2
141         7         7 -1.000000 -1.000000  1.000000  1.000000  1.000000   2
151         8         8  1.000000  1.000000  1.000000  1.000000  1.000000   2
111         9         9  0.000000  0.000000  0.000000  0.000000  0.000000   2
21         10        10  0.000000  0.000000  0.000000  0.000000  0.000000   2
31         11        11  0.000000  0.000000  0.000000  0.000000  0.000000   2
411        12        12  0.000000  0.000000  0.000000  0.000000  0.000000   2
16          1         1 -2.160247  0.000000  0.000000  0.000000  0.000000   3
22          2         2  2.160247  0.000000  0.000000  0.000000  0.000000   3
32          3         3  0.000000 -2.160247  0.000000  0.000000  0.000000   3
43          4         4  0.000000  2.160247  0.000000  0.000000  0.000000   3
5           5         5  0.000000  0.000000 -2.160247  0.000000  0.000000   3
62          6         6  0.000000  0.000000  2.160247  0.000000  0.000000   3
72          7         7  0.000000  0.000000  0.000000 -2.160247  0.000000   3
8           8         8  0.000000  0.000000  0.000000  2.160247  0.000000   3
92          9         9  0.000000  0.000000  0.000000  0.000000 -2.160247   3
10         10        10  0.000000  0.000000  0.000000  0.000000  2.160247   3
112        11        11  0.000000  0.000000  0.000000  0.000000  0.000000   3
122        12        12  0.000000  0.000000  0.000000  0.000000  0.000000   3
131        13        13  0.000000  0.000000  0.000000  0.000000  0.000000   3
142        14        14  0.000000  0.000000  0.000000  0.000000  0.000000   3

Data are stored in coded form using these coding formulas ...
A ~ A.as.is
B ~ B.as.is
C ~ C.as.is
D ~ D.as.is
E ~ E.as.is

Sometimes it’s difficult to pick the numbers of center points. We really want orthogonal blocking, but rotatable is also nice. Can we get an orthogonal design that is nearly rotatable?

ccd.pick() evaluates different numbers of center points for different blocks, computes the alpha for orthogonal, and tries to find combinations that are nearly orthogonal.

Here we just ask for three factors in cube and axial blocks. None of the designs is both rotatable and orthogonal, but the one that is closest has 9 center points in the cube block and 6 in the axial block. That could be too many, or the block size could be too large. The next one has only three total center points (2 and 1), which is probably too few.

> ccd.pick(3)
   n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s  N alpha.rot alpha.orth
1    8    9      1   6    6     1     1     1 29  1.681793   1.680336
2    8    2      1   6    1     1     1     1 17  1.681793   1.673320
3    8    6      1   6    4     1     1     1 24  1.681793   1.690309
4    8    5      1   6    3     1     1     1 22  1.681793   1.664101
5    8   10      1   6    7     1     1     1 31  1.681793   1.699673
6    8    8      1   6    5     1     1     1 27  1.681793   1.658312
7    8    3      1   6    2     1     1     1 19  1.681793   1.705606
8    8    7      1   6    5     1     1     1 26  1.681793   1.712698
9    8    4      1   6    2     1     1     1 20  1.681793   1.632993
10   8    4      1   6    3     1     1     1 21  1.681793   1.732051

You can get fancier and ask for the cube portion to be blocked further. Here was ask for two blocks of size 4 in the cube portion.

> ccd.pick(3,n.c=4,blks.c=2)
   n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s  N alpha.rot alpha.orth
1    4    1      2   6    1     1     1     1 17  1.681793   1.673320
2    4    6      2   6    8     1     1     1 34  1.681793   1.673320
3    4    3      2   6    4     1     1     1 24  1.681793   1.690309
4    4    5      2   6    7     1     1     1 31  1.681793   1.699673
5    4    4      2   6    5     1     1     1 27  1.681793   1.658312
6    4    7      2   6   10     1     1     1 38  1.681793   1.705606
7    4    7      2   6    9     1     1     1 37  1.681793   1.651446
8    4    2      2   6    2     1     1     1 20  1.681793   1.632993
9    4    2      2   6    3     1     1     1 21  1.681793   1.732051
10   4    4      2   6    6     1     1     1 28  1.681793   1.732051

rsm also has a bbd() function that will generate Box-Behnken designs.

> bbd(4,block=FALSE,randomize=FALSE)
   run.order std.order x1.as.is x2.as.is x3.as.is x4.as.is
1          1         1       -1       -1        0        0
2          2         2        1       -1        0        0
3          3         3       -1        1        0        0
4          4         4        1        1        0        0
5          5         5        0        0       -1       -1
6          6         6        0        0        1       -1
7          7         7        0        0       -1        1
8          8         8        0        0        1        1
9          9         9       -1        0        0       -1
10        10        10        1        0        0       -1
11        11        11       -1        0        0        1
12        12        12        1        0        0        1
13        13        13        0       -1       -1        0
14        14        14        0        1       -1        0
15        15        15        0       -1        1        0
16        16        16        0        1        1        0
17        17        17       -1        0       -1        0
18        18        18        1        0       -1        0
19        19        19       -1        0        1        0
20        20        20        1        0        1        0
21        21        21        0       -1        0       -1
22        22        22        0        1        0       -1
23        23        23        0       -1        0        1
24        24        24        0        1        0        1
25        25        25        0        0        0        0
26        26        26        0        0        0        0
27        27        27        0        0        0        0
28        28        28        0        0        0        0

Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is
x4 ~ x4.as.is

50.2 Potatoes data

From Eren and Kaymak-Erktekin (2007) via Lye (2019). We want optimum conditions to dehydrate potatoes via osmotic dehydration. This is a rotatable central composite design. There are several response variables, but we would like to maximize water loss (WL).

Variable Unit -2 -1 0 1 2
Temperature \(^o\)C 20 30 40 50 60
Sucrose conc. % 40 45 50 55 60
Salt conc. % 0 3.75 7.5 11.25 15
Time min 29.5 142 254.5 367 479.5

Read the data in. There are 16 factorial points, 8 axial points, and 7 center points for a total of 31 runs. They are entered in that order in the data set, although they would have been randomized when the experiment was run.

Then use coded.data() to create coded versions of the predictor variables within the data frame.

> library(emmeans)
> potatoes <- read.csv("Potatoes.csv",header=TRUE)
> cpotatoes <- coded.data(potatoes,cTemp~(Temp-40)/10,cSucrose~(Sucrose-50)/5,cSalt~(Salt-7.5)/3.75,cTime~(Time-254.4)/112.5)
> cpotatoes
   Temp Sucrose  Salt  Time   WL  SG   WR    aw
1    30      45  3.75 142.0 40.0 3.6 36.4 0.954
2    50      45  3.75 142.0 46.9 4.5 42.5 0.931
3    30      55  3.75 142.0 46.2 4.0 42.2 0.942
4    50      55  3.75 142.0 54.6 5.5 49.0 0.919
5    30      45 11.25 142.0 48.6 5.0 43.5 0.878
6    50      45 11.25 142.0 56.0 6.6 49.4 0.855
7    30      55 11.25 142.0 54.2 5.9 48.3 0.861
8    50      55 11.25 142.0 60.5 7.1 53.4 0.828
9    30      45  3.75 367.0 48.9 5.8 43.1 0.929
10   50      45  3.75 367.0 52.0 7.4 44.6 0.919
11   30      55  3.75 367.0 55.9 6.5 49.5 0.911
12   50      55  3.75 367.0 60.5 8.0 52.5 0.896
13   30      45 11.25 367.0 56.9 7.0 49.8 0.849
14   50      45 11.25 367.0 58.6 8.2 50.4 0.838
15   30      55 11.25 367.0 61.4 7.4 54.0 0.816
16   50      55 11.25 367.0 64.8 8.6 56.2 0.798
17   20      50  7.50 254.5 54.6 4.3 50.3 0.897
18   60      50  7.50 254.5 62.2 7.6 54.6 0.864
19   40      40  7.50 254.5 50.3 7.0 43.3 0.891
20   40      60  7.50 254.5 62.8 8.1 54.7 0.846
21   40      50  0.00 254.5 43.4 5.1 38.3 0.957
22   40      50 15.00 254.5 61.9 6.7 55.1 0.778
23   40      50  7.50  29.5 40.9 3.7 37.2 0.941
24   40      50  7.50 479.5 60.6 9.3 51.3 0.871
25   40      50  7.50 254.5 60.5 7.0 53.4 0.878
26   40      50  7.50 254.5 60.5 6.5 53.9 0.869
27   40      50  7.50 254.5 61.5 7.1 54.4 0.874
28   40      50  7.50 254.5 60.6 6.8 53.7 0.876
29   40      50  7.50 254.5 61.9 6.4 55.5 0.878
30   40      50  7.50 254.5 59.1 6.4 52.8 0.880
31   40      50  7.50 254.5 63.7 6.8 56.9 0.881

Data are stored in coded form using these coding formulas ...
cTemp ~ (Temp - 40)/10
cSucrose ~ (Sucrose - 50)/5
cSalt ~ (Salt - 7.5)/3.75
cTime ~ (Time - 254.4)/112.5

This is a central composite design, but let’s pretend for a moment that we only have the factorial points and some of the center points. That would be a good design for fitting a first order model.

Fit the first order model to this subset of the data. The FO(a,b,c) creates the first order model for the variables a, b, and c.

The lack of fit test is significant, meaning that there is enough curvature that the first order model is not sufficient. There are 12 df for lack of fit because there are 11 for interactions from the factorial portion of the data, plus 1 df for the contrast between the center points and the factorial points.

> use1 <- rep(TRUE,31)
> use1[17:24]<-FALSE # don't use axial points
> use1[29:31]<-FALSE # don't use last three center points
> WL.fo <- rsm(WL~FO(cTemp,cSucrose,cSalt,cTime),
+                        data=cpotatoes,subset=use1)
> summary(WL.fo)

Call:
rsm(formula = WL ~ FO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes, 
    subset = use1)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 55.45211    0.75314 73.6282 < 2.2e-16 ***
cTemp        2.61250    0.84203  3.1026 0.0072783 ** 
cSucrose     3.13750    0.84203  3.7261 0.0020278 ** 
cSalt        3.50000    0.84203  4.1566 0.0008438 ***
cTime        3.25000    0.84203  3.8597 0.0015430 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.7878,    Adjusted R-squared:  0.7312 
F-statistic: 13.92 on 4 and 15 DF,  p-value: 6.168e-05

Analysis of Variance Table

Response: WL
                                  Df Sum Sq Mean Sq F value    Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime)  4 631.71 157.926  13.921 6.168e-05
Residuals                         15 170.16  11.344                  
Lack of fit                       12 169.46  14.121  59.879  0.003106
Pure error                         3   0.71   0.236                  

Direction of steepest ascent (at radius 1):
    cTemp  cSucrose     cSalt     cTime 
0.4157757 0.4993287 0.5570201 0.5172330 

Corresponding increment in original units:
     Temp   Sucrose      Salt      Time 
 4.157757  2.496644  2.088825 58.188708 

Now we’ve added more runs to fill out the central composite design, so we can fit the second order model (OK, the data were already there, we’re just using all of them now).

The SO(a,b,c) form gives the full second order model: linear effects, pure quadratic effects, and two way interaction effects (the cross products). An equivalent, but less compact, form of the same model is FO(a,b,c)+TWI(a,b,c)+PQ(a,b,c) yielding first order, two way interaction, and pure quadratic effects.

> WL.so <- rsm(WL~SO(cTemp,cSucrose,cSalt,cTime),data=cpotatoes)

The summary() gives us tons of information. The first block is the usual summary for regression terms. The first order terms are all highly significant. One of the pure quadratics is marginally significant, but the others are highly significant. Only one of the two way interactions (cross product) terms is even marginally significant.

The next block of output groups the terms into FO, TWI, and PQ, testing them as groups. It also does the lack of fit test, which is not significant for these data.

The third block prints the stationary point for the model surface in both coded and original variables. This point is reasonably within the experimental region, so it’s not an extrapolation.

The final block of output is the canonical analysis. The eigenvalues here are all negative, indicating a maximum. The eigenvectors (canonical directions) are each dominated by a single variable. This means that the canonical directions are fairly similar to the variable axes, which happens when the TWI terms are relatively small.

> summary(WL.so)

Call:
rsm(formula = WL ~ SO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes)

               Estimate Std. Error  t value  Pr(>|t|)    
(Intercept)    61.11090    0.59429 102.8303 < 2.2e-16 ***
cTemp           2.37590    0.32095   7.4026 1.495e-06 ***
cSucrose        3.13321    0.32095   9.7622 3.835e-08 ***
cSalt           3.87540    0.32095  12.0747 1.881e-09 ***
cTime           3.81308    0.32095  11.8805 2.378e-09 ***
cTemp:cSucrose  0.22500    0.39309   0.5724 0.5750098    
cTemp:cSalt    -0.26250    0.39309  -0.6678 0.5137829    
cTemp:cTime    -1.01250    0.39309  -2.5758 0.0203150 *  
cSucrose:cSalt -0.53750    0.39309  -1.3674 0.1904068    
cSucrose:cTime  0.13750    0.39309   0.3498 0.7310529    
cSalt:cTime    -0.45000    0.39309  -1.1448 0.2691328    
cTemp^2        -0.75565    0.29403  -2.5700 0.0205556 *  
cSucrose^2     -1.21815    0.29403  -4.1429 0.0007646 ***
cSalt^2        -2.19315    0.29403  -7.4589 1.359e-06 ***
cTime^2        -2.66815    0.29403  -9.0743 1.042e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9731,    Adjusted R-squared:  0.9496 
F-statistic:  41.4 on 14 and 16 DF,  p-value: 7.048e-10

Analysis of Variance Table

Response: WL
                                   Df  Sum Sq Mean Sq  F value    Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime)   4 1079.46 269.865 109.1573 2.125e-11
TWI(cTemp, cSucrose, cSalt, cTime)  6   26.48   4.413   1.7851    0.1657
PQ(cTemp, cSucrose, cSalt, cTime)   4  326.92  81.730  33.0589 1.499e-07
Residuals                          16   39.56   2.472                   
Lack of fit                        10   27.03   2.703   1.2944    0.3911
Pure error                          6   12.53   2.088                   

Stationary point of response surface:
    cTemp  cSucrose     cSalt     cTime 
1.3693434 1.3056026 0.5966421 0.4380645 

Stationary point in original units:
      Temp    Sucrose       Salt       Time 
 53.693434  56.528013   9.737408 303.682260 

Eigenanalysis:
eigen() decomposition
$values
[1] -0.6064319 -1.1566542 -2.1874299 -2.8846030

$vectors
                [,1]       [,2]       [,3]        [,4]
cTemp     0.95448498  0.1724550  0.0542147 0.237230786
cSucrose  0.18569220 -0.9420455 -0.2794036 0.001550546
cSalt    -0.07927487  0.2524180 -0.9018507 0.341563955
cTime    -0.21952749 -0.1381839  0.3250533 0.909424661

We can do predictions using emmeans(). For example, here are the predictions for various salt and sucrose coded values when time and temperature are at the level of the stationary point.

> emmeans(WL.so,~cSalt+cSucrose,at=list(cSalt=c(-1,0,1),cSucrose=c(-1,0,1),cTime=1.37,cTemp=.44))
NOTE: Results may be misleading due to involvement in interactions
 cSalt cSucrose emmean    SE df lower.CL upper.CL
    -1       -1   51.1 1.218 16     48.5     53.7
     0       -1   57.0 0.970 16     54.9     59.0
     1       -1   58.5 1.218 16     55.9     61.0
    -1        0   56.3 0.970 16     54.2     58.3
     0        0   61.6 0.755 16     60.0     63.2
     1        0   62.6 0.970 16     60.5     64.6
    -1        1   59.0 1.218 16     56.4     61.6
     0        1   63.8 0.970 16     61.8     65.9
     1        1   64.2 1.218 16     61.6     66.8

Confidence level used: 0.95 

Most of the two way interactions are not significant. Does it make sense to fit a reduced model including only the significant terms? Yes, because including non-significant terms in a model makes predictions worse. If you want more precise prediction, you don’t want the noise terms in your model.

See below how we can select just the terms we want. If you look at the help for rsm() you’ll see that the TWI() form allows a lot of flexibility (that we don’t need here with just a single cross product term).

> WL.red <- rsm(WL~FO(cTemp,cSucrose,cSalt,cTime)+PQ(cTemp,cSucrose,cSalt,cTime)+TWI(cTime,cTemp),data=cpotatoes)
> summary(WL.red)

Call:
rsm(formula = WL ~ FO(cTemp, cSucrose, cSalt, cTime) + PQ(cTemp, 
    cSucrose, cSalt, cTime) + TWI(cTime, cTemp), data = cpotatoes)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 61.11090    0.58107 105.1695 < 2.2e-16 ***
cTemp        2.37590    0.31381   7.5710 1.970e-07 ***
cSucrose     3.13333    0.31381   9.9847 1.990e-09 ***
cSalt        3.87500    0.31381  12.3481 4.290e-11 ***
cTime        3.81308    0.31381  12.1507 5.782e-11 ***
cTemp^2     -0.75565    0.28749  -2.6284 0.0157053 *  
cSucrose^2  -1.21815    0.28749  -4.2372 0.0003685 ***
cSalt^2     -2.19315    0.28749  -7.6285 1.750e-07 ***
cTime^2     -2.66815    0.28749  -9.2808 7.042e-09 ***
cTime:cTemp -1.01250    0.38434  -2.6344 0.0155014 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9663,    Adjusted R-squared:  0.9518 
F-statistic: 66.89 on 9 and 21 DF,  p-value: 1.964e-13

Analysis of Variance Table

Response: WL
                                  Df  Sum Sq Mean Sq  F value    Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime)  4 1079.46 269.865 114.1799 6.236e-14
PQ(cTemp, cSucrose, cSalt, cTime)  4  326.92  81.730  34.5800 5.814e-09
TWI(cTime, cTemp)                  1   16.40  16.402   6.9399    0.0155
Residuals                         21   49.63   2.364                   
Lack of fit                       15   37.10   2.474   1.1847    0.4444
Pure error                         6   12.53   2.088                   

Stationary point of response surface:
    cTemp  cSucrose     cSalt     cTime 
1.2525887 1.2860982 0.8834306 0.4768896 

Stationary point in original units:
     Temp   Sucrose      Salt      Time 
 52.52589  56.43049  10.81286 308.05008 

Eigenanalysis:
eigen() decomposition
$values
[1] -0.6299144 -1.2181548 -2.1931548 -2.7938951

$vectors
               [,1] [,2] [,3]      [,4]
cTemp     0.9705122    0    0 0.2410519
cSucrose  0.0000000   -1    0 0.0000000
cSalt     0.0000000    0   -1 0.0000000
cTime    -0.2410519    0    0 0.9705122

The results are very similar (the same coefficients, for example) in the reduced model. But the proof of the pudding is in the predictions. If you compare the predictions from the reduced model, they are somewhat different and have smaller standard errors.

> emmeans(WL.red,~cSalt+cSucrose,at=list(cSalt=c(-1,0,1),cSucrose=c(-1,0,1),cTime=1.37,cTemp=.44))
NOTE: Results may be misleading due to involvement in interactions
 cSalt cSucrose emmean    SE df lower.CL upper.CL
    -1       -1   51.2 0.813 21     49.5     52.9
     0       -1   57.3 0.771 21     55.7     58.9
     1       -1   58.9 0.813 21     57.3     60.6
    -1        0   55.5 0.771 21     53.9     57.2
     0        0   61.6 0.739 21     60.1     63.2
     1        0   63.3 0.771 21     61.7     64.9
    -1        1   57.5 0.813 21     55.8     59.2
     0        1   63.5 0.771 21     61.9     65.1
     1        1   65.2 0.813 21     63.5     66.9

Confidence level used: 0.95 

The lack of fit test says everything is OK, but what haven’t we done? We haven’t checked residuals!

So we do, and it looks bad. We are over estimating the high and low values and under estimating the middle values. (You can see a much weaker form of the same pattern in the full model residuals.)

> plot(WL.red,which=1)
Residual plot for reducted water loss model

Figure 50.1: Residual plot for reducted water loss model

Box-Cox picks the power 2, but to my great surprise the power 1 is just inside the interval. Things are a little better on the squared scale, but we’re going to continue with the original data to control the length of the example.

50.2.1 But wait! There’s more!

I hadn’t given you the full problem. The full problem is to maximize water loss while minimizing solid gain (SL). Let’s interpret “minimizing sold gain” to mean “keep solid gain under 6.5”.

We need to model solid gain. First fit the full model. Then fit a reduced model including only terms that are reasonably significant.

> SG.so <-  rsm(SG~SO(cTemp,cSucrose,cSalt,cTime),data=cpotatoes)
> summary(SG.so)

Call:
rsm(formula = SG ~ SO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes)

                Estimate Std. Error t value  Pr(>|t|)    
(Intercept)     6.713252   0.133778 50.1822 < 2.2e-16 ***
cTemp           0.720817   0.072248  9.9769 2.836e-08 ***
cSucrose        0.295872   0.072248  4.0952  0.000845 ***
cSalt           0.571028   0.072248  7.9037 6.487e-07 ***
cTime           1.162642   0.072248 16.0923 2.651e-11 ***
cTemp:cSucrose  0.006250   0.088486  0.0706  0.944565    
cTemp:cSalt    -0.018750   0.088486 -0.2119  0.834862    
cTemp:cTime     0.018750   0.088486  0.2119  0.834862    
cSucrose:cSalt -0.031250   0.088486 -0.3532  0.728575    
cSucrose:cTime -0.043750   0.088486 -0.4944  0.627722    
cSalt:cTime    -0.218750   0.088486 -2.4722  0.025036 *  
cTemp^2        -0.217113   0.066188 -3.2802  0.004711 ** 
cSucrose^2      0.182887   0.066188  2.7631  0.013855 *  
cSalt^2        -0.229613   0.066188 -3.4691  0.003164 ** 
cTime^2        -0.079613   0.066188 -1.2028  0.246542    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9676,    Adjusted R-squared:  0.9392 
F-statistic:  34.1 on 14 and 16 DF,  p-value: 3.079e-09

Analysis of Variance Table

Response: SG
                                   Df Sum Sq Mean Sq  F value    Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime)   4 54.825 13.7062 109.4090 2.088e-11
TWI(cTemp, cSucrose, cSalt, cTime)  6  0.824  0.1373   1.0959 0.4064144
PQ(cTemp, cSucrose, cSalt, cTime)   4  4.164  1.0411   8.3102 0.0007954
Residuals                          16  2.004  0.1253                   
Lack of fit                        10  1.516  0.1516   1.8615 0.2307223
Pure error                          6  0.489  0.0814                   

Stationary point of response surface:
     cTemp   cSucrose      cSalt      cTime 
 2.7205417  0.5983824 -7.1207671 17.2405071 

Stationary point in original units:
      Temp    Sucrose       Salt       Time 
  67.20542   52.99191  -19.20288 2193.95705 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.18483027 -0.02185009 -0.21777830 -0.28865425

$vectors
                [,1]        [,2]        [,3]        [,4]
cTemp    -0.00642643 -0.06564489  0.99643763 -0.05254991
cSucrose -0.99701518 -0.05701335 -0.01284677 -0.05044965
cSalt     0.01799251  0.46973620 -0.01547838 -0.88248773
cTime     0.07480440 -0.87851473 -0.08189887 -0.46465982
> SG.red <- rsm(SG~FO(cTemp,cSucrose,cSalt,cTime)+PQ(cTemp,cSucrose,cSalt)+TWI(cSalt,cTime),data=cpotatoes)
> summary(SG.red)

Call:
rsm(formula = SG ~ FO(cTemp, cSucrose, cSalt, cTime) + PQ(cTemp, 
    cSucrose, cSalt) + TWI(cSalt, cTime), data = cpotatoes)

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  6.631945   0.104163 63.6688 < 2.2e-16 ***
cTemp        0.720833   0.065189 11.0576 1.880e-10 ***
cSucrose     0.295833   0.065189  4.5381 0.0001621 ***
cSalt        0.571028   0.065189  8.7596 1.268e-08 ***
cTime        1.162500   0.065189 17.8328 1.443e-14 ***
cTemp^2     -0.208644   0.059382 -3.5136 0.0019592 ** 
cSucrose^2   0.191356   0.059382  3.2225 0.0039190 ** 
cSalt^2     -0.221144   0.059382 -3.7241 0.0011793 ** 
cSalt:cTime -0.218750   0.079840 -2.7399 0.0119567 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9637,    Adjusted R-squared:  0.9505 
F-statistic: 73.01 on 8 and 22 DF,  p-value: 4.738e-14

Analysis of Variance Table

Response: SG
                                  Df Sum Sq Mean Sq  F value    Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime)  4 54.825 13.7062 134.3884 4.015e-15
PQ(cTemp, cSucrose, cSalt)         3  3.983  1.3277  13.0177 4.210e-05
TWI(cSalt, cTime)                  1  0.766  0.7656   7.5069   0.01196
Residuals                         22  2.244  0.1020                   
Lack of fit                       16  1.755  0.1097   1.3472   0.37622
Pure error                         6  0.489  0.0814                   

Stationary point of response surface:
     cTemp   cSucrose      cSalt      cTime 
 1.7274272 -0.7729905  5.3142857 -8.1344592 

Stationary point in original units:
      Temp    Sucrose       Salt       Time 
  57.27427   46.13505   27.42857 -660.72666 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.19135638  0.04495637 -0.20864362 -0.26609999

$vectors
         [,1]       [,2] [,3]      [,4]
cTemp       0  0.0000000    1 0.0000000
cSucrose    1  0.0000000    0 0.0000000
cSalt       0 -0.3801685    0 0.9249173
cTime       0  0.9249173    0 0.3801685

What we’re going to do is make some contour plots of our WL predictions and overlay contour plots of our SG predictions. Then we look at how well we can do. If we’re lucky, high values of WL are in the same regions as low values of SG, but we’re not lucky.

In this first plot we look at the two variables as time and temperature vary, with sucrose and salt set to the stationary point values for WL. The black countours are WL; the red contours are SG. The SL=6.5 contour curves from the southeast up to the northwest. At one point it gets between the 60 and 65 contours of WL. The best we can do in this slice is about 61 for WL.

> contour(WL.red,cTime~cTemp,at=c(cSalt=.88,cSucrose=1.29))
> contour(SG.red,cTime~cTemp,at=c(cSalt=.88,cSucrose=1.29),add=TRUE,col="red")
Contour plot of SG overlaid on WL; by Time and Temperature

Figure 50.2: Contour plot of SG overlaid on WL; by Time and Temperature

Now flip around and look at the results by salt and sucrose with time and temperature set to the stationary point levels.

We can clearly see the saddle point nature of SG in this slice. The 6.5 contour for SG is way over in the west southwest. It does get up between 50 and 55 on WL, but no higher. Clearly, the other point we found is better.

> contour(WL.red,cSucrose~cSalt,at=c(cTime=.48,cTemp=1.25))
> contour(SG.red,cSucrose~cSalt,at=c(cTime=.48,cTemp=1.25),add=TRUE,col="red")
Contour plot of SG overlaid on WL; by Sucrose and Salt

Figure 50.3: Contour plot of SG overlaid on WL; by Sucrose and Salt

50.3 Conversion data

These are data taken from the examples in the Design-Expert manual. We are making a product that should have an activity of 63; we could accept the range 60 to 66, but 63 is more desirable. There are three process variables: time (minutes), temperature (degrees C), and catalyst (percent). When we do our experiment, we have two responses: the conversion factor (essentially the percent yield) and the activity.

Our goal is to get process variables so that we hit an activity of 63 (or close enough) and maximize the conversion. In addition, we want to operate within the “cube” part of the experimental region.

Read in the data and covert to coded units.

> conversion <- read.csv("Conversion.csv",header=TRUE)
> conversion
    Time Temperature Catalyst Conversion Activity
1  45.00       85.00     2.50         75     60.4
2  45.00       85.00     2.50         83     60.6
3  45.00       85.00     2.50         76     59.1
4  45.00       85.00     2.50         81     59.2
5  45.00       85.00     2.50         91     58.9
6  45.00       85.00     2.50         80     60.8
7  53.41       85.00     2.50         79     65.9
8  36.59       85.00     2.50         76     53.6
9  45.00       85.00     1.66         55     57.4
10 45.00       85.00     3.34         81     63.2
11 45.00       93.41     2.50         97     60.7
12 45.00       76.59     2.50         85     60.0
13 40.00       80.00     2.00         74     53.2
14 40.00       90.00     3.00         66     59.8
15 50.00       90.00     2.00         70     62.6
16 40.00       90.00     2.00         88     53.4
17 50.00       90.00     3.00         97     67.8
18 50.00       80.00     3.00         90     67.9
19 40.00       80.00     3.00         71     57.3
20 50.00       80.00     2.00         51     62.9
> c.conversion <- coded.data(conversion,cTime~(Time-45)/5,cTemp~(Temperature-85)/5,
+   cCat~(Catalyst-2.5)/.5)

Let’s fit conversion starting with a second order model. Looking at the summary, one pure quadratic is not even marginally significant; one cross product is not marginally significant; and one first order term is not marginally significant. However, all variables need to be retained due to hierarchy.

> conversion.so <- rsm(Conversion~SO(cTime,cTemp,cCat),data=c.conversion)
> summary(conversion.so)

Call:
rsm(formula = Conversion ~ SO(cTime, cTemp, cCat), data = c.conversion)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  81.0917     1.9233 42.1632 1.351e-12 ***
cTime         1.0284     1.2760  0.8059 0.4390217    
cTemp         4.0403     1.2760  3.1664 0.0100492 *  
cCat          6.2060     1.2766  4.8612 0.0006602 ***
cTime:cTemp   2.1250     1.6673  1.2745 0.2312930    
cTime:cCat   11.3750     1.6673  6.8225 4.611e-05 ***
cTemp:cCat   -3.8750     1.6673 -2.3242 0.0424712 *  
cTime^2      -1.8311     1.2419 -1.4744 0.1711487    
cTemp^2       2.9407     1.2419  2.3678 0.0394226 *  
cCat^2       -5.2027     1.2443 -4.1811 0.0018843 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:   0.92, Adjusted R-squared:  0.8479 
F-statistic: 12.77 on 9 and 10 DF,  p-value: 0.0002207

Analysis of Variance Table

Response: Conversion
                        Df  Sum Sq Mean Sq F value   Pr(>F)
FO(cTime, cTemp, cCat)   3  762.93  254.31 11.4357 0.001428
TWI(cTime, cTemp, cCat)  3 1191.37  397.12 17.8576 0.000243
PQ(cTime, cTemp, cCat)   3  601.51  200.50  9.0160 0.003414
Residuals               10  222.38   22.24                 
Lack of fit              5   56.38   11.28  0.3397 0.869479
Pure error               5  166.00   33.20                 

Stationary point of response surface:
     cTime      cTemp       cCat 
-1.0179845 -0.5294920 -0.3192329 

Stationary point in original units:
       Time Temperature    Catalyst 
  39.910077   82.352540    2.340384 

Eigenanalysis:
eigen() decomposition
$values
[1]  3.408437  2.322758 -9.824322

$vectors
            [,1]      [,2]       [,3]
cTime  0.1380685 0.7982920  0.5862312
cTemp -0.9428308 0.2872127 -0.1690534
cCat   0.3033270 0.5293758 -0.7923093

Refit omitting the noise terms. Notice how you can add just certain two way interactions in the TWI() form. Residuals seem OK.

> conversion.red <- rsm(Conversion~FO(cTime,cTemp,cCat)+PQ(cTemp,cCat)+TWI(formula=~cTime*cCat+cTemp*cCat),data=c.conversion)
> summary(conversion.red)

Call:
rsm(formula = Conversion ~ FO(cTime, cTemp, cCat) + PQ(cTemp, 
    cCat) + TWI(formula = ~cTime * cCat + cTemp * cCat), data = c.conversion)

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  79.5933     1.7509 45.4589 8.375e-15 ***
cTime         1.0284     1.3683  0.7516 0.4667826    
cTemp         4.0403     1.3683  2.9529 0.0120800 *  
cCat          6.2060     1.3690  4.5334 0.0006857 ***
cTemp^2       3.1226     1.3252  2.3564 0.0362905 *  
cCat^2       -5.0213     1.3278 -3.7817 0.0026153 ** 
cTime:cCat   11.3750     1.7878  6.3624 3.597e-05 ***
cCat:cTemp   -3.8750     1.7878 -2.1674 0.0510298 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.8896,    Adjusted R-squared:  0.8251 
F-statistic: 13.81 on 7 and 12 DF,  p-value: 7.061e-05

Analysis of Variance Table

Response: Conversion
                                            Df  Sum Sq Mean Sq F value
FO(cTime, cTemp, cCat)                       3  762.93  254.31  9.9453
PQ(cTemp, cCat)                              2  553.16  276.58 10.8163
TWI(formula = ~cTime * cCat + cTemp * cCat)  2 1155.25  577.62 22.5891
Residuals                                   12  306.85   25.57        
Lack of fit                                  7  140.85   20.12  0.6061
Pure error                                   5  166.00   33.20        
                                               Pr(>F)
FO(cTime, cTemp, cCat)                       0.001418
PQ(cTemp, cCat)                              0.002063
TWI(formula = ~cTime * cCat + cTemp * cCat) 8.545e-05
Residuals                                            
Lack of fit                                  0.736262
Pure error                                           

Stationary point of response surface:
      cTime       cTemp        cCat 
-0.86490534 -0.70305722 -0.09040788 

Stationary point in original units:
       Time Temperature    Catalyst 
  40.675473   81.484714    2.454796 

Eigenanalysis:
eigen() decomposition
$values
[1]  4.588664  2.460161 -8.947536

$vectors
            [,1]       [,2]       [,3]
cTime  0.5989293 -0.5989112  0.5315910
cTemp -0.6385820 -0.7577546 -0.1342427
cCat   0.4832150 -0.2590625 -0.8362953
> plot(conversion.red,which=1)
Residual plot of reduced conversion model

Figure 50.4: Residual plot of reduced conversion model

Now fit activity. The pure quadratics and interactions are not significant.

> activity.so <- rsm(Activity~SO(cTime,cTemp,cCat),data=c.conversion)
> summary(activity.so)

Call:
rsm(formula = Activity ~ SO(cTime, cTemp, cCat), data = c.conversion)

            Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 59.81833    0.42544 140.6029 < 2.2e-16 ***
cTime        4.26033    0.28226  15.0935 3.293e-08 ***
cTemp        0.25460    0.28226   0.9020    0.3883    
cCat         2.23118    0.28240   7.9008 1.313e-05 ***
cTime:cTemp -0.38750    0.36881  -1.0507    0.3181    
cTime:cCat  -0.03750    0.36881  -0.1017    0.9210    
cTemp:cCat   0.31250    0.36881   0.8473    0.4166    
cTime^2      0.06769    0.27472   0.2464    0.8104    
cTemp^2      0.27977    0.27472   1.0184    0.3325    
cCat^2       0.26294    0.27525   0.9553    0.3620    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multiple R-squared:  0.9672,    Adjusted R-squared:  0.9376 
F-statistic: 32.74 on 9 and 10 DF,  p-value: 2.955e-06

Analysis of Variance Table

Response: Activity
                        Df Sum Sq Mean Sq F value    Pr(>F)
FO(cTime, cTemp, cCat)   3 316.71 105.571 97.0169 1.079e-07
TWI(cTime, cTemp, cCat)  3   1.99   0.665  0.6107    0.6232
PQ(cTime, cTemp, cCat)   3   1.94   0.645  0.5931    0.6336
Residuals               10  10.88   1.088                  
Lack of fit              5   7.23   1.446  1.9786    0.2359
Pure error               5   3.65   0.731                  

Stationary point of response surface:
    cTime     cTemp      cCat 
 15.05168  17.57008 -13.61031 

Stationary point in original units:
       Time Temperature    Catalyst 
 120.258424  172.850398   -4.305155 

Eigenanalysis:
eigen() decomposition
$values
[1]  0.48931661  0.18268588 -0.06160297

$vectors
            [,1]       [,2]       [,3]
cTime  0.3691679  0.4710760  0.8011258
cTemp -0.7502855 -0.3576299  0.5560328
cCat  -0.5484402  0.8063425 -0.2214160

Refit with the first order and also omitting unneeded first order term. The reduced model is acceptable from an ANOVA perspective, and residuals look OK.

> activity.fo <- rsm(Activity~FO(cTime,cTemp,cCat),data=c.conversion)
> activity.red <- rsm(Activity~FO(cTime,cCat),data=c.conversion)
> anova(activity.red,activity.fo,activity.so)
Analysis of Variance Table

Model 1: Activity ~ FO(cTime, cCat)
Model 2: Activity ~ FO(cTime, cTemp, cCat)
Model 3: Activity ~ FO(cTime, cTemp, cCat) + TWI(cTime, cTemp, cCat) + 
    PQ(cTime, cTemp, cCat)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     17 15.697                           
2     16 14.812  1    0.8853 0.8136 0.3883
3     10 10.882  6    3.9299 0.6019 0.7240
> plot(activity.red,which=1)
Residual plot for first order activity model

Figure 50.5: Residual plot for first order activity model

We have our models, now we work on the experimental goal. We are going to approach this graphically. For each level of cTemp from -1 to 1 by increments of .25 we will:

  1. Make a contour plot of activity with contours at levels 60, 63, and 66 (the minimum, desired, and maximum levels);
  2. Overlay a contour plot of conversion (red, to make it visible);
  3. Overlay a box (in green) to show the experimental area.

The goal is to walk through the images and find where we want to operate the system. This will involve some tradeoff between being as close as possible to 63 for activity and having the highest conversion.

Note: you can do all these within a loop, but I’ve kept them separate to get better alt text for the graphics.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-1),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at -1.

Figure 50.6: Overlaid contours of activity and conversion with coded temperature at -1.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.75),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at -.75.

Figure 50.7: Overlaid contours of activity and conversion with coded temperature at -.75.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.5),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at -.5.

Figure 50.8: Overlaid contours of activity and conversion with coded temperature at -.5.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.25),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at -.25.

Figure 50.9: Overlaid contours of activity and conversion with coded temperature at -.25.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=0),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at 0.

Figure 50.10: Overlaid contours of activity and conversion with coded temperature at 0.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.25),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at .25.

Figure 50.11: Overlaid contours of activity and conversion with coded temperature at .25.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.5),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at .5.

Figure 50.12: Overlaid contours of activity and conversion with coded temperature at .5.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.75),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at .75.

Figure 50.13: Overlaid contours of activity and conversion with coded temperature at .75.

> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=1),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")
Overlaid contours of activity and conversion with coded temperature at 1.

Figure 50.14: Overlaid contours of activity and conversion with coded temperature at 1.

In general, we can achieve higher conversions and stay within our activity bounds as we move to higher levels of temperature. If we stick to activity right at 63, we can almost hit conversion of 90. If we allow ourselves an activity of 66, we can hit a conversion of about 95.

50.4 Constrained maximization

Optimizing a quadratic function subject to linear equality and/or inequality constraints is called quadratic programming. If you want to maximize the function, use -2 times matrix of second order terms and the linear terms as the first two arguments. If you want to minimize, use 2B and -b instead.

The constraints are “greater than or equal to” constraints. The third argument has a column of coefficients for every constraint, and the four argument has the lower bound. Here we’re saying that the first argument must be at least -100, which seems safe enough. The constrained solution is the same as the unconstrained solution, because the unconstrained solution meets the constraint.

Since we are maximizing, the value at the optimum is the Intercept from the model minus the value reported here, in this case, 61.1- (-6.1) = 67.2.

> library(quadprog)
> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,0,0,0),ncol=1),-100)
$solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896

$value
[1] -6.123755

$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896

$iterations
[1] 1 0

$Lagrangian
[1] 0

$iact
[1] 0

OK, let’s make it bite a little bit. Let’s constrain the third variable to be at least 1, and the first variable to be no more than one (its negative no more than -1).

> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(0,0,1,0,-1,0,0,0),ncol=2),c(1,-1))
$solution
[1] 1.0000000 1.2860982 1.0000000 0.5248153

$value
[1] -6.05187

$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896

$iterations
[1] 3 0

$Lagrangian
[1] 0.5113095 0.3332150

$iact
[1] 2 1

The first meq of the constraints are equality constraints. It defaults to zero, but here we say that the first variable must equal 1.1.

> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,0,0,0),ncol=1),1.1,meq=1)
$solution
[1] 1.1000000 1.2860982 0.8834306 0.5058415

$value
[1] -6.108397

$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896

$iterations
[1] 2 0

$Lagrangian
[1] 0.201295

$iact
[1] 1

You can use more interesting constraints like the first variable minus the second variable must be at least .5.

> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,-1,0,0),ncol=1),.5)
$solution
[1] 1.5986920 1.0986920 0.8834306 0.4112207

$value
[1] -6.00196

$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896

$iterations
[1] 2 0

$Lagrangian
[1] 0.4565795

$iact
[1] 1