Quadratic Response Surface Regression

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: RSM                                                 */
 /*   TITLE: Quadratic Response Surface Regression               */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX  REGR    SUGI6                               */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: RHD                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/


proc iml;

 /*----------QUADRATIC RESPONSE SURFACE REGRESSION---------------*
  | This matrix routine reads in the factor variables and the    |
  | response, forms the quadratic regression model and estimates |
  | the parameters, then solves for the optimal response, prints |
  | the optimal factors and response, and then displays the      |
  | the eigenvalues and eigenvectors of the matrix of quadratic  |
  | parameter estimates to determine if the solution is a maximum|
  | or minimum, or saddlepoint, and which direction has the      |
  | steepest and gentlest slopes.                                |
  *--------------------------------------------------------------*
  | Given that D contains the factor variables,                  |
  |        and Y contains the response.                          |
  *--------------------------------------------------------------*/
START RSM;
  N=NROW(D);  K=NCOL(D);               /* dimensions */
  X=J(N,1,1)||D;                       /* set up design matrix */
  DO I=1 TO K; DO J=1 TO I; X=X||D[,I]#D[,J]; END; END;

  BETA=SOLVE(X`*X,X`*Y);          /* solve parameter estimates */
  PRINT "Parameter Estimates" , BETA;

  C=BETA[1];                      /* intercept estimate */
  B=BETA[2:(K+1)];               /* linear estiamtes */
  A=J(K,K,0); L=K+1;             /* form quadratics into matrix */
  DO I=1 TO K; DO J=1 TO I; L=L+1; A[I,J]=BETA[L,]; END; END;
  A=(A+A`)*.5;                   /* symmetrize */

  XX=-.5*SOLVE(A,B);             /* solve for critical value */
  PRINT , "Critical Factor Values" , XX;
  YOPT=C + B`*XX + XX`*A*XX; /*compute response at critical value*/
  PRINT , "Response at Critical Value" YOPT;

  CALL EIGEN(EVAL,EVEC,A);
  PRINT , "Eigenvalues and Eigenvectors", EVAL, EVEC;

  IF MIN(EVAL)>0 THEN PRINT , "Solution Was a Minimum";
  IF MAX(EVAL)<0 THEN PRINT , "Solution Was a Maximum";
FINISH;




 /* Sample Problem with Two Factors */

D={-1 -1,
   -1  0,
   -1  1,
    0 -1,
    0  0,
    0  1,
    1 -1,
    1  0,
    1  1};

Y={ 71.7, 75.2, 76.3, 79.2, 81.5, 80.2, 80.1, 79.1, 75.8};

RUN RSM;