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;

/* repeated measures analysis */
proc glm;	/* historical method; approximate for unbalanced data */
   class tree soil reps;
   model score1 score2 score3 score4 score5 score6 = tree soil(tree);
   repeated time 6 (1 2 3 4 5 6) polynomial / short printe summary;

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;	/* split plot with compound symmetry in mixed model */
   class tree soil reps time;
   model score = tree soil(tree) time time*tree time*soil(tree) / solution;
   repeated / type=cs sub = reps(tree*soil);
   lsmeans time*soil(tree);

proc glm data = tree;	/* split plot approximation */
   class tree soil reps time;
   model score = tree soil(tree) reps(tree*soil)
      time time*tree time*soil(tree) / ss3;
   test h = tree soil(tree) e = reps(tree*soil);
*  lsmeans tree soil(tree) / stderr pdiff e=reps(tree*soil);
*  lsmeans time*tree / stderr;
   lsmeans time*soil(tree) / stderr out=lsm;

/* interaction plots by tree, soil and time */
data control; set lsm;
   if tree = 0;
   lsmctrl = lsmean;
proc sort; by tree time;
data trees; set lsm;
   if tree = 1;
   lsmtree = lsmean;
proc sort; by tree time;
data; merge trees control; by tree time;
proc plot hprecent = 50 vpercent = 50;
   plot (lsmtree lsmctrl) * time = soil;

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

proc sort data = soil; by tree;
proc anova; by tree;	/* historical setup for balanced repeated meas */
   class soil reps;
   model score1 score2 score3 score4 score5 score6 = soil / nouni;
   repeated time 6 (1 2 3 4 5 6) polynomial / summary;

proc sort data = tree; by tree;
proc mixed; by tree;	/* split plot with compound symmetry in mixed model */
   class soil reps time;
   model score = soil time time*soil;
   repeated / type = cs sub = reps(soil);
   lsmeans time*soil;

proc anova; by tree;	/* split plot approximation */
   class soil reps time;
   model score = soil reps(soil) time time*soil;
   test h = soil e = reps(soil);
run;