Output from imlreg1.sas	Source 0 Graphs

NOTE: Capture of log output started.

NOTE: %INCLUDE (level 1) file n:\psy6140\examples\iml\imlreg1.sas is file
      n:\psy6140\examples\iml\imlreg1.sas.

972  + /****************************************************************/
973  + /*    NAME: IMLREG1                                             */
974  + /*   TITLE: Regression with IML                                 */
975  + /*    DATA: THERAPY data                                        */
976  + /****************************************************************/
977  +
978  +data THERAPY;
979  +     input NAME $    SEX $  PERSTEST  THERAPY  INTEXT;
980  +
981  +     If SEX = 'F' then SX = 1;    /* Dummy variable for sex */
982  +                  else SX = 0;
983  +cards;

NOTE: The data set WORK.THERAPY has 10 observations and 6 variables.
NOTE: The DATA statement used 0.22 seconds.

994  +;
995  +proc iml;

IML Ready

996  +
997  + /*---------REGRESSION ROUTINE------------*/
998  + /* given X, and y, this fits y = X b + e */
999  + /* by least squares.                     */
1000 +
1001 +start reg;
1002 +  n = nrow(X);
1002 +                                /* number of observations   */
1003 +  k = ncol(X);
1003 +                                /* number of variables      */
1004 +  xpx = X`*X;
1004 +                                /* cross-products, X'X      */
1005 +  xpy = X`*y;
1006 +  xpxi = inv(xpx);
1006 +                                /* inverse crossproducts    */
1007 +  b = xpxi*xpy;
1007 +                                /* parameter estimates      */
1008 +
1009 +  yhat = X*b;
1009 +                                /* predicted values         */
1010 +  resid = y-yhat;
1010 +                                /* residuals                */
1011 +  sse = resid`*resid;
1011 +                                /* sum of squared errors    */
1012 +  dfe = n-k;
1012 +                                /* degrees of freedom error */
1013 +  mse = sse/dfe;
1013 +                                /* mean squared error       */
1014 +  rmse = sqrt(mse);
1014 +                                /* root mean squared error  */

1016 +
1017 +  covb = xpxi#mse;
1017 +                                      /* covariance of estimates  */
1018 +  stdb = sqrt(vecdiag(covb));
1018 +                                      /* standard errors          */
1019 +  t = b/stdb;
1019 +                                      /* ttest for estimates=0    */
1020 +  probt = 1-probf(t#t,1,dfe);
1020 +                                      /* significance probability */
1021 +  print vnames b stdb t[f=7.3]
1022 +                 probt[f=7.4];
1023 +  s = diag(1/stdb);
1024 +  corrb = s*covb*s;
1024 +                                      /* correlation of estimates */
1025 +  print "Covariance of Estimates",  covb[r=vnames c=vnames],
1026 +        "Correlation of Estimates",corrb[r=vnames c=vnames f=7.3];
1027 +
1028 +  if nrow(alpha)= 0 then return;
1028 +                                      /* is alpha specified?      */
1029 +  H = X*xpxi*X`;
1029 +                                      /* hat matrix               */
1030 +  vresid = (I(n)-H)*mse;
1030 +                                      /* covariance of residuals  */
1031 +  vpred  = H#mse;
1031 +                                     /* covariance of yhat values */
1032 +  hat= vecdiag(H);
1032 +                                      /* hat leverage values      */
1033 +
1034 +  alpha=.05;
1035 +  t = tinv(1-alpha/2,dfe);
1036 +  lowerm=yhat-t#sqrt(hat*mse);
1036 +                                    /* lower conf. limit for mean */
1037 +  upperm=yhat+t#sqrt(hat*mse);
1037 +                                    /* upper */
1038 +  lower =yhat-t#sqrt(hat*mse+mse);
1038 +                                    /* lower conf. limit for indiv */
1039 +  upper =yhat+t#sqrt(hat*mse+mse);
1039 +                                    /* upper */
1040 +  print "Predicted Values, Residuals, Leverage, and Conf. Limits" ,
1041 +        id y yhat resid[f=6.2] hat[f=6.2]
1042 +        lowerm[f=6.2] upperm[f=6.2]
1043 +    /*  lower[f=6.2]  upper[f=6.2] */ ;
1044 +finish;

NOTE: Module REG defined.

1045 +
1046 +
1047 + /*---routine to test a linear combination of the estimates---*/
1048 + /* given L, this routine tests hypothesis that L b = 0.      */
1049 +
1050 +start test;
1051 +  dfh=nrow(L);
1051 +                                  /* degrees of freedom of test */
1052 +  Lb=L*b;
1052 +                                  /* lin. comb of coeffs tested */
1053 +  vlb=L*xpxi*L`;
1053 +                                  /* variance of LB             */
1054 +  q=Lb`*inv(vlb)*lb /dfh;
1054 +                                  /* Mean square for H0: LB=0   */
1055 +  f=q/mse;
1056 +  prob=1-probf(f,dfh,dfe);
1057 +  print  l[c=vnames], lb f dfh dfe prob;
1058 +finish;

NOTE: Module TEST defined.

1059 +
1060 +  /* Read the SAS data set THERAPY into IML variables */
1061 +  use therapy;
1062 +  read all var{THERAPY} into y;
1063 +  read all var{PERSTEST INTEXT SX} into X[r=NAME c=VNAMES];
1064 +
1065 +  X = J( nrow(X),1) || X;
1065 +                                  /* add column of 1s           */
1066 +  vnames = 'Intercept' // vnames`;
1067 +  id = name;
1068 +
1069 +  alpha=.05;
1069 +                         /* Error rate for confidence intervals */
1070 +  run reg;
1071 +
1072 +  print "Test Coef for X1";
1072 +                                L= {0 1 0 0};
1072 +                                                        run test;
1073 +* print "TEST Coef for X1, X2";
1073 +                                L= {0 1 0 0,0 0 1 0} ;
1073 +                                                       *run test;
1074 +quit;

Exiting IML.

NOTE: The PROCEDURE IML used 0.59 seconds.

1075 +

NOTE: %INCLUDE (level 1) ending.