/****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: REGDEMO2                                            */
 /*   TITLE: REGRESSION OF SUBSETS OF VARIABLES                  */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX  REGR    SUGI6                               */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: LWB                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/

proc iml;

 /*-------INITIALIZATION-------------------------------*
 | C,CSAVE THE CROSSPRODUCTS MATRIX                   |
 | N       NUMBER OF OBSERVATIONS                     |
 | K       TOTAL NUMBER OF VARIABLES TO CONSIDER      |
 | L       NUMBER OF VARIABLES CURRENTLY IN MODEL     |
 | IN      A 0-1 VECTOR OF WHETHER VARIABLE IS IN     |
 | B       PRINT COLLECTS RESULTS (L MSE RSQ BETAS )  |
 *----------------------------------------------------*/
START INITIAL;
  N=NROW(X); K=NCOL(X); K1=K+1; IK=1:K;
  BNAMES={NPARM MSE RSQUARE} ||VARNAMES;

  /*---CORRECT BY MEAN, ADJUST OUT INTERCEPT PARAMETER---*/
  Y=Y-Y[+,]/N;                   /* CORRECT Y BY MEAN */
  X=X-REPEAT(X[+,]/N,N,1);       /* CORRECT X BY MEAN */
  XPY=X`*Y;                      /* CROSSPRODUCTS     */
  YPY=Y`*Y;
  XPX=X`*X;
  FREE X Y;                      /* NO LONGER NEED THE DATA*/
  CSAVE=(XPX || XPY) //
        (XPY`|| YPY);            /* SAVED COPY OF CROSSPRODUCTS*/
FINISH;


 /*-----FORWARD METHOD------------------------------------------*/
START FORWARD;
  PRINT / "FORWARD SELECTION METHOD";
  FREE BPRINT;
  C=CSAVE; IN=REPEAT(0,K,1); L=0;      /* NO VARIABLES ARE IN */
  DFE=N-1; MSE=YPY/DFE;
  SPROB=0;

  DO WHILE(SPROB<.15 & L<K);
     INDX=LOC(^IN);           /* WHERE ARE THE VARIABLES NOT IN?*/
     CD=VECDIAG(C)[INDX,];    /* XPX DIAGONALS                  */
     CB=C[INDX,K1];           /* ADJUSTED XPY                   */
     TSQR=CB#CB/(CD#MSE);     /* SQUARES OF T TESTS             */
     IMAX=TSQR[<:>,];         /* LOCATION OF MAXIMUM IN INDX    */
     SPROB=(1-PROBT(SQRT(TSQR[IMAX,]),DFE))*2;
     IF SPROB<.15 THEN DO;    /* IF T-TEST SIGNIFICANT          */
        II=INDX[,IMAX];       /* PICK MOST SIGNIFICANT          */
        RUN SWP;              /* ROUTINE TO SWEEP               */
        RUN BPR;              /* ROUTINE TO COLLECT RESULTS     */
        END;
     END;

  PRINT  BPRINT[COLNAME=BNAMES] ;
FINISH;



 /*-----BACKWARD METHOD----------------------------------------*/
START BACKWARD;
  PRINT / "BACKWARD ELIMINATION ";
  FREE BPRINT;
  C=CSAVE; IN=REPEAT(0,K,1);
  II=1:K; RUN SWP; RUN BPR;       /* START WITH ALL VARIABLES IN*/
  SPROB=1;

  DO WHILE(SPROB>.15 & L>0);
     INDX=LOC(IN);                /* WHERE ARE THE VARIABLES IN? */
     CD=VECDIAG(C)[INDX,];        /* XPX DIAGONALS               */
     CB=C[INDX,K1];               /* BVALUES                     */
     TSQR=CB#CB/(CD#MSE);         /* SQUARES OF T TESTS          */
     IMIN=TSQR[>:<,];             /* LOCATION OF MINIMUM IN INDX */
     SPROB=(1-PROBT(SQRT(TSQR[IMIN,]),DFE))*2;
     IF SPROB>.15 THEN DO;        /* IF T-TEST NONSIGNIFICANT    */
        II=INDX[,IMIN];           /* PICK LEAST SIGNIFICANT      */
        RUN SWP;                  /* ROUTINE TO SWEEP IN VARIABLE*/
        RUN BPR;                  /* ROUTINE TO COLLECT RESULTS  */
        END;
     END;

  PRINT  BPRINT[COLNAME=BNAMES] ;
FINISH;



 /*-----STEPWISE METHOD-----------------------------------------*/
