* Maestre_biomass.sas;
options pageno=1 linesize=80;
goptions reset=all;
title "Two-way ANOVA for biomass";
title2 "Data from Maestre and Reynolds (2007)";
data maestre;
	input nitrogen water biomass;
	* Apply transformations here;
	y = log10(biomass);
	datalines;
 40  125   4.372
 40  125   4.482
 40  125   4.221
 40  125   3.977
 40  250   7.400
 40  250   8.027
 40  250   7.883
 40  250   7.769
 40  375   7.226
 40  375   8.126
 40  375   6.840
 40  375   7.901
 80  125   5.140
 80  125   3.913
 80  125   4.669
 80  125   4.306
 80  250   9.099
 80  250   9.711
 80  250   9.123
 80  250   9.709
 80  375  10.701
 80  375  11.552
 80  375  11.356
 80  375   9.759
120  125   5.021
120  125   4.970
120  125   5.055
120  125   4.862
120  250   9.029
120  250  10.791
120  250   9.115
120  250  10.319
120  375  12.189
120  375  14.381
120  375  13.153
120  375  14.066
;
run;
* Print data set;
proc print data=maestre;
run;
* Plot means, standard errors, and observations;
proc gplot data=maestre;
	plot y*nitrogen=water / vaxis=axis1 haxis=axis1 legend=legend1;
	symbol1 i=std1mjt v=star height=2 width=3;
	axis1 label=(height=2) value=(height=2) width=3 major=(width=2) minor=none;
	legend1 label=(height=2) value=(height=2);
run;
* Two-way ANOVA with all fixed effects;
proc glm data=maestre;
	class nitrogen water;
	model y = nitrogen water nitrogen*water;
    lsmeans nitrogen water / adjust=tukey cl lines;
	output out=resids p=pred r=resid;
run;
goptions reset=all;
title "Diagnostic plots to check anova assumptions";
* Plot residuals vs. predicted values;
proc gplot data=resids;
	plot resid*pred=1 / vaxis=axis1 haxis=axis1;
	symbol1 v=star height=2 width=3;
	axis1 label=(height=2) value=(height=2) width=3 major=(width=2) minor=none;run;
* Normal quantile plot of residuals;
proc univariate noprint data=resids;
	qqplot resid / normal waxis=3 height=4;
run;
quit;
