Resources

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;