* Myers and Montgomery, Example 10.5, p. 511;

* First, model mean and ln of variance separately, as suggested;
* on p. 519 of text;

data mm_eg105a;
 input x1 x2 y1-y4;
 ybar=mean(of y1-y4);
 s2=var(of y1-y4);
 lns2=log(s2);
 x1x2=x1*x2;
 x1x1=x1*x1;
 x2x2=x2*x2;
cards;
-1 -1 33.5021 41.2268 25.2683 31.9930
-1  0 35.8234 38.0689 32.7928 34.0383
-1  1 33.0773 31.8435 36.2500 34.0162
 0 -1 30.4481 41.2870 15.1493 23.9883
 0  0 34.8679 40.2276 27.7724 31.1321
 0  1 35.2202 37.1008 33.3280 35.2085
 1 -1 21.1553 34.1086  0.7917 15.7450
 1  0 27.6736 38.1477 15.5132 25.9873
 1  1 32.1245 38.1193 26.1673 32.1622
;
run;

proc rsreg data=mm_eg105a;
 model ybar=x1 x2;
run;

proc rsreg data=mm_eg105a;
 model lns2=x1 x2;
 ridge min radius=0 to 2 by .1;
run;

* Next, model the mean and variance as on p. 512;

data mm_eg105;
 set mm_eg105a;
 y=y1; z1=-1; z2=-1; x1z1=x1*z1; x1z2=x1*z2; x2z1=x2*z1; x2z2=x2*z2; output;
 y=y2; z1=-1; z2= 1; x1z1=x1*z1; x1z2=x1*z2; x2z1=x2*z1; x2z2=x2*z2; output;
 y=y3; z1= 1; z2=-1; x1z1=x1*z1; x1z2=x1*z2; x2z1=x2*z1; x2z2=x2*z2; output;
 y=y4; z1= 1; z2= 1; x1z1=x1*z1; x1z2=x1*z2; x2z1=x2*z1; x2z2=x2*z2; output;
 keep y x1 x2 x1x2 x1x1 x2x2 z1 z2 x1z1 x1z2 x2z1 x2z2;
run;

proc reg data=mm_eg105;
 model y = x1 x2 x1x2 x1x1 x2x2 z1 z2 x1z1 x1z2 x2z1 x2z2;
run;