START STEPWISE;
  PRINT /"STEPWISE METHOD";
  FREE BPRINT;
  C=CSAVE; IN=REPEAT(0,K,1); L=0;
  DFE=N-1; MSE=YPY/DFE;
  SPROB=0;

  DO WHILE(SPROB<.15 & L<K);
     INDX=LOC(^IN);           /* WHERE ARE THE VARIABLES NOT IN?*/
     NINDX=LOC(IN);           /* WHERE ARE THE VARIABLES IN?    */
     CD=VECDIAG(C)[INDX,];    /* XPX DIAGONALS                 */
     CB=C[INDX,K1];           /* ADJUSTED XPY                  */
     TSQR=CB#CB/CD/MSE;       /* SQUARES OF T TESTS            */
     IMAX=TSQR[<:>,];         /* LOCATION OF MAXIMUM IN INDX   */
     SPROB=(1-PROBT(SQRT(TSQR[IMAX,]),DFE))*2;
     IF SPROB<.15 THEN DO;    /* IF T-TEST SIGNIFICANT         */
        II=INDX[,IMAX];       /* FIND INDEX INTO C             */
        RUN SWP;              /* ROUTINE TO SWEEP              */
        RUN BACKSTEP;         /* CHECK IF REMOVE ANY TERMS     */
        RUN BPR;              /* ROUTINE TO COLLECT RESULTS    */
        END;
     END;

  PRINT  BPRINT[COLNAME=BNAMES] ;
FINISH;


 /*----ROUTINE TO BACKWARDS-ELIMINATE FOR STEPWISE--*/
START BACKSTEP;
  IF NROW(NINDX)=0 THEN RETURN;
  BPROB=1;
  DO WHILE(BPROB>.15 & L<K);
     CD=VECDIAG(C)[NINDX,];      /* XPX DIAGONALS         */
     CB=C[NINDX,K1];             /* BVALUES               */
     TSQR=CB#CB/(CD#MSE);        /* SQUARES OF T TESTS    */
     IMIN=TSQR[>:<,];            /* LOCATION OF MINIMUM IN NINDX*/
     BPROB=(1-PROBT(SQRT(TSQR[IMIN,]),DFE))*2;
     IF BPROB>.15 THEN DO;
        II=NINDX[,IMIN];
        RUN SWP;
        RUN BPR;
        END;
     END;
FINISH;


 /*-----SEARCH ALL POSSIBLE MODELS----------------------------*/
START ALL;
     /*---USE METHOD OF SCHATZOFF ET AL. FOR SEARCH TECHNIQUE--*/
  BETAK=REPEAT(0,K,K); /* RECORD ESTIMATES FOR BEST L-PARAM MODEL*/
  MSEK=REPEAT(1E50,K,1);/* RECORD BEST MSE PER # PARMS            */
  RSQK=REPEAT(0,K,1);  /* RECORD BEST RSQUARE                    */
  INK=REPEAT(0,K,K);   /* RECORD BEST SET PER # PARMS            */
  LIMIT=2##K-1;        /* NUMBER OF MODELS TO EXAMINE            */

C=CSAVE; IN=REPEAT(0,K,1);/* START OUT WITH NO VARIABLES IN MODEL*/

     DO KK=1 TO LIMIT;
        RUN ZTRAIL;              /* FIND WHICH ONE TO SWEEP    */
        RUN SWP;                 /* SWEEP IT IN                */
        BB=BB//(L||MSE||RSQ||(C[IK,K1]#IN)`);
        IF MSE<MSEK[L,] THEN DO; /* WAS THIS BEST FOR L PARMS? */
           MSEK[L,]=MSE;         /* RECORD MSE                 */
           RSQK[L,]=RSQ;         /* RECORD RSQUARE             */
           INK[,L]=IN;           /* RECORD WHICH PARMS IN MODEL*/
           BETAK[L,]=(C[IK,K1]#IN)`;/* RECORD ESTIMATES        */
           END;
        END;

     PRINT / "ALL POSSIBLE MODELS" " IN SEARCH ORDER";
     PRINT BB[COLNAME=BNAMES]; FREE BB;
     BPRINT=IK`||MSEK||RSQK||BETAK;
     PRINT ,"THE BEST MODEL FOR EACH NUMBER OF PARAMETERS";
     PRINT BPRINT[COLNAME=BNAMES];

     /*---MALLOWS CP PLOT---*/
     CP=MSEK#(N-IK`-1)/MIN(MSEK)-(N-2#IK`);
     CP=IK`||CP;
     CPNAME={"NPARM" "CP"};
     /*OUTPUT CP OUT=CP COLNAME=CPNAME;*/
FINISH;



 /*-SUBROUTINE TO FIND NUMBER OF TRAILING ZEROS IN BINARY NUMBER*/
 /* ON ENTRY: KK IS THE NUMBER TO EXAMINE                       */
 /* ON EXIT:  II HAS THE RESULT                                 */
 /*-------------------------------------------------------------*/
START ZTRAIL;
    II=1; ZZ=KK;
    DO WHILE(MOD(ZZ,2)=0); II=II+1; ZZ=ZZ/2; END;
FINISH;


 /*-----SUBROUTINE TO SWEEP IN A PIVOT--------------------------*/
 /* ON ENTRY: II HAS THE POSITION(S) TO PIVOT                   */
 /* ON EXIT:  IN, L, DFE, MSE, RSQ RECALCULATED                 */
 /*-------------------------------------------------------------*/
START SWP;
     IF ABS(C[II,II])<1E-9 THEN DO; PRINT , "FAILURE", C;STOP;END;
     C=SWEEP(C,II);
     IN[II,]=^IN[II,];
     L=SUM(IN); DFE=N-1-L;
     SSE=C[K1,K1];
     MSE=SSE/DFE;
     RSQ=1-SSE/YPY;
FINISH;


 /*-----SUBROUTINE TO COLLECT BPRINT RESULTS--------------------*/
 /* ON ENTRY: L,MSE,RSQ, AND C SET UP TO COLLECT                */
 /* ON EXIT:  BPRINT HAS ANOTHER ROW                            */
 /*-------------------------------------------------------------*/
START BPR;
     BPRINT=BPRINT//(L||MSE||RSQ||(C[IK,K1]#IN)`);
