# Maestre_biomass_3way_new.R
# Three-way ANOVA for biomass (Maestre and Reynolds 2007)

# Load necessary libraries
library(car)
library(ggplot2)
library(emmeans)
library(multcomp)
library(phia)

# Read in data set
biomdata <- read.table(header=T,colClasses=c(rep("factor",3),"numeric"),text="
nitrohet nitrogen water biomass
N  040  125   4.372 
N  040  125   4.482 
N  040  125   4.221 
N  040  125   3.977 
N  040  250   7.400 
N  040  250   8.027 
N  040  250   7.883 
N  040  250   7.769 
N  040  375   7.226 
N  040  375   8.126 
N  040  375   6.840 
N  040  375   7.901 
N  080  125   5.140 
N  080  125   3.913 
N  080  125   4.669 
N  080  125   4.306 
N  080  250   9.099 
N  080  250   9.711 
N  080  250   9.123 
N  080  250   9.709 
N  080  375  10.701
N  080  375  11.552
N  080  375  11.356
N  080  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  040  125   5.458 
Y  040  125   5.017 
Y  040  125   5.479 
Y  040  125   5.714 
Y  040  250   8.972 
Y  040  250   9.234 
Y  040  250   8.032 
Y  040  250   8.372 
Y  040  375   9.464 
Y  040  375   9.563 
Y  040  375   9.385 
Y  040  375   8.226 
Y  080  125   6.616 
Y  080  125   6.909 
Y  080  125   6.851 
Y  080  125   6.098 
Y  080  250  10.792
Y  080  250  10.164
Y  080  250  10.947
Y  080  250   9.582 
Y  080  375  14.936
Y  080  375  13.607
Y  080  375  14.231
Y  080  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
")


# Apply transformations here
biomdata <- transform(biomdata,y=log10(biomass))

# Print data
biomdata

# Graphics using ggplot2
ggplot(biomdata,aes(nitrogen,y,group=water,color=water))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="line",size=1)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)
ggplot(biomdata,aes(nitrogen,y,group=nitrohet,color=nitrohet))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="line",size=1)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)
ggplot(biomdata,aes(water,y,group=nitrohet,color=nitrohet))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="line",size=1)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)

# Plots to show three-way interaction
ggplot(biomdata,aes(nitrogen,y,group=water,color=water))+
geom_jitter(size=2,position=position_jitter(width=0.1))+
stat_summary(fun.y="mean",geom="line",size=1)+
stat_summary(fun.data="mean_cl_normal",geom="errorbar",width=0.2,size=1)+
facet_grid(.~nitrohet)

# MODEL WITH ALL FOUR INTERACTIONS
aovout <- aov(y~nitrohet*nitrogen*water,data=biomdata)

# ANOVA using Type II SS
Anova(aovout)

# MODEL WITH ONLY SIGNIFICANT INTERACTIONS
aovout_nonsint <- aov(y~nitrohet+nitrogen+water+nitrohet:water+nitrogen:water,data=biomdata) 

# ANOVA using Type II SS
Anova(aovout_nonsint)

# Calculate least squares means
emmeans(aovout_nonsint,~nitrohet:water)

# Tests for simple effects 
testInteractions(aovout_nonsint,fixed="water",across="nitrohet",adjustment="none")
testInteractions(aovout_nonsint,fixed="nitrohet",across="water",adjustment="none")

# Diagnostic plots to check ANOVA assumptions
p <- predict(aovout)
r <- resid(aovout)
plot(p,r)
qqnorm(r)
