* Maestre_biomass_3way.sas;
options pageno=1 linesize=80;
goptions reset=all;
title "Three-way ANOVA for biomass";
title2 "Data from Maestre and Reynolds (2007)";
data maestre;
	input nitrohet $ nitrogen water biomass;
	* Apply transformations here;
	y = log10(biomass);
	datalines;
N   40  125   4.372 
N   40  125   4.482 
N   40  125   4.221 
N   40  125   3.977 
N   40  250   7.400 
N   40  250   8.027 
N   40  250   7.883 
N   40  250   7.769 
N   40  375   7.226 
N   40  375   8.126 
N   40  375   6.840 
N   40  375   7.901 
N   80  125   5.140 
N   80  125   3.913 
N   80  125   4.669 
N   80  125   4.306 
N   80  250   9.099 
N   80  250   9.711 
N   80  250   9.123 
N   80  250   9.709 
N   80  375  10.701
N   80  375  11.552
N   80  375  11.356
N   80  375   9.759 
N  120  125   5.021
N  120  125   4.970
N  120  125   5.055
N  120  125   4.862
N  120  250   9.029
N  120  250  10.791
N  120  250   9.115 
N  120  250  10.319
N  120  375  12.189
N  120  375  14.381
N  120  375  13.153
N  120  375  14.066
Y   40  125   5.458 
Y   40  125   5.017 
Y   40  125   5.479 
Y   40  125   5.714 
Y   40  250   8.972 
Y   40  250   9.234 
Y   40  250   8.032 
Y   40  250   8.372 
Y   40  375   9.464 
Y   40  375   9.563 
Y   40  375   9.385 
Y   40  375   8.226 
Y   80  125   6.616 
Y   80  125   6.909 
Y   80  125   6.851 
Y   80  125   6.098 
Y   80  250  10.792
Y   80  250  10.164
Y   80  250  10.947
Y   80  250   9.582 
Y   80  375  14.936
Y   80  375  13.607
Y   80  375  14.231
Y   80  375  12.038
Y  120  125   7.389 
Y  120  125   6.683 
Y  120  125   7.759 
Y  120  125   6.752 
Y  120  250  10.731
Y  120  250  12.640
Y  120  250  10.350
Y  120  250  11.550
Y  120  375  14.697
Y  120  375  17.826
Y  120  375  14.711
Y  120  375  13.614
;
run;
* Print data set;
proc print data=maestre;
run;
proc gplot data=maestre;
	plot y*nitrohet=nitrogen y*nitrogen=water y*nitrohet=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;
* Sort data by nitrohet levels;
proc sort data=maestre;
	by nitrohet;
run;
* Plots to show three-way interaction;
proc gplot data=maestre;
	by nitrohet;
	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;
* Three-way ANOVA with interaction;
title3 "MODEL WITH ALL FOUR INTERACTIONS";
proc glm data=maestre;
	class nitrohet nitrogen water;
	model y = nitrohet|nitrogen|water / ss2;
	output out=resids p=pred r=resid;
run;
* Three-way ANOVA dropping ns interactions;
title3 "MODEL WITH ONLY SIGNIFICANT INTERACTIONS";
proc glm data=maestre;
	class nitrohet nitrogen water;
	model y = nitrohet nitrogen water nitrohet*water nitrogen*water / ss2;
	lsmeans nitrohet*water / slice=water slice=nitrohet;
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;
