options ls=80 ps=64;
data a;
   infile 'tomato.dat' firstobs=2;
   input entry yr mfwlog tg430 $;

/* one-factor anova for each year */
proc sort; by yr;
proc glm; by yr;
   class tg430;
   model mfwlog = tg430;
   /* note for contrasts that alphabetic order is A,B,H */
   contrast 'additive' tg430 1 -1 0;
   estimate 'additive' tg430 1 -1 0;
   contrast 'dominance' tg430 -1 -1 2;
   estimate 'dominance' tg430 -1 -1 2 / divisor = 2;
   means tg430;
   lsmeans tg430 / stderr pdiff;

/* histogram type summaries */
/* caution: this probably produces more output than you need */
/*          be selective! */
proc sort; by yr tg430;
proc univariate plot; by yr tg430;
   var mfwlog;

/* full split plot model */
proc mixed;
   class tg430 yr entry;
   model mfwlog = tg430 yr tg430*yr;
   random entry(tg430);

/* reduced split plot model over two years */
proc mixed;
    class tg430 yr entry;
    model mfwlog = tg430 yr / outp = blup outpm = blue;
    random entry(tg430) / solution;
    contrast 'additive' tg430 1 -1;
    estimate 'additive' tg430 1 -1;
    contrast 'dominance' tg430 -1 -1 2;
    estimate 'dominance' tg430 -1 -1 2 / divisor = 2;
    lsmeans tg430 / pdiff;
    ods output SolutionR=random;
    /* Output data set random contains the random effects */
proc print data = random;

/* sort out whole plot and split plot residuals */
data blup; set blup; /* best linear unbiased predictors */
    blup = pred;
    spresid = resid; /* split plot residual */
    drop Pred L95 U95 Resid;
data blue; set blue; /* best linear unbiased estimators */
    blue = pred; /* estimate of fixed effects */
    wpresid = resid; 
    drop Pred L95 U95 Resid;
proc sort data = blup; by tg430 yr entry;
proc sort data = blue; by tg430 yr entry;
data both; merge blue blup; by tg430 yr entry;
    wpresid = wpresid - spresid; /* whole plot residual */
proc print;
    /* notice that wpresid agree with Estimate in random dataset */
proc sort data = both; by yr;
proc plot data = both; by yr;
    plot spresid*wpresid = tg430;