options ls=79 ps=50 nocenter;

/* read data (level is score as a factor) */
data dat;
   infile 'forage.dat' firstobs=2;
   input cow $ hc $ trt per $ dmi code;
   ldmi = log10(dmi);

proc glm;				/* proc glm works OK */
   class cow hc trt per;
   model ldmi = hc|trt cow(hc*trt) per per*hc per*trt per*hc*trt;
   test h = hc|trt e = cow(hc*trt);
   contrast 'linear trt' trt -2 -1 0 1 2 / e = cow(hc*trt);

/* split plot analysis */
proc mixed;				/* save LS means for whole plot */
   class cow hc trt per;
   model ldmi = hc|trt|per;
   random cow(hc*trt);
   contrast 'linear trt' trt -2 -1 0 1 2;
   lsmeans hc*trt;
   make 'lsmeans' out=lsm;		/* save LS means dataset */

data plot; set lsm;			/* pull out substrings with levels */
   hc = substr(level,8,1);
   trt = substr(level,10,1);
proc plot data=plot;			/* interaction plot */
   plot lsmean*trt=hc;

proc mixed data=dat;			/* save LS means for split plot */
   class cow hc trt per;
   model ldmi = hc|trt|per;
   random cow(hc*trt);
   lsmeans hc*trt*per;
   make 'lsmeans' out=lsmsub;

data plotsub; set lsmsub;
   hc = substr(level,12,1);
   trt = substr(level,14,1);
   per = substr(level,16,6);
proc sort data=plotsub; by per;
proc plot; by per;
   plot lsmean*trt=hc;
run;