/* scab.sas */
/* potato scab disease example */
options ls=72 ps=60;
data scab;
   infile 'scab.dat';
   input time $ sulfur $ index;
/* get the standard two-way ANOVA table */
proc glm data=scab;
   class time sulfur;
   model index = time sulfur time*sulfur / ss1 ss2;
title 'scab ANOVA standard form';   
/* fit the cell means model */
proc glm data=scab;
   class time sulfur;
   model index = time*sulfur / noint solution clparm ss1 ss2;
/* estimate the time and level means */
      estimate 'fall'   time*sulfur 1 1 1 0 0 0 / divisor=3;
      estimate 'spring' time*sulfur 0 0 0 1 1 1 / divisor=3;
      estimate 'low'    time*sulfur 0 1 0 0 1 0 / divisor=2;
      estimate 'medium' time*sulfur 0 0 1 0 0 1 / divisor=2;
      estimate 'high'   time*sulfur 1 0 0 1 0 0 / divisor=2;
/* estimate the main effect and interaction contrasts */
      estimate 'fall vs spring'     time*sulfur 1  1  1 -1 -1 -1 / divisor=3;
      estimate 'low vs medium'      time*sulfur 0  1 -1  0  1 -1 / divisor=2;
      estimate 'high vs others'     time*sulfur 2 -1 -1  2 -1 -1 / divisor=4;
      estimate 'int1: low vs med'   time*sulfur 0  1 -1  0 -1  1;
      estimate 'int2: hi vs others' time*sulfur 2 -1 -1 -2  1  1 / divisor=2;
/* get the single degree of freedom sums of squares */
      contrast 'fall vs spring'     time*sulfur 1  1  1 -1 -1 -1;
      contrast 'low vs medium'      time*sulfur 0  1 -1  0  1 -1;
      contrast 'high vs others'     time*sulfur 2 -1 -1  2 -1 -1;
      contrast 'int1: low vs med'   time*sulfur 0  1 -1  0 -1  1;
      contrast 'int2: hi vs others' time*sulfur 2 -1 -1 -2  1  1;
/* get the ANOVA table sums of squares */
      contrast 'time main'   time*sulfur 1  1  1 -1 -1 -1;
      contrast 'level main'  time*sulfur 0  1 -1  0  1 -1,
                             time*sulfur 2 -1 -1  2 -1 -1;
      contrast 'interaction' time*sulfur 0  1 -1  0 -1  1,
                             time*sulfur 2 -1 -1 -2  1  1;
title 'scab analysis: cell means model (version 1)';
/* input cell means and create interaction plots */
data means;
   input time $ sulfur $ mean;
cards;
fall   high        5.75
fall   low         9.50
fall   medium     15.50
spring high       14.25
spring low        16.75
spring medium     18.25
;
proc plot data=means;
   plot mean*time=sulfur;
   plot mean*sulfur=time;
title 'interaction plots';
/* analysis with orthogonal polynomial contrasts */
proc glm data=scab;
   class time sulfur;
   model index = time*sulfur / noint solution clparm ss1 ss2;
/* estimate the main effect and interaction contrasts */
      estimate 'fall vs spring'     time*sulfur  1  1  1 -1 -1 -1 / divisor=3;
      estimate 'level: linear'      time*sulfur  5 -4 -1  5 -4 -1 / divisor=28;
      estimate 'level: quadratic'   time*sulfur  3  6 -9  3  6 -9 / divisor=36;
      estimate 'int1: linear'       time*sulfur  5 -4 -1 -5  4  1 / divisor=14;
      estimate 'int2: quadratic'    time*sulfur  3  6 -9 -3 -6  9 / divisor=18;
/* get the single degree of freedom sums of squares */
      contrast 'fall vs spring'     time*sulfur  1  1  1 -1 -1 -1;
      contrast 'level: linear'      time*sulfur  5 -4 -1  5 -4 -1;
      contrast 'level: quadratic'   time*sulfur  3  6 -9  3  6 -9;
      contrast 'int1: linear'       time*sulfur  5 -4 -1 -5  4  1;
      contrast 'int2: quadratic'    time*sulfur  3  6 -9 -3 -6  9;
/* get the ANOVA table sums of squares */
      contrast 'time main'   time*sulfur  1  1  1 -1 -1 -1;
      contrast 'level main'  time*sulfur  5 -4 -1  5 -4 -1,
                             time*sulfur  3  6 -9  3  6 -9;
      contrast 'interaction' time*sulfur  5 -4 -1 -5  4  1,
                             time*sulfur  3  6 -9 -3 -6  9;
title 'scab analysis: cell means model (version 2)';
/* input values and find orthogonal polynomial contrast coefficients */
data ortho;
   input x0 x1 x2;
cards;
1 1 1
1 2 4
1 4 16
;
proc reg data=ortho;
   model x1 = x0 / noint seqb;
      output out=ortho1 r=linear;
   model x2 = x0 x1 / noint seqb;
      output out=ortho2 r=quad;
title 'orthogonal polynomial contrasts';
data orthoc;
   merge ortho1 ortho2;
   linear2 = 3*linear;
   quad2 = 7*quad;      
proc print data=orthoc;
   var x0 linear quad linear2 quad2;
run;