title 'Brassica family tolerance field study analysis';
options ls = 79 ps = 56;
data a;
   infile 'brassica.dat' missover firstobs = 2;
   input crop pyrrat block inj7dat in14dat yldkga pconyld pcnblk year trt $;
   logyld = log10(yldkga);
proc sort; by year;

proc glm; by year;
   classes block crop pyrrat;
   model yldkga logyld = block crop crop * block pyrrat crop * pyrrat;
   test h = crop e = crop * block;
   random block crop * block;

proc mixed; by year;
   classes block crop pyrrat;
   model yldkga = crop | pyrrat;
   random block crop * block;
proc mixed; by year;
   classes block crop pyrrat;
   model logyld = crop | pyrrat;
   random block crop * block;

/* combined full analysis over years */

proc glm;
   classes year block crop pyrrat;
   model yldkga logyld = year block(year) crop crop * year crop * block(year) 
      pyrrat pyrrat * year pyrrat * crop pyrrat * year * crop;
   test h = year e = block(year);
   test h = crop crop * year e = crop * block(year);
   random block(year) crop * block(year);

proc mixed;
   classes year block crop pyrrat;
   model yldkga = crop | year | pyrrat;
   random block(year) crop * block(year);
proc mixed;
   classes year block crop pyrrat;
   model logyld = crop | year | pyrrat;
   random block(year) crop * block(year);

proc glm;
   classes year block crop pyrrat;
   model yldkga logyld = year block(year) crop crop * year crop * block(year) 
      pyrrat crop * pyrrat;
   test h = year e = block(year);
   test h = crop crop * year e = crop * block(year);
   random block(year) crop * block(year);
   lsmeans crop * year / e = crop * block(year) stderr out = lsmcy;
   lsmeans crop * pyrrat / stderr out = lsmcp;

proc print data = lsmcp;
proc plot data = lsmcp;
   plot lsmean * pyrrat = crop;
proc plot data = lsmcy;
   plot lsmean * year = crop;

proc mixed data = a;
   classes year block crop pyrrat;
   model yldkga = crop | pyrrat year crop * year;
   random block(year) crop * block(year);
   lsmeans crop * year;
   lsmeans crop * pyrrat;
proc mixed data = a;
   classes year block crop pyrrat;
   model logyld = crop | pyrrat year crop * year;
   random block(year) crop * block(year);
   lsmeans crop * year;
   lsmeans crop * pyrrat;