/* hald.sas */
options ls=72 ps=60;
data hald;
   input y x1 x2 x3 x4;
cards;
78.5 7 26 6 60
74.3 1 29 15 52
104.3 11 56 8 20
87.6 11 31 8 47
95.9 7 52 6 33
109.2 11 55 9 22
102.7 3 71 17 6
72.5 1 31 22 44
93.1 2 54 18 22
115.9 21 47 4 26
83.8 1 40 23 34
113.3 11 66 9 12
109.4 10 68 8 12
;
proc print;
title 'hald cement data';
proc plot data=hald;
   plot y*x1 y*x2 y*x3 y*x4 x1*x2 x1*x3 x1*x4 x2*x3 x2*x4 x3*x4;
proc corr data=hald;
   var x1 x2 x3 x4;
   with x1 x2 x3 x4;   
proc reg data=hald;
   model y = x1-x4;
      output out=fullfit p=yhat r=resid;
title 'full model with all variables';      
proc plot data=fullfit;
   plot resid*yhat / vref=0;
proc reg data=hald;
   model y = x1-x4 / selection=rsquare adjrsq cp aic sbc details;
title 'all possible regressions';
proc reg data=hald;
   model y = /noprint;
      output out=new0 press=pressresid;
   model y = x1 /noprint;
      output out=new1 press=pressresid;
   model y = x2 /noprint;
      output out=new2 press=pressresid;
   model y = x3 /noprint;
      output out=new3 press=pressresid;
   model y = x4 /noprint;
      output out=new4 press=pressresid;
   model y = x1 x2 /noprint;
      output out=new12 press=pressresid;
   model y = x1 x3 /noprint;
      output out=new13 press=pressresid;
   model y = x1 x4 /noprint;
      output out=new14 press=pressresid;
   model y = x2 x3 /noprint;
      output out=new23 press=pressresid;
   model y = x2 x4 /noprint;
      output out=new24 press=pressresid;
   model y = x3 x4 /noprint;
      output out=new34 press=pressresid;
   model y = x1 x2 x3 /noprint;
      output out=new123 press=pressresid;
   model y = x1 x2 x4 /noprint;
      output out=new124 press=pressresid;
   model y = x1 x3 x4 /noprint;
      output out=new134 press=pressresid;
   model y = x2 x3 x4 /noprint;
      output out=new234 press=pressresid;
   model y = x1 x2 x3 x4 /noprint;
      output out=new1234 press=pressresid;
proc means data=new0 noprint;
   var pressresid;
   output out=new0a uss=press;
proc means data=new1 noprint;
   var pressresid;
   output out=new1a uss=press;
proc means data=new2 noprint;
   var pressresid;
   output out=new2a uss=press;
proc means data=new3 noprint;
   var pressresid;
   output out=new3a uss=press;
proc means data=new4 noprint;
   var pressresid;
   output out=new4a uss=press;
proc means data=new12 noprint;
   var pressresid;
   output out=new12a uss=press;
proc means data=new13 noprint;
   var pressresid;
   output out=new13a uss=press;
proc means data=new14 noprint;
   var pressresid;
   output out=new14a uss=press;
proc means data=new23 noprint;
   var pressresid;
   output out=new23a uss=press;
proc means data=new24 noprint;
   var pressresid;
   output out=new24a uss=press;
proc means data=new34 noprint;
   var pressresid;
   output out=new34a uss=press;
proc means data=new123 noprint;
   var pressresid;
   output out=new123a uss=press;
proc means data=new124 noprint;
   var pressresid;
   output out=new124a uss=press;
proc means data=new134 noprint;
   var pressresid;
   output out=new134a uss=press;
proc means data=new234 noprint;
   var pressresid;
   output out=new234a uss=press;
proc means data=new1234 noprint;
   var pressresid;
   output out=new1234a uss=press;
data press1;
   set new0a
       new1a   new2a   new3a   new4a
       new12a  new13a  new14a  new23a  new24a new34a
       new123a new124a new134a new234a
       new1234a;
data labels;
   input model @@;
cards;
0 1 2 3 4 12 13 14 23 24 34 123 124 134 234 1234
;
data press2;
   merge labels press1;
proc sort data=press2 out=press3;
   by press;
proc print data=press3;
   var model press;
title 'press selection';
data cpplot;
   input number R2 R2adj Cp model $;
   p = number + 1;
cards;
         2       0.9787      0.9744      2.6782    e
         2       0.9725      0.9670      5.4959    f
         3       0.9823      0.9764      3.0182    k
         3       0.9823      0.9764      3.0413    l
         3       0.9813      0.9750      3.4968    m
         3       0.9728      0.9638      7.3375    n
         4       0.9824      0.9736      5.0000    o
         1        .           .          2         +
         2        .           .          2.5       +
         3        .           .          3         +
         4        .           .          3.5       +
         5        .           .          4         +
         6        .           .          4.5       +
         7        .           .          5         +
;
proc plot data=cpplot;
   plot Cp*p=model;
title 'Cp plot for the seven best models';
proc reg data=hald;
   model y = x4 x1 / ss1 ss2 seqb;
      output out=final1 p=yhat14 r=resid14;
title 'analysis for the model with x0 x1 and x4';
proc print data=final1;
   var y x1 x4 yhat14 resid14;
proc plot data=final1;
   plot resid14*yhat14 / vref=0;
proc univariate data=final1 plot normal;
     var resid14;
proc rank data=final1 normal=blom out=newf1;
     var resid14;
     ranks normsc;
proc print data=newf1;
   var resid14 normsc;
proc plot data=newf1;
     plot resid14*normsc;
proc reg data=hald;
   model y = x2 x1 / ss1 ss2 seqb;
      output out=final2 p=yhat12 r=resid12;
title 'analysis for the model with x0 x1 and x2';
proc print data=final2;
   var y x1 x2 yhat12 resid12;
proc plot data=final2;
   plot resid12*yhat12 / vref=0;
proc univariate data=final2 plot normal;
     var resid12;
proc rank data=final2 normal=blom out=newf2;
     var resid12;
     ranks normsc;
proc print data=newf2;
   var resid12 normsc;
proc plot data=newf2;
     plot resid12*normsc;
run;