Response Surfaces
Goal: Fit full quadratic model in N continuous variables
N=3:
1 intercept + N linear + N pure quadratic +
N(N-1)/2 crossproduct terms =(N+1)(N+2)/2 parameters.
Central Composite Design
2N factorial + 2N axial points + m center points
Assume coded factorial coordinates -1, 1
Example: N=3
8 factorial points (-1,-1,-1),(-1,-1,1),…,(1,1,1)
6 axial points (,0,0),(,0,0),(0,,0),…(0,0,).
m center points.
Example N=2 X=
Surface:
Clearly (2,3) is the critical (stationary) point and the predicted response there is 25. Is that a maximum, minimum or saddle point?
Canonical Analysis:
(1) Get eigenvalues:
A = Z Z’
Columns of Z are “eigenvectors”
Let D’=(row) vector of deviations=(x1-2,x2-3)
= c + D’AD = c + D’ZZ’D=W’W where this defines a new set of axes, “canonical axes” by W’=(w1,w2)=D’Z
[Note that Z’Z=ZZ’=I so that W=Z’DZW=D]
Specifically we have and
Finally, we have , a saddle point. The critical point is (0,0) and if we move from there along the w1 axis, the predictions decrease while movement along the w2 axis causes increasing predictions.
**Quadsurface.sas**
Illustrates the ideas of eigenvalues and eigenvectors, giving 3 illustrations of the surfaces:
Minimum Maximum Saddle Point
Note: (X1,X2) axis, black square (rug) shows (W1,W2) axes.
Contour plot axes: (W1,W2)
**DOE10.SAS**
Runs this analysis in PROC RSREG, including graphics
Example: Chicken Broth Flavor as judged by 12 judges (twice each). Y=average rating.
Control variables: MSG, NaCl amounts.
** Run DOE11.sas **
Ridge Analysis
Several of our examples have explored the response surface through something called “ridge analysis” not to be confused with ridge regression. Imagine a probe point on the surface. Now in the floor, draw concentric circles with that point as center and increasing radii. For each circle, look at the function values and find the coordinates that maximize (minimize) the predictions. That is the idea of ridge analysis.
Ridge analysis uses a mathematical optimization technique called the LaGrange multiplier technique.
Example 1: Try to maximize Y=10-5x +2x2
subject to the restriction (X-1)2=4
Solution: x is -1 (Y=17) or 3 (Y=13) so x=-1 is answer.
Example 2: Try to maximize
Y=25 +(x1-2)2 + 2(x2-3)2+2(x1-2)(x2-3) for (x1,x2) on a circle with radius r and center (5,3).
Solution: Find critical points of
f(x1,x2,)=25 +(x1-2)2 + 2(x2-3)2+2(x1-2)(x2-3)-
[(x1-5)2 + (x2-3)2-r2].
(LaGrange multiplier)
(1) f(x1,x2,)=(x1-5)2 + (x2-3)2-r2
Setting this to 0 forces solution on the circle.
(2) f(x1,x2,)=constant+linear +
= constant + linear +
This is quadratic for all values and for any , this quadratic surface must have a critical point (allnontrivial quadratic surfaces do). As before, find the eigenvalues or “roots”
A =
Roots found using |A-I|=0. (1-)(2-)-6=
2-l(3-2)+(1-)(2-)-4 =(+1+)(-4+)=0
Critical point is maximum if >4.
Critical point is minimum if <-1.
Otherwise, it is a saddle point.
On circle of radius r, (x1-5)2 + (x2-3)2-r2 =0 so
Y=f(x1,x2,) on circle of radius r. Set all derivatives of f(x1,x2,) to 0.
(x1-5)2 + (x2-3)2-r2=0
2(x1-2) +2(x2-3)-2(x1-5)=0
4(x2-3)+2(x1-2) -4 (x2-3)=0
For some M>4, f(x1,x2,) must have a global maximum and for some m1 a minimum. The global max and min for f(x1,x2,) must be on the circle (why?) and thus must be a max or min over the circle. Intuitively, the surface above the circle is like a necklace or rubber band lying on the surface which, by its connected nature must have a maximum and minimum.
** DOE12.sas** gives the picture:
Most analyses, such as that in PROC RSREG uses circle in the coded data (it makes a difference!)
Example 3: Conversion of 1,2 propanediol to 2,5 dimethylpiperaznie (Meyers: (1971 pg. 85)
Pct is % conversion.
** DOE13.sas **
Example 4:A more complex example: 8 factors, central composite design. Box: ¼ replicate of 28, run in 4 blocks of size 16 plus center points. Axial points (16) and center points in block 5.
** DOE14.sas**
RSREG – shows only 2 active factors
RSREG – shows 8 factor surface stationary point is saddle point.
RSREG on 2 active factors (x1, x6) has maximum.
Example 5: Eleven design factors affecting noise in submarines. Run a 211 first to identify active factors. Add axial points for active factors in a second block.
** DOE15.sas**
DATA LIN; BLOCK=1;
INPUT A B C D E F G H I J Y CENTER;
LABEL A="SHAFT DIAMETER";
LABEL B="SHAFT METAL";
LABEL C="SHAFT LENGTH";
LABEL D="PROP PITCH";
LABEL E="PROP METAL";
LABEL F="PROP BLADE THICKNESS";
LABEL G="PROP BLADE WIDTH";
LABEL H="PROP BLADE LENGTH";
LABEL I="RUDDER HEIGHT";
LABEL J="RUDDER WIDTH";
LABEL Y= "NOISE";
DATALINES;
(data)
PROCRSREGdata=all;
MODEL Y = block A D H/covar=1lackfit;
covar=1: first variable comes in linearly.
lackfit: test if more than quadratic terms are needed (need pure replicates at some points, e.g. center).
Three factors (A, D, H) are active, critical point is a minimum but far from reasonable settings - run ridge in reasonable settings area:
RIDGECENTER=(0,0,0) radius=0 to 3 by .1outr=quiet
minimum; run;