options nocenter ps=64 ls=80;

data season;
   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 */
   lscc = log10(scc);
proc sort; by age;

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

proc plot data=resid; by age;	
   plot relscc*prlscc=trt;
   plot relscc*stage=trt;

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

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

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

proc glm data=season; by age;
   classes cow;
   model lscc = cow;
   output out=rlscc r=rlscc;
proc glm; by age;
   classes trt cow season;
   model rlscc = cow(trt) trt stage season trt*season / ss3;
   lsmeans season trt*season;
   means season trt*season;
   output out=resids r=rrlscc p=plscc;

proc plot; by age;
   plot rrlscc*plscc=trt;
   plot rrlscc*stage=trt;