options nocenter ls=80 ps=100;

data soil;
   infile 'tree.dat' firstobs=2;
   input soil$ reps acc$ score1 score2 score3 ht1 score4 ht2 score5 ht3 score6;
   if acc = 'notauto';
   tree = 0;
   if soil = 'dtac2' or soil = 'nc11004' or soil = 'ne-332' or soil = 'nm-6'
      then tree = 1;

data tree; set soil;
   score = score1; time = 1; output;
   score = score2; time = 2; output;
   score = score3; time = 3; output;
   score = score4; time = 4; output;
   score = score5; time = 5; output;
   score = score6; time = 6; output;
   drop score1--score6;

proc mixed;
   title2 'split plot as mixed model';
   class tree soil reps time;
   model score = tree soil(tree) time time*tree time*soil(tree) / solution;
   random reps(tree*soil);
   lsmeans time*soil(tree);

proc mixed;
   title2 'split plot with compound symmetry';
   class tree soil reps time;
   model score = tree soil(tree) time time*tree time*soil(tree) / solution;
   repeated / type = cs subject = reps(tree*soil);
   lsmeans time*soil(tree);

proc mixed;
   title2 'repeated measures with multivariate analysis';
   /* 
    * Hotelling-Lawley-McKeon (HLM) has better F approximation
    * than Hotelling-Lawley-Pillai-Samson (HLPS).
    * Both perform better than F for small samples.
    * Other MV test (Roy's greatest root, Wilk's lambda, Pillai trace)
    * are apparently not available.
    */
   class tree soil reps time;
   model score = tree soil(tree) time time*tree time*soil(tree) / solution;
   repeated / hlm hlps type = un subject = reps(tree*soil);
   lsmeans time*soil(tree);

/*** repeated measures analysis by tree ***/

proc sort data = tree; by tree;
proc mixed; by tree;
   title2 'mixed model as split plot';
   class soil reps time;
   model score = soil time time*soil;
   random reps(soil);
   lsmeans time*soil;

proc mixed; by tree;
   title2 'repeated measures with compound symmetry';
   class soil reps time;
   model score = soil time time*soil;
   repeated / type = cs subject = reps(soil);
   lsmeans time*soil;

proc mixed; by tree;
   title2 'repeated measures with multivariate analysis';
   class soil reps time;
   model score = soil time time*soil;
   repeated / hlm hlps type = un subject = reps(soil);
   lsmeans time*soil;
run;