/* cuzn.sas */
/* copper and zinc in minnow larvae example */
options ls=72 ps=60;
/* input the data, transform the response, and
   create new groupings */
data minnow;
   infile 'cuzn.dat';
   input cu zn protein;
   logprot = log(protein);
if cu = 0 then delete;
if cu = 0 then delete;
copper = 'high';
if cu = 37.5 or cu = 75 then copper = 'low';
if zn = 0 then delete;
proc print data=minnow;
title 'minnow data';
/* get the standard two-way ANOVA table */
proc glm data=minnow;
   class copper zn;
   model logprot = copper zn copper*zn / ss1 ss2;
title 'minnow example: summary ANOVA';
/* input values and find orthogonal polynomial contrast coefficients */
data ortho;
   input x0 x1 x2 x3;
cards;
1 1  1  1
1 2  4  8
1 3  9 27
1 4 16 64
;
proc reg data=ortho;
   model x1 = x0 / noint noprint;
      output out=ortho1 r=linear;
   model x2 = x0 x1 / noint noprint;
      output out=ortho2 r=quad;
   model x3 = x0 x1 x2 / noint noprint;
      output out=ortho3 r=cubic;
title 'orthogonal polynomial contrasts';
data orthoc;
   merge ortho1 ortho2 ortho3;
   linear2 = 2*linear;
   quad2 = 1*quad;
   cubic2 = 10*cubic;
proc print data=orthoc;
   var x0 linear quad cubic linear2 quad2 cubic2;
/* fit the cell means model */
proc glm data=minnow;
   class copper zn;
   model logprot = copper*zn / noint solution clparm ss1 ss2;
/* estimate the marginal means */
   estimate 'high copper mean'  copper*zn  1 1 1 1 0 0 0 0 /divisor=4;
   estimate 'low copper mean'   copper*zn  0 0 0 0 1 1 1 1 /divisor=4;
   estimate '375ppm zn mean'    copper*zn  1 0 0 0 1 0 0 0 /divisor=2;
   estimate '750ppm zn mean'    copper*zn  0 1 0 0 0 1 0 0 /divisor=2;
   estimate '1125ppm zn mean'   copper*zn  0 0 1 0 0 0 1 0 /divisor=2;
   estimate '1500ppm zn mean'   copper*zn  0 0 0 1 0 0 0 1 /divisor=2;
/* estimate the main effect and interaction contrasts */
   estimate 'high vs low'       copper*zn  1  1  1 1 -1 -1 -1 -1 /divisor=4;
   estimate 'zn linear'         copper*zn -3 -1  1 3 -3 -1  1  3 /divisor=20;
   estimate 'zn quadratic'      copper*zn  1 -1 -1 1  1 -1 -1  1 /divisor=8;
   estimate 'zn cubic'          copper*zn -3  9 -9 3 -3  9 -9  3 /divisor=36;
   estimate 'int1: linear'      copper*zn -3 -1  1 3  3  1 -1 -3 /divisor=10;
   estimate 'int2: quadratic'   copper*zn  1 -1 -1 1 -1  1  1 -1 /divisor=4;
   estimate 'int3: cubic'       copper*zn -3  9 -9 3  3 -9  9 -3 /divisor=18;
/* get the single degree of freedom sums of squares */
   contrast 'high vs low'       copper*zn  1  1  1  1 -1 -1 -1 -1;
   contrast 'zn linear'         copper*zn -3 -1  1  3 -3 -1  1  3;
   contrast 'zn quadratic'      copper*zn  1 -1 -1  1  1 -1 -1  1;
   contrast 'zn cubic'          copper*zn -3  9 -9  3 -3  9 -9  3;
   contrast 'int1: linear'      copper*zn -3 -1  1  3  3  1 -1 -3;
   contrast 'int2: quadratic'   copper*zn  1 -1 -1  1 -1  1  1 -1;
   contrast 'int3: cubic'       copper*zn -3  9 -9  3  3 -9  9 -3;
/* get the ANOVA table sums of squares */
   contrast 'copper main'       copper*zn  1  1  1  1 -1 -1 -1 -1;
   contrast 'zn main'           copper*zn -3 -1  1  3 -3 -1  1  3,
                                copper*zn  1 -1 -1  1  1 -1 -1  1,
                                copper*zn -3  9 -9  3 -3  9 -9  3;
   contrast 'interaction'       copper*zn -3 -1  1  3  3  1 -1 -3,
                                copper*zn  1 -1 -1  1 -1  1  1 -1,
                                copper*zn -3  9 -9  3  3 -9  9 -3;
   output out=new1 p=cellmean r=resid;
title 'minnow example: cell means model';
proc plot data=new1;
   plot resid*cellmean/vref=0;
   plot cellmean*copper=zn;
   plot cellmean*zn=copper;
proc univariate data=new1 plot;
   var resid;
proc rank data=new1 normal=blom out=new1a;
   var resid;
   ranks nscore;
proc plot data=new1a;
   plot nscore*resid;
run;