FINISH;



 /*--------------STEPWISE METHODS---------------------*/
 /* AFTER A RUN TO THE INITIAL ROUTINE, WHICH SETS UP */
 /* THE DATA, FOUR DIFFERENT ROUTINES CAN BE CALLED   */
 /* TO DO FOUR DIFFERENT MODEL-SELECTION METHODS.     */
 /*---------------------------------------------------*/
START SEQ;
  RUN INITIAL;          /* INITIALIZATION             */
  RUN ALL;              /* ALL POSSIBLE MODELS        */
  RUN FORWARD;          /* FOREWARD SELECTION METHOD  */
  RUN BACKWARD;         /* BACKWARD ELIMINATION METHOD*/
  RUN STEPWISE;         /* STEPWISE METHOD            */
  FINISH;


 /*------------------------DATA ON PHYSICAL FITNESS--------------*
| THESE MEASUREMENTS WERE MADE ON MEN INVOLVED IN A PHYSICAL     |
| FITNESS COURSE AT N.C.STATE UNIV.  THE VARIABLES ARE AGE(YEARS)|
| WEIGHT(KG), OXYGEN UPTAKE RATE(ML PER KG BODY WEIGHT PER       |
| MINUTE), TIME TO RUN 1.5 MILES(MINUTES), HEART RATE WHILE      |
| RESTING, HEART RATE WHILE RUNNING (SAME TIME OXYGEN RATE       |
| MEASURED), AND MAXIMUM HEART RATE RECORDED WHILE RUNNING.      |
| CERTAIN VALUES OF MAXPULSE WERE MODIFIED FOR CONSISTENCY.      |
|   DATA COURTESY DR. A.C. LINNERUD                              |
*---------------------------------------------------------------*/
data =
   { 44 89.47  44.609 11.37 62 178 182  ,
     40 75.07  45.313 10.07 62 185 185  ,
     44 85.84  54.297  8.65 45 156 168  ,
     42 68.15  59.571  8.17 40 166 172  ,
     38 89.02  49.874  9.22 55 178 180  ,
     47 77.45  44.811 11.63 58 176 176  ,
     40 75.98  45.681 11.95 70 176 180  ,
     43 81.19  49.091 10.85 64 162 170  ,
     44 81.42  39.442 13.08 63 174 176  ,
     38 81.87  60.055  8.63 48 170 186  ,
     44 73.03  50.541 10.13 45 168 168  ,
     45 87.66  37.388 14.03 56 186 192  ,
     45 66.45  44.754 11.12 51 176 176  ,
     47 79.15  47.273 10.60 47 162 164  ,
     54 83.12  51.855 10.33 50 166 170  ,
     49 81.42  49.156  8.95 44 180 185  ,
     51 69.63  40.836 10.95 57 168 172  ,
     51 77.91  46.672 10.00 48 162 168  ,
     48 91.63  46.774 10.25 48 162 164  ,
     49 73.37  50.388 10.08 67 168 168  ,
     57 73.37  39.407 12.63 58 174 176  ,
     54 79.38  46.080 11.17 62 156 165  ,
     52 76.32  45.441  9.63 48 164 166  ,
     50 70.87  54.625  8.92 48 146 155  ,
     51 67.25  45.118 11.08 48 172 172  ,
     54 91.63  39.203 12.88 44 168 172  ,
     51 73.71  45.790 10.47 59 186 188  ,
     57 59.08  50.545  9.93 49 148 155  ,
     49 76.32  48.673  9.40 56 186 188  ,
     48 61.24  47.920 11.50 52 170 176  ,
     52 82.78  47.467 10.50 53 170 172  };


x=data[,{1 2 4 5 6 7 }];
y=data[,3];
free data;
varnames={age weight runtime rstpulse runpulse maxpulse};
reset fw=8 linesize=90;
run seq;