options nodate nonumber ps=150 ls=80;

title 'RPD -- Example 8.5 -- Salmon Pentachlorophenate Resistance';

data salmon;
input rep salinity temperature oxygen lethal_time;
datalines;
1    -1       -1       -1         53 
1    -1       -1        1         54 
1    -1        1       -1         40 
1    -1        1        1         37 
1     1       -1       -1         84 
1     1       -1        1         76 
1     1        1       -1         40 
1     1        1        1         50 
1     0        0        0         50 
1   1.215      0        0         61
1  -1.215      0        0         54
1     0      1.215      0         39
1     0     -1.215      0         67
1     0        0      1.215       44
1     0        0     -1.215       61
2    -1       -1       -1         50
2    -1       -1        1         42
2    -1        1       -1         31
2    -1        1        1         28
2     1       -1       -1         57
2     1       -1        1         78
2     1        1       -1         49
2     1        1        1         54
2     0        0        0         50
2   1.215      0        0         76
2  -1.215      0        0         45
2     0      1.215      0         33
2     0     -1.215      0         54
2     0        0      1.215       45
2     0        0     -1.215       38
3  -1.2500  -1.8867  -0.6350      46
3   0.8600  -2.2200  -0.4250      66
3   1.0000  -2.2400  -0.3100      68
3   2.1165  -2.4167  -0.1450      75
3   2.5825  -2.4900  -0.0800      75
3   3.2475  -2.6667   0.0800      68
3   1.1760  -1.3333     0         78
3   1.4700  -1.6667     0         93
3   1.7640  -2.0000     0         96
3   2.0580  -2.3333     0         66
;


proc iml;

use salmon;
read all;

Y = lethal_time;

   /* Quadratic polynomial response surface fit,
      allowing different intercepts for different reps: */

X = design(rep)||salinity||temperature||oxygen||
                 salinity##2||temperature##2||oxygen##2||
                 salinity#temperature||salinity#oxygen||temperature#oxygen;

   /* Note: design(rep) implicitly includes the intercept */

print X;

XTXINV = inv(X`*X);

BETAHAT = XTXINV*(X`*Y);

SSRES = Y`*Y - BETAHAT`*X`*Y;

S2 = SSRES/(nrow(X)-ncol(X));

SE_BETAHAT = sqrt(vecdiag(XTXINV)#S2);

T_STATS = BETAHAT/SE_BETAHAT;

print BETAHAT SE_BETAHAT T_STATS S2;


   /* After performing selection of terms, as discussed in the textbook: */

X_R = design(rep)||salinity||temperature||temperature##2;

XTXINV_R = inv(X_R`*X_R);

BETAHAT_R = XTXINV_R*(X_R`*Y);

SSRES_R = Y`*Y - BETAHAT_R`*X_R`*Y;

S2_R = SSRES_R/(nrow(X_R)-ncol(X_R));

SE_BETAHAT_R = sqrt(vecdiag(XTXINV_R)#S2_R);

print BETAHAT_R SE_BETAHAT_R S2_R;

quit;


/* PROC RSREG is specialized to fit quadratic response surfaces.  It cannot
   directly handle class variables, so indicator variables for rep must be 
   added to the data set: */

data salmon;
set salmon;
rep1 = (rep=1);
rep2 = (rep=2);
rep3 = (rep=3);
run;


proc rsreg;
model lethal_time = rep1 rep2 rep3 salinity temperature oxygen / 
                    covar=3 nocode nooptimal;
   /* 'covar=3':    treat first 3 independent variables as "covariates"
                      (so they do not generate quadratic terms)
      'nocode':     suppress recoding of variables
      'nooptimal':  suppress display of canonical analysis */
run;

   /* No easy way to reduce the model within PROC RSREG. */