/* Fit random effects models to
/* the milk protein data */

options ls=65;

data set2;
  infile 'c:\mydocuments\courses\st565\data\dhlz.example1_4.data';
  input diet cow week protein;
  if(protein ge 9.99) then protein='.';
  run;


proc sort data=set2; by diet week; run;


/* Fit cubic trends across time with a compound
   symmetry covariance structure  */

proc mixed data=set2;
   class diet cow;
   model protein = diet diet*week diet*week*week 
                   diet*week*week*week  
                 /  s outpm=means dfm=kr;
   random intercept / subject=cow g s;
run;



/* plot the fitted curves */


axis1 label=(f=swiss h=1.8 a=90 r=0 "Protein (percent)")
      order = 2.5 to 4.5 by .5
      value=(f=swiss h=1.8)  w=3.0
      length= 4.0in;

axis2 label=(f=swiss h=2.0 "Time(weeks)")
      order = 0 to 20 by 5
      value=(f=swiss h=1.8) w= 3.0
      length = 6.5 in;

 SYMBOL1 V=CIRCLE H=1.7 w=3 l=1 i=join ;
 SYMBOL2 V=DIAMOND H=1.7 w=3 l=3 i=join ;
 SYMBOL3 V=square H=1.7 w=3 l=9 i=join ;

 PROC GPLOT DATA=means; 
   PLOT pred*week=diet /
             vaxis=axis1 haxis=axis2;
   TITLE1 ls=0.01in H=2.0 F=swiss "Estimated Protein Content";
           footnote ls=0.01in;
   RUN;




/* Fit a model with random coefficients */

proc mixed data=set2;
   class diet cow;
   model protein = diet diet*week diet*week*week 
                   diet*week*week*week
                 /  s outpm=means outp=pred dfm=kr;
   random int week week*week / subject=cow s g type=un v vcorr;
run;


/* Plot predicted values for some cows */

data predt; set pred; if(cow<6); run;

proc sort data=predt; by cow week; run;

proc print data=predt; run;


axis1 label=(f=swiss h=1.8 a=90 r=0 "Protein (percent)")
      order = 2.5 to 4.5 by .5
      value=(f=swiss h=1.8)  w=3.0
      length= 4.0in;

axis2 label=(f=swiss h=2.0 "Time(weeks)")
      order = 0 to 20 by 5
      value=(f=swiss h=1.8) w= 3.0
      length = 6.5 in;

 SYMBOL1 V=none H=1.7 w=3 l=1 i=join ;
 SYMBOL2 V=none H=1.7 w=3 l=3 i=join ;
 SYMBOL3 V=none H=1.7 w=3 l=9 i=join ;
 SYMBOL4 V=none H=1.7 w=3 l=17 i=join ;
 SYMBOL5 V=none H=1.7 w=3 l=22 i=join ;
 SYMBOL6 V=none H=1.7 w=3 l=12 i=join ;

 PROC GPLOT DATA=predt; 
   PLOT pred*week=cow /
             vaxis=axis1 haxis=axis2;
   TITLE1 ls=0.01in H=2.0 F=swiss "Predicted Protein Content";
           footnote ls=0.01in;
   RUN;