options nocenter ls=80 ps=50;
/* hc = 1 for young cows, 2 for older cows */
/* trt = 1..5 is assigned to cow(hc) */
/* diet = 1..10 changes over time according to trt */
/* alfa = 38,48,58,68,78,88,98 is associated with diet */
/* alfa can increase for given trt from time 1 to time 2 */
data a;
   infile 'chew2.dat';
   input cow hc trt time diet alfa eat rum tot eatp rump totp;
   alfalin = alfa;
/* interaction plots */
proc glm noprint;
   class hc trt time;
   model eatp = hc*trt*time;
   lsmeans hc*trt*time / out=lsm;
proc sort data=lsm; by time;
proc plot; by time;
   plot lsmean*trt=hc;
/* split plot with trt */
proc glm data=a;
   class trt time cow hc alfa;
   model eatp = trt hc hc*trt trt*cow(hc)
      time hc*time trt trt*hc time*trt*hc;
   test h = trt hc hc*trt e = trt*cow(hc);
proc glm data=a;
   class trt time cow hc alfa;
   model eatp = trt hc trt*cow(hc)
      time hc*time trt trt*hc;
   test h = trt hc e = trt*cow(hc);
/* naive split plot with alfa covariate instead of diet */
proc glm data=a;
   class trt time cow hc;
   model eatp = trt hc cow(hc) time hc*time alfa alfa*hc;
   test h = alfa hc e = cow(hc);
/* separate covariate into whole- and split- plot portions */
proc sort data=a; by cow hc;
proc means noprint; by cow hc;
   var alfa;
   output out=mean mean=malfa;
data b; merge a mean; by cow hc;
   dalfa = alfa - malfa;
proc glm data=b;
   class time cow hc;
   model eatp = malfa hc malfa*hc malfa*cow(hc) time hc*time dalfa dalfa*hc;
   test h = hc malfa malfa*hc e = malfa*cow(hc);
proc glm data=b;
   class time cow hc;
   model eatp = malfa hc malfa*cow(hc) time hc*time dalfa;
   test h = hc malfa e = malfa*cow(hc);
   /* notice that least squares means not computed properly -- annoying! */
   lsmeans hc*time / stderr pdiff;