options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 4.10 -- Physical Fitness Data';

data physfit;
input age weight oxygen time hr_rest hr_run hr_max;
datalines;
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 184
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 176
52 76.32 45.441  9.63 48 164 166
50 70.87 54.625  8.92 48 146 186
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 160
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
;


proc iml;

use physfit;
read all;

Y = oxygen;

X0 = j(nrow(Y),1);   /* the ONES vector */
X1 = time;
X2 = hr_rest;
X3 = hr_run;
X4 = hr_max;

R0 = Y`*X0*inv(X0`*X0)*X0`*Y;   /* SS(mu) */

X01 = X0||X1;

R01 = Y`*X01*inv(X01`*X01)*X01`*Y;

X013 = X0||X1||X3;

R013 = Y`*X013*inv(X013`*X013)*X013`*Y;

X0132 = X0||X1||X3||X2;

R0132 = Y`*X0132*inv(X0132`*X0132)*X0132`*Y;

X01324 = X0||X1||X3||X2||X4;

R01324 = Y`*X01324*inv(X01324`*X01324)*X01324`*Y;

X0124 = X0||X1||X2||X4;

R0124 = Y`*X0124*inv(X0124`*X0124)*X0124`*Y;

X0134 = X0||X1||X3||X4;

R0134 = Y`*X0134*inv(X0134`*X0134)*X0134`*Y;

X0324 = X0||X3||X2||X4;

R0324 = Y`*X0324*inv(X0324`*X0324)*X0324`*Y;

R1324_0 = R01324-R0;   /* SS(Regr) */


   /* Sequential Sums of Squares in order 1,3,2,4 */

R1_0 = R01-R0;

R3_01 = R013-R01;

R2_013 = R0132-R013;

R4_0132 = R01324-R0132;

SSRESID = Y`*Y - Y`*X01324*inv(X01324`*X01324)*X01324`*Y;

print R1324_0 R1_0 R3_01 R2_013 R4_0132 SSRESID;


   /* Partial Sums of Squares */

R1_0324 = R01324-R0324;

R3_0124 = R01324-R0124;

R2_0134 = R01324-R0134;

print R1_0324 R3_0124 R2_0134 R4_0132;

quit;

   /* Note:  Practical computer algorithms for regression use much 
      more efficient methods for computing these sums of squares. */


proc reg;
model oxygen = time hr_run hr_rest hr_max / SS1 SS2;

/* Note the order in which the independent variables are listed! */
/* SS1 gives sequential sums of squares, SS2 gives partial sums of squares */

run;