/* soybean.sas */
/* soybean covariate example from Allen and Cady */
options ls=72 ps=60;
data soybean1;
   infile 'soybean.dat';
   input treat $ yield height @@;
proc sort data=soybean1 out=soybean;
   by treat;
proc means data=soybean;
   var yield height;
   by treat;
title 'soybean example';
proc plot data=soybean;
   plot yield*height=treat;
/* get the adjusted height variable for the adjusted means parameterization */
proc glm data=soybean noprint;
   model height=;
      output out=new1 r=adjhgt1 p=hmean;
/* model with adjusted means parameterization */
proc glm data=new1;
   class treat;
   model yield = treat adjhgt1/noint ss1 ss2 solution clparm;
      estimate 'control mean: unadj' treat 1 0 0 adjhgt1 -1.86666667;
      estimate 'light mean: unadj'   treat 0 1 0 adjhgt1 -2.06666667;
      estimate 'shade mean: unadj'   treat 0 0 1 adjhgt1  3.93333333;
      estimate 'control mean: adj' treat 1 0 0 adjhgt1 0;
      estimate 'light mean: adj'   treat 0 1 0 adjhgt1 0;
      estimate 'shade mean: adj'   treat 0 0 1 adjhgt1 0;
      estimate 'light vs shade'    treat 0 1 -1;
      estimate 'control vs others' treat 2 -1 -1  / divisor=2;
      estimate 'light vs control'  treat -1 1  0;
      estimate 'control vs shade'  treat  1 0 -1;
         contrast 'light vs shade'    treat 0 1 -1;
         contrast 'control vs others' treat 2 -1 -1;
         contrast 'among (adjusted)'  treat  0  1 -1,
                                      treat  2 -1 -1;
      means treat / deponly clm;
      lsmeans treat / stderr;
title 'adjusted means parameterization';
/* get the adjusted height variable for the unadjusted means parameterization */
proc glm data=soybean noprint;
   class treat;
   model height=treat;
      output out=new2 r=adjhgt2 p=hmeans;
/* model with unadjusted means parameterization */
proc glm data=new2;
   class treat;
   model yield = treat adjhgt2/noint ss1 ss2 solution clparm;
      estimate 'control mean: unadj' treat 1 0 0 adjhgt2  0;
      estimate 'light mean: unadj'   treat 0 1 0 adjhgt2  0;
      estimate 'shade mean: unadj'   treat 0 0 1 adjhgt2  0;
      estimate 'control mean: adj' treat 1 0 0 adjhgt2  1.86666667;
      estimate 'light mean: adj'   treat 0 1 0 adjhgt2  2.06666667;
      estimate 'shade mean: adj'   treat 0 0 1 adjhgt2 -3.93333333;
      estimate 'light vs shade'     treat 0 1 -1;
      estimate 'control vs others'  treat 2 -1 -1 / divisor=2;
         contrast 'light vs shade'     treat 0  1 -1;
         contrast 'control vs others'  treat 2 -1 -1;
         contrast 'among (unadjusted)' treat 0  1 -1,
                                       treat 2 -1 -1;
      means treat / deponly clm;
      lsmeans treat / stderr;
title 'unadjusted means parameterization';
/* alternate version of the model with original height variable */
proc glm data=soybean;
   class treat;
   model yield = treat height / ss1 ss2 clparm;
      estimate 'light vs shade'     treat 0  1 -1;
      estimate 'control vs others'  treat 2 -1 -1 / divisor=2;
         contrast 'light vs shade'  treat 0  1 -1;
         contrast 'shade vs others' treat 2 -1 -1;
      means treat / deponly clm;
      lsmeans treat / stderr;
title 'alternate approach without adjusting height variable';
/* get the Scheffe intervals for pairwise differences for adjusted means*/
data interval;
   input differ $ estimate stderr;
   multi = sqrt(2*finv(.95,2,41));
   lowerCL = estimate - multi*stderr;
   upperCL = estimate + multi*stderr;
cards;
L-S  6.7435425     0.04900394
L-C  3.6116736     0.04714147
C-S  3.1318688     0.04888388
;
proc print data=interval;
   var differ estimate stderr lowerCL upperCL;
title 'Scheffe intervals for pairwise differences for adjusted means';
run;