options ls=79 ps=100;

   /* B=BARN M=MIXED A=ALFALFA */
/*
data a;
   infile 'depies.dat' firstobs=2;
   input cow trt $ pas barn pad parity $ wk kmy kfat pfat kpro ppro scc fcm;
proc sort; by parity trt pad wk cow;

proc glm;
   class barn pas pad parity wk cow;
   model kmy = barn pas(barn) pad(barn*pas)
      parity parity*barn parity*pas(barn) cow(parity*pad*barn*pas)
      wk wk*barn wk*pas(barn) wk*parity wk*barn*parity wk*parity*pas(barn);
   random pad(barn*pas) cow(parity*pad*barn*pas);
   test h = barn pas(barn) e = pad(barn*pas);
   test h = parity parity*barn parity*pas(barn) e = cow(parity*pad*barn*pas);
   lsmeans barn pas(barn) / e = pad(barn*pas) pdiff stderr;

proc mixed;
   class barn pas pad parity wk cow;
   model kmy = barn|pas(barn)|parity|wk;
   random pad(barn*pas) cow(parity*pad*barn*pas);
   lsmeans pas(barn) barn;

proc mixed;
   class barn pas pad parity wk cow;
   model kmy = barn|pas(barn)|parity|wk;
   random pad(barn*pas) cow(parity*pad*barn*pas);
   lsmeans pas(barn) barn;

proc mixed;
   class barn pas pad parity wk cow;
   model kmy = barn|pas(barn)|parity|wk;
   random pad(barn*pas);
   repeated / type=cs subject = cow(parity*pad);
   lsmeans barn pas(barn);
*/

data b;
   infile 'depies3.dat' firstobs=2;
   input cow trt $ pad parity $
      kmy1 kmy2 kmy3 kmy4 kmy5 kmy6 kmy7 kmy8 kmy9 kmy10 kmy11 kmy12 kmy13;
   /* trt: M=MIXED A=ALFALFA */

proc glm;
   class trt pad parity cow;
   model kmy1-kmy13 = trt pad(trt) parity parity*trt;
   test h = trt e = pad(trt);
   lsmeans parity*trt / out=lsm;

data c; set lsm;
   week = substr(_NAME_,4,2);
proc sort; by trt;
proc plot; by trt;
   plot lsmean*week=parity;

proc glm data=b;
   class trt pad parity cow;
   model kmy1-kmy13 = trt pad(trt) parity parity*trt / nouni;
   repeated week 13 (1 2 3 4 5 6 7 8 9 10 11 12 13) polynomial / summary;