options ls=78 ps=50;

data frost;
   infile 'frost.dat' firstobs=2;
   input catrt accn hardy rep y frost sps;
   if hardy = 1;		/* only consider hardy plants */

/* average across subsamples *
/* since hardy*catrt*accn(sps) is the experimental unit */
proc sort; by hardy sps accn catrt;
proc means noprint; by hardy sps accn catrt;
   var y frost;
   output out=subsamp mean=y frost;

/* overall analysis to justify separating by catrt */
proc glm;
   class sps catrt;
   model frost y = sps catrt sps * catrt;
proc glm;
   class sps catrt;
   model frost = sps catrt;
   lsmean sps catrt / pdiff stderr;
   output out = diag p = pred r = resid;

data jitter; set diag;		/* 4% jitter of pred for plot */
   jitter = pred + .04 * 3 * (ranuni(0) - .5);
proc plot;
   plot resid * jitter = catrt;

/* now investigate variation among accessions */

proc sort data = frost; by catrt hardy accn sps rep;
proc means noprint; by catrt hardy accn sps rep;
   var y frost;
   output out = rep mean = y frost;

proc glm;
   class rep accn sps catrt;
   model y = rep sps | catrt accn(sps) catrt*accn(sps);
   random accn(sps) catrt*accn(sps) / test;

proc sort; by catrt;
proc mixed; by catrt;
   class rep accn sps;
   model y = rep sps / solution;
   random accn(sps) / solution;

/* additional investigation: adjust y for frost in ancova */
/* but this can get very messy! */