High Resolution Plot of Response Surface and Ridge
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: RSREGPLT */
/* TITLE: High Resolution Plot of Response Surface and Ridge */
/* PRODUCT: STAT GRAPH */
/* SYSTEM: ALL */
/* KEYS: response surface methods, graphics */
/* PROCS: RSREG GCONTOUR */
/* DATA: Frankel (1961), reported in Myers (1976), p. 158. */
/* */
/* SUPPORT: RDT UPDATE: January 2010 */
/* REF: Myers, Raymond H. (1976), "Response Surface Metho- */
/* dology", Blacksburg, Virginia: Virginia Poly- */
/* technic Institute and State University. */
/* MISC: */
/* */
/****************************************************************/
/*
/ The data is from a rotatable central composite design with two
/ factor variables: the object is to find the settings of process
/ time and temperature which maximize the yield of a chemical.
/---------------------------------------------------------------*/
data design;
length Position $ 1;
input Time Temp MBT Position;
label time = "Reaction Time (Hours)"
temp = "Temperature (Degrees Centigrade)"
mbt = "Percent Yield Mercaptobenzothiazole";
datalines;
4.0 250 83.8 3
4.0 250 82.0 9
6.3 229 81.3 2
6.3 271 83.1 2
12.0 220 84.7 2
12.0 250 82.4 1
12.0 250 82.9 7
12.0 250 81.2 6
12.0 280 57.9 8
17.7 229 85.3 8
17.7 271 72.7 2
20.0 250 81.7 2
;
/*
/ We'll use RSREG to compute the predicted surface directly by
/ appending to the original data a grid of points over the
/ region of interest (with MBT set to missing so these points
/ don't get involved in the fit.)
/---------------------------------------------------------------*/
data desnsurf;
set design end=eof; /* The original data */
output;
if eof then do;
mbt = . ; /* Missing response flags no fit */
do time=3 to 21 by 1; /* a grid for time and temp */
do temp=215 to 285 by 5;
output;
end;
end;
end;
run;
/*
/ Fit the full quadratic response surface, outputting the
/ predicted surface over the grid in FIT and the ridge of
/ maximum predicted yield in RIDGE.
/---------------------------------------------------------------*/
ods graphics on;
proc rsreg data=desnsurf out=fit plots=surface;
model mbt = time temp / lackfit predict;
ridge max outr=ridge;
run;
ods graphics off;
/*
/ Form data sets for plotting: the predicted surface.
/ DESNRDG - The DATA 'n the RIDGE.
/---------------------------------------------------------------*/
data surface;
set fit;
if _n_ > 12;
run;
/*
/ An ANNOTATE= data set to plot the design points.
/---------------------------------------------------------------*/
data annodesn;
set design;
length text $ 4;
length function $ 8;
x = time; y = temp;
xsys = '2' ; ysys = '2' ;
size = 3.00;
function = 'label'; text = left(mbt); /* the value of mbt */
if (_n_ = 8) then x = x + .5;
output;
if (_n_ = 8) then x = x - .5;
function = 'symbol'; text = '+'; /* the design point */
position = '5';
output;
run;
/*
/ An ANNOTATE= data set to plot the ridge of maximum response.
/---------------------------------------------------------------*/
data annordg;
/*
/ Variables we won't need after this.
/-------------------------------------------------------------*/
drop _depvar_ _type_ _radius_ _pred_ _stderr_ pi twopi thisx
thisy dx dy d t t1 t2 lastx lasty;
/*
/ Constants.
/-------------------------------------------------------------*/
pi = constant('pi'); twopi = 2*pi;
/*
/ All through, we interpret X and Y in terms of data values.
/-------------------------------------------------------------*/
xsys = '2'; ysys = '2';
do n = 1 to 11;
set ridge point = n;
x = time; y = temp;
if n = 1 then do;
/*
/ Move to the first point.
/-------------------------------------------------------*/
function = 'move';
output;
lastx = x; lasty = y;
end;
/*
/ Draw a line from the last point to this point.
/----------------------------------------------------------*/
else do;
function = 'draw';
output;
/*
/ Now, a little code to put an arrow-head at the end of
/ the last line drawn. First, preserve the coordinates
/ of the present point.
/-------------------------------------------------------*/
thisx = x; thisy = y;
/*
/ Compute the (scaled) x-distance, the (scaled) y-distance,
/ and the overall (scaled) distance since the last point.
/-------------------------------------------------------*/
dx = (lastx - thisx)/8;
dy = (lasty - thisy)/30;
d = sqrt(dx**2 + dy**2);
/*
/ The angle of the line drawn from the last point.
/-------------------------------------------------------*/
t = arcos(dx/d);
if (arsin(dy/d) <= 0) then t = pi + (pi - t);
/*
/ A point 2/5 of the way back toward the last one and at
/ a slightly smaller angle: draw a line out to it, then
/ move back to the present point.
/-------------------------------------------------------*/
t1 = t - .25;
if (t1 < 0) then t1 = t1 + twopi;
x = thisx + 8*(.4*d)*cos(t1);
y = thisy + 30*(.4*d)*sin(t1);
output;
function = 'move';
x = thisx;
y = thisy;
output;
/*
/ A point 2/5 of the way back toward the last one and at
/ a slightly greater angle: once again, draw a line out
/ to it, then move back to the present point.
/-------------------------------------------------------*/
function = 'draw';
t2 = t + .25;
if (t1 > twopi) then t1 = t1 - twopi;
x = thisx + 8*(.4*d)*cos(t2);
y = thisy + 30*(.4*d)*sin(t2);
output;
function = 'move';
x = thisx;
y = thisy;
output;
/*
/ Finally, remember this point for the next time around.
/-------------------------------------------------------*/
lastx = x; lasty = y;
end;
end;
stop;
run;
data anno;
set annodesn annordg;
run;
/*
/ Plot the surface with a contour plot, adding the design and
/ the ridge with an ANNOTATE = option.
/---------------------------------------------------------------*/
proc gcontour data=surface;
plot temp * time = mbt / annotate = anno levels = 54 to 90 by 3;
run;