# Maestre_biomass_new.R
# Two-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("factor","factor","numeric"),text="
nitrogen water biomass
040  125   4.372
040  125   4.482
040  125   4.221
040  125   3.977
040  250   7.400
040  250   8.027
040  250   7.883
040  250   7.769
040  375   7.226
040  375   8.126
040  375   6.840
040  375   7.901
080  125   5.140
080  125   3.913
080  125   4.669
080  125   4.306
080  250   9.099
080  250   9.711
080  250   9.123
080  250   9.709
080  375  10.701
080  375  11.552
080  375  11.356
080  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
")

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

# Print data set
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)

# MODEL WITH INTERACTION - USE THIS OUTPUT IF INTERACTION SIGNIFICANT
aovout <- aov(y~nitrogen+water+nitrogen:water,data=biomdata)  

# ANOVA using Type I SS
anova(aovout)

# ANOVA using Type II SS
Anova(aovout)

# Calculate least squares means
emmeans(aovout,~nitrogen:water)

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

# MODEL WITHOUT INTERACTION - USE THIS OUTPUT IF INTERACTION NS
aovout_noint <- aov(y~nitrogen+water,data=biomdata)  # Model without interaction

# Type I SS
anova(aovout_noint)

# Type II SS
Anova(aovout_noint)

# Calculate least squares means
emmeans(aovout_noint,~nitrogen)
emmeans(aovout_noint,~water)

# All pairwise comparisons - Tukey method
compout <- glht(aovout_noint,linfct=mcp(nitrogen="Tukey"))
summary(compout)
confint(compout)
cld(compout)

# All pairwise comparisons - Tukey method
compout <- glht(aovout_noint,linfct=mcp(water="Tukey"))
summary(compout)
confint(compout)
cld(compout)

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