* Myers and Montgomery, Example 11.2, Table 11.4, p. 571;

* MM example 11.2, p. 554;
data mm_eg112;
 input x1 x2 x3 y;
 if x1>0 then x1=1/x1;
 if x2>0 then x2=1/x2;
 if x3>0 then x3=1/x3;
 x12=x1*x2;
 x13=x1*x3;
 x23=x2*x3;
 x1212=x12*(x1-x2);
 x1313=x13*(x1-x3);
 x2323=x23*(x2-x3);
 x123=x12*x3;
 group=100*x1+10*x2+x3;
cards;
1   0   0   540
1   0   0   560
0   1   0   330
0   1   0   350
0   0   1   295
0   0   1   260
2   2   0   610
0   2   2   330
2   0   2   425
1.5 6   6   710
6   1.5 6   640
6   6   1.5 460
3   3   3   800
3   3   3   850
;
run;

/*
Here, we successively fit linear, quadratic, special cubic and full
cubic models.  Root MSE values from fits 1-4 show up in bottom part of
Table 11.5, p. 573. Note that the R^2 values do not match up (the
noint option messes these up). The top part of the table is obtained
by the tests denoted linear_term, quadratic_term, spec_cubic_term, and
full_cubic_term; the residual is from fit 4. In the middle part of the
table, SSPE and MSPE are obtained from the proc glm code. The lof SS
and MS come from each individual fit and the PE glm fit. This is
illustrated by the results in Table 11.6, p. 574, which come from fit
3.
*/

proc reg data=mm_eg112;
* Fit 1;
 model y=x1 x2 x3/noint;
 linear_term: test x1=x2=x3;
* Fit 2;
 model y=x1 x2 x3 x12 x13 x23/noint;
 quad_term: test x12=0, x13=0, x23=0;
 quad_model: test x1=x2=x3, x12=0, x13=0, x23=0;
* Fit 3;
 model y=x1 x2 x3 x12 x13 x23 x123/noint;
 spec_cubic_term: test x123=0;
 special_cubic_model: test x1=x2=x3, x12=0, x13=0, x23=0, x123=0;
* Fit 4;
 model y=x1 x2 x3 x12 x13 x23 x123 x1212 x1313/noint;
 full_cubic_term: test x1212=0, x1313=0;
run;
quit;

proc glm data=mm_eg112;
 class group;
 model y=group;
run;

/*
The results show that the special cubic is the model of choice.
We fit it and generate predicted values on the simplex.
*/

data pred;
 do x1=0 to 1 by .1;
  do x2=0 to 1-x1 by .1;
   x3=1-x1-x2;
   x12=x1*x2;
   x13=x1*x3;
   x23=x2*x3;
   x123=x12*x3;
   output;
  end;
 end;
run;

data mm_eg112a;
 set mm_eg112 pred;
run;

proc reg data=mm_eg112a;
 model y=x1 x2 x3 x12 x13 x23 x123/noint;
 output out=mm_eg112_out p=pred r=resid rstudent=studres;
run;
quit;