options nocenter ps=64 ls=80;

data phoffman;
   infile 'season.dat' firstobs=2;
   input age herd trt cow season milk fat prot scc lscore stage;
   /* season: 1=spring 2=summer 3=fall 4=winter */
proc sort; by age;

proc glm; by age;
   classes herd trt cow season;
   model milk = herd trt cow(herd*trt) season trt*season / ss3;
   test h=herd trt e=cow(herd*trt);
   means season trt*season;
   lsmeans trt / e=cow(herd*trt) pdiff stderr;
   lsmeans trt*season / pdiff stderr out=lsm;
   output out=resid r=remilk p=prmilk;

proc plot data=resid; by age;		/* diagnostic plots */
   plot remilk*prmilk=trt;
   plot remilk*stage=trt;

proc plot data=lsm; by age;		/* interaction plots */
   plot lsmean*season=trt;

proc mixed data=phoffman; by age;
   class herd trt cow season;
   model milk = trt|season;
   random herd cow(herd*trt);
   lsmeans trt*season;

proc mixed data=phoffman; by age;
   class herd trt cow season;
   model milk = trt|season;
   random herd;
   repeated / type=cs subject=cow(herd*trt);
   lsmeans trt*season;

proc glm data=phoffman; by age;
   classes cow;
   model milk = cow;
   output out=rmilk r=rmilk;
proc glm; by age;
   classes trt cow season;
   model rmilk = cow(trt) trt stage season trt*season / ss3;
   lsmeans season trt*season;
   means season trt*season;
   output out=resids r=rrmilk p=pmilk;

proc plot; by age;
   plot rrmilk*pmilk=trt;
   plot rrmilk*stage=trt;