library(ggplot2)
library(ggpubr)
library(tidyverse)
library(readxl)
library(multcompView)
library(rstatix)
library(multcomp)
library(plyr)
library(moments)
library(MASS)

#set workspace
setwd('/Users/kianspeck/Desktop/R_WD')
#load data
biomass <- readxl::read_xlsx('Conyza_harvest_10.1_R.xlsx', sheet = 5)

#start with shoot biomass#############
#create data table containing means, sd, se, and N
sbiomassSE <- ddply(biomass, c("Inoculation","Watering"), summarise,
                   N = length(s_biomass_DW),
                   mean = mean(s_biomass_DW),
                   sd = sd(s_biomass_DW),
                   se = sd/ sqrt(N)
)
sbiomassSE


#create theme for future use
my_theme = theme(
  axis.title = element_text(size = 12, face = 'bold'),
  axis.text.x = element_text(size = 12),
  axis.text.y = element_text(size = 12)
)

#plot
tiff("s_biomass.tiff", units="in", width=7, height=5, res=300)

sb <- ggplot(data = sbiomassSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("Shoot Dry Weight (g)")+
  theme_bw() + 
  scale_color_manual(values = c("C"="#88d8b0",
                                "M"="#ffcc5c",
                                "S"="#ff6f69")) +
  my_theme + 
  theme(legend.position = 'none')
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(biomass$s_biomass_DW)
#p > 0.05 so it IS normal
#lm
mod1 <- lm(s_biomass_DW ~ Inoculation*Watering, biomass)
#ANOVA
sb_anova <- aov(mod1)
summary(sb_anova)

#look at root biomass#############
#create data table containing means, sd, se, and N
rbiomassSE <- ddply(biomass, c("Inoculation","Watering"), summarise,
                    N = length(r_biomass_DW),
                    mean = mean(r_biomass_DW),
                    sd = sd(r_biomass_DW),
                    se = sd/ sqrt(N)
)
rbiomassSE

#plot
tiff("r_biomass.tiff", units="in", width=7, height=5, res=300)
rb <- ggplot(data = rbiomassSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("Root Dry Weight (g)")+
  theme_bw() + scale_color_manual(values = c("C"="#88d8b0",
                              "M"="#ffcc5c",
                              "S"="#ff6f69"))+
  my_theme

dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(biomass$r_biomass_DW)
#p < 0.05 it IS NOT normal
skewness(biomass$r_biomass_DW, na.rm = TRUE)
#positively skewed
#transform
biomass$r_biomass_DW <- log10(biomass$r_biomass_DW)
shapiro.test(biomass$r_biomass_DW)
#p >0.05 it IS normal
mod2 <- lm(r_biomass_DW ~ Inoculation*Watering, biomass)
rb_anova <- aov(mod2)
summary(rb_anova)

#look at root shoot ratio###############
#create data table containing means, sd, se, and N
rsratioSE <- ddply(biomass, c("Inoculation","Watering"), summarise,
                   N = length(rs_ratio),
                   mean = mean(rs_ratio),
                   sd = sd(rs_ratio),
                   se = sd/ sqrt(N)
)
rsratioSE

#plot
tiff("rs_ratio.tiff", units="in", width=7, height=5, res=300)
rs <- ggplot(data = rsratioSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black")+
                  ylab("Root:Shoot Ratio")+
                  theme_bw()

rs + scale_color_manual(values = c("C"="#88d8b0",
                                  "M"="#ffcc5c",
                                  "S"="#ff6f69"))+
  my_theme
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(biomass$rs_ratio)
#p < 0.05 so it IS NOT normal
skewness(biomass$rs_ratio)
biomass$rs_ratio <- log10(biomass$rs_ratio)
shapiro.test(biomass$rs_ratio)
#p > 0.05 so it IS normal
#lm
mod3 <- lm(rs_ratio ~ Inoculation*Watering, biomass)
#ANOVA
rsr_anova <- aov(mod3)
summary(rsr_anova)

#look at shoot water concentration###############
#create data table containing means, sd, se, and N
swaterconcentrationSE <- ddply(biomass, c("Inoculation","Watering"), summarise,
                     N = length(s_water_concen),
                     mean = mean(s_water_concen),
                     sd = sd(s_water_concen),
                     se = sd/ sqrt(N)
)
swaterconcentrationSE
#plot
tiff("s_content.tiff", units="in", width=7, height=5, res=300)
swconcen <- ggplot(data = swaterconcentrationSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("Shoot Water Concentration (g water/g biomass)")+
  theme_bw() + scale_color_manual(values = c("C"="#88d8b0",
                                     "M"="#ffcc5c",
                                     "S"="#ff6f69"))+
  my_theme + theme(legend.position = 'none')
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(biomass$s_water_concen)
#p < 0.05 so it IS NOT normal
skewness(biomass$s_water_concen, na.rm = TRUE)
biomass$s_water_concen <- log10(biomass$s_water_concen)

hist(biomass$s_water_concen); abline(v=mean(biomass$s_water_concen))


mod4 <- lm(s_water_concen ~ Inoculation*Watering, biomass)
#ANOVA
swconcen_anova <- aov(mod4)
summary(swconcen_anova)

#look at root water concentration##################
#create data table containing means, sd, se, and N
rwaterconcentrationSE <- ddply(biomass, c("Inoculation","Watering"), summarise,
                             N = length(r_water_concen),
                             mean = mean(r_water_concen),
                             sd = sd(r_water_concen),
                             se = sd/ sqrt(N)
)
rwaterconcentrationSE
#plot
tiff("r_content.tiff", units="in", width=7, height=5, res=300)

rwconcen <- ggplot(data = rwaterconcentrationSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("Root Water Concentration (g water/g biomass)")+
  theme_bw()+ 
  scale_color_manual(values = c("C"="#88d8b0",
                                  "M"="#ffcc5c",
                                  "S"="#ff6f69"))+
  my_theme
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(biomass$r_water_concen)
#p < 0.05 so it IS NOT normal
skewness(biomass$r_water_concen)
#positively skewed
biomass$r_water_concen <- log10(biomass$r_water_concen)
shapiro.test(biomass$r_water_concen) 
#p > 0.05 so it IS normal
#ANOVA
mod5 <- lm(r_water_concen ~ Inoculation*Watering, biomass)
rwconcen_anova <- aov(mod5)
summary(rwconcen_anova)

#Leaf water potential and Relative Leaf Water Content########
##first look at LWP
#look at data 9.27 data first
potential27 <- readxl::read_xlsx('LWP_measurements_R.xlsx', sheet = 3)

#creat data table including mean, sd, se, and N
leafpot27SE <- ddply(potential27, c("Inoculation","Watering"), summarise,
                   N = length(LWP),
                   mean = mean(LWP),
                   sd = sd(LWP),
                   se = sd/ sqrt(N)
)
leafpot27SE
#plot
tiff("leaf_pot.tiff", units="in", width=7, height=5, res=300)

leafpot27 <- ggplot(data = leafpot27SE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black")+ 
  ylab("Leaf Water Potential 9.27 (MPa)")+
  theme_bw()

leafpot27 + scale_color_manual(values = c("C"="#88d8b0",
                                       "M"="#ffcc5c",
                                       "S"="#ff6f69"))+
  my_theme
dev.off()

#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(potential27$LWP)
#p < 0.05 it IS NOT normal
#transform
skewness(potential27$LWP, na.rm = TRUE)
#negative skew
#add a constant to the negative values so we can transform with log
potential27$LWP <- potential27$LWP + 0.001 - min(potential27$LWP)
#log transformation                    
potential27$LWP <- log10(max(potential27$LWP+1) - potential27$LWP)
#re-test
shapiro.test(potential27$LWP)
#p > 0.05 it IS normal
mod6 <- lm(LWP ~ Inoculation*Watering, potential27)
lwp27_aov <- aov(mod6)
summary(lwp27_aov)

#load data-relative leaf water content
rlwc <- readxl::read_xlsx('LWP_measurements_R.xlsx', sheet = 5)
#Now, examine relative leaf water content##########
rlwcSE <- ddply(rlwc, c("Inoculation","Watering"), summarise,
                   N = length(RLWC_per),
                   mean = mean(RLWC_per),
                   sd = sd(RLWC_per),
                   se = sd/ sqrt(N)
)
rlwcSE
#plot
rel.lwc <- ggplot(data = rlwcSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "grey") +
  ylab("Relative Leaf Water Content (%)")+
  theme_grey()
#add colors and fix scale
rel.lwc + scale_color_manual(values = c("C"="#88d8b0",
                                       "M"="#ffcc5c",
                                       "S"="#ff6f69"))+
  scale_y_continuous(n.breaks = 3, labels = c("0%","50%","100%"),
                     limits = c(0, 100))
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(rlwc$RLWC)
#p < 0.05 so it IS NOT normal
skewness(rlwc$RLWC, na.rm = TRUE)
#negative skew
rlwc$RLWC <- 1/(max(rlwc$RLWC+1)-rlwc$RLWC)
rlwc$RLWC <- log10(max(rlwc$RLWC+1)-rlwc$RLWC)
shapiro.test(rlwc$RLWC)
#p < 0.05 so it IS NOT normal
#can't due further analysis with ANOVA. May try with different model.
#mod7 <- lm(RLWC ~ Inoculation*Watering, rlwc)
#ANOVA
#rlwc_anova <- aov(mod7)
#summary(rlwc_anova)
#photosynthesis##############################

#load data
licor <- readxl::read_xlsx('Conyza_Pn.xlsx', sheet = 1)

#create data table containing means, sd, se, and N
photoSE <- ddply(licor, c("Inoculation", "Watering"),summarise,
                    N = length(Photo),
                    mean = mean(Photo),
                    sd = sd(Photo),
                    se = sd/ sqrt(N)
)
photoSE

#plot
tiff("photosynthesis.tiff", units="in", width=7, height=5, res=300)

pn <- ggplot(data = photoSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("Photosynthetic Rate")+
  theme_bw() + 
  scale_color_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69"))+
  my_theme

dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
mod8 <- lm(Photo ~ Inoculation*Watering, licor)
#Shapiro-Wilk test
shapiro.test(licor$Photo)
#p > 0.05 so it IS normal
photo_anova <- aov(mod8)
summary(photo_anova)

#does the order of the measurement affect the data?
modtest <- lm(Photo ~ Obs, licor)
test_anova <- aov(modtest)
summary(test_anova)

#stomatal conductance###########################
condSE <- ddply(licor, c("Inoculation", "Watering"),summarise,
              N = length(Cond),
              mean = mean(Cond),
              sd = sd(Cond),
              se = sd/ sqrt(N)
)
condSE

#plot
tiff("conductance.tiff", units="in", width=7, height=5, res=300)

sc <- ggplot(data = condSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black")  +
  ylab("Stomatal Conductance")+
  theme_bw() + scale_color_manual(values = c("C"="#88d8b0",
                                  "M"="#ffcc5c",
                                  "S"="#ff6f69"))+
  my_theme + theme(legend.position = 'none')
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(licor$Cond)
#p > 0.05 so it IS normal
mod9 <- lm(Cond ~ Inoculation*Watering, licor)
cond_anova <- aov(mod9)
summary(cond_anova)


#nutrient analysis#########################
#start with nitrogen
nutrient <- readxl::read_xlsx('NP_measurement.xlsx', sheet = 2)

#create data table containing means, sd, se, and N
nitrogenSE <- ddply(nutrient, c("Inoculation", "Watering"),summarise,
                 N = length(nitrogen),
                 mean = mean(nitrogen),
                 sd = sd(nitrogen),
                 se = sd/ sqrt(N)
)
nitrogenSE

#plot
tiff("nitrogen_concen.tiff", units="in", width=7, height=5, res=300)

nit <- ggplot(data = nitrogenSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black") +
  ylab("% Nitrogen")+
  theme_bw() + scale_color_manual(values = c("C"="#88d8b0",
                                  "M"="#ffcc5c",
                                  "S"="#ff6f69"))+
  my_theme + theme(legend.position = 'none')
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(nutrient$nitrogen_raw)
#p < 0.05 so it IS NOT normal
skewness(nutrient$nitrogen_raw, na.rm = TRUE)
#slight positive skew
nutrient$nitrogen_raw <- sqrt(nutrient$nitrogen_raw)
shapiro.test(nutrient$nitrogen_raw)
#p > 0.05 so it IS normal
mod10 <- lm(nitrogen_raw ~ Inoculation*Watering, nutrient)
nit_anova <- aov(mod10)
summary(nit_anova)

#look phophorous############
#create data table containing means, sd, se, and N
phosphorousSE <- ddply(nutrient, c("Inoculation", "Watering"),summarise,
                    N = length(phosphorous),
                    mean = mean(phosphorous),
                    sd = sd(phosphorous),
                    se = sd/ sqrt(N)
)
phosphorousSE

#plot
tiff("p_concentration.tiff", units="in", width=7, height=5, res=300)

phos <- ggplot(data = phosphorousSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "black")  +
  ylab("Phosphorous mg/g")+
  theme_bw() + scale_color_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69"))+
  my_theme
dev.off()
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(nutrient$phosphorous)
#p < 0.05 it IS NOT normal
skewness(nutrient$phosphorous)
#slightly positive
nutrient$phosphorous <- sqrt(nutrient$phosphorous)
shapiro.test(nutrient$phosphorous)
#p > 0.05 it IS normal
mod11 <- lm(phosphorous ~ Inoculation*Watering, nutrient)
phos_anova <- aov(mod11)
summary(phos_anova)

#check moisture treatments
moisture <- readxl::read_xlsx('moisture_9.30.xlsx', sheet = 2)

moistSE <- ddply(moisture, c("Inoculation", "Watering"),summarise,
                       N = length(Moisture),
                       mean = mean(Moisture),
                       sd = sd(Moisture),
                       se = sd/ sqrt(N)
)
moistSE

#plot
tiff("moisture.tiff", units="in", width=7, height=5, res=300)

moist <- ggplot(moistSE, aes(x=Inoculation, y=mean, fill=Watering)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  ylab("Moisture %")+
  theme_bw()
moist + scale_fill_manual(values = c("C"="#88d8b0",
                                    "M"="#ffcc5c",
                                    "S"="#ff6f69"))+
  my_theme
dev.off()

#multiplots
tiff("biomass_multi.tiff", units="in", width=14, height=5, res=300)

figure <- ggarrange(sb, rb,
                    labels = c("A", "B"),
                    ncol = 2, nrow = 1)
figure

dev.off()

tiff("water_content_multi.tiff", units="in", width=14, height=5, res=300)

figure2 <- ggarrange(swconcen, rwconcen,
                    labels = c("A", "B"),
                    ncol = 2, nrow = 1)
figure2

dev.off()

tiff("pn_sc.tiff", units="in", width=14, height=5, res=300)

figure3 <- ggarrange(sc, pn,
                    labels = c("A", "B"),
                    ncol = 2, nrow = 1)
figure3

dev.off()

tiff("nutrient_multi.tiff", units="in", width=14, height=5, res=300)

figure4 <- ggarrange(nit, phos,
                    labels = c("A", "B"),
                    ncol = 2, nrow = 1)
figure4

dev.off()

#look at colonization######################
#load data
colonization <- readxl::read_xlsx('Colonization2.xlsx', sheet = 2)

#create data table containing means, sd, se, and N
colonizationSE <- ddply(colonization, c("Inoculation","Watering"), summarise,
                    N = length(Colonization),
                    mean = mean(Colonization),
                    sd = sd(Colonization),
                    se = sd/ sqrt(N)
)
colonizationSE


#plot
col <- ggplot(colonizationSE, aes(x=Inoculation, y=mean, fill=Watering)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  ylab("% Colonization")+
  theme_bw()
col + scale_fill_manual(values = c("C"="#88d8b0",
                                     "M"="#ffcc5c",
                                     "S"="#ff6f69"))
#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(colonization$raw)
#normal
col_mod <- lm(raw ~ Watering, colonization)
col_anova <- aov(col_mod)
summary(col_anova)

col_anova[["model"]][["Watering"]] = as.factor(col_anova[["model"]][["Watering"]])
tuk <- glht(col_anova, mcp(Watering = "Tukey")) 
#create cld object so we can pull out letters for better looking graph
tuk.cld <- cld(tuk)
#check it out
plot(tuk.cld)
print(tuk.cld)
#now to put it into ggplot format...
#create annotation
annotate <- data.frame(
  x = c(.7,1,1.3),
  y = 40,
  label = c("a" , "a", "a"))

col + scale_fill_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69")) + 
      geom_text(data = annotate, aes(x = x, y = y, label = label), inherit.aes = FALSE)

#vesicles#########
#create data table containing means, sd, se, and N
vesiclesSE <- ddply(colonization, c("Inoculation","Watering"), summarise,
                        N = length(v_per),
                        mean = mean(v_per),
                        sd = sd(v_per),
                        se = sd/ sqrt(N)
)
vesiclesSE


#plot
ves <- ggplot(vesiclesSE, aes(x=Inoculation, y=mean, fill=Watering)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  ylab("% Vesicles")+
  theme_bw()

ves + scale_fill_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69"))
#statistical analysis ~ LM, equal variance, and ANOVA
#lm
mod12 <- lm(v_raw ~ Watering, colonization)
plot(mod12)
#normal
#ANOVA
v_anova <- aov(mod12)
summary(v_anova)

v_anova[["model"]][["Watering"]] = as.factor(v_anova[["model"]][["Watering"]])
tuk <- glht(v_anova, mcp(Watering = "Tukey")) 
#create cld object so we can pull out letters for better looking graph
tuk.cld <- cld(tuk)
#check it out
plot(tuk.cld)
print(tuk.cld)
#now to put it into ggplot format...
#create annotation
annotate <- data.frame(
  x = c(.7,1,1.3),
  y = 40,
  label = c("a" , "a", "a"))

ves + scale_fill_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69")) + 
  geom_text(data = annotate, aes(x = x, y = y, label = label), inherit.aes = FALSE)

#arbuscules###########
#create data table containing means, sd, se, and N
arbusculesSE <- ddply(colonization, c("Inoculation","Watering"), summarise,
                    N = length(a_per),
                    mean = mean(a_per),
                    sd = sd(a_per),
                    se = sd/ sqrt(N)
)
arbusculesSE


#plot
arb <- ggplot(arbusculesSE, aes(x=Inoculation, y=mean, fill=Watering)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  ylab("% Arbuscules")+
  theme_bw()

arb + scale_fill_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69"))
#statistical analysis ~ LM, equal variance, and ANOVA
#lm
mod13 <- lm(a_raw ~ Watering, colonization)
plot(mod13)
#normal
#ANOVA
a_anova <- aov(mod13)
summary(a_anova)

a_anova[["model"]][["Watering"]] = as.factor(a_anova[["model"]][["Watering"]])
tuk <- glht(a_anova, mcp(Watering = "Tukey")) 
summary(tuk)
#create cld object so we can pull out letters for better looking graph
tuk.cld <- cld(tuk)
#check it out
plot(tuk.cld)
print(tuk.cld)
#now to put it into ggplot format...
#create annotation
annotate <- data.frame(
  x = c(.7,1,1.3),
  y = 4,
  label = c("a" , "a", "a"))

arb + scale_fill_manual(values = c("C"="#88d8b0",
                                   "M"="#ffcc5c",
                                   "S"="#ff6f69")) + 
  geom_text(data = annotate, aes(x = x, y = y, label = label), inherit.aes = FALSE)

nutrients3 <- readxl::read_xlsx('NP_measurement.xlsx', sheet = 3)

#total nitrogen###################
#create data table containing means, sd, se, and N
nitrogenSE <- ddply(nutrients3, c("Inoculation", "Watering"),summarise,
                    N = length(T_N),
                    mean = mean(T_N),
                    sd = sd(T_N),
                    se = sd/ sqrt(N)
)
nitrogenSE

#plot
nit <- ggplot(data = nitrogenSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "grey") +
  ylab("Total Nitrogen")+
  theme_grey()
nit + scale_color_manual(values = c("C"="#88d8b0",
                                    "M"="#ffcc5c",
                                    "S"="#ff6f69"))


#statistical analysis ~ LM, equal variance, and ANOVA
#Shapiro-Wilk test
shapiro.test(nutrients3$T_N)
#normal
mod14 <- lm(T_N ~ Inoculation*Watering, nutrients3)
nit_anova <- aov(mod14)
summary(nit_anova)

#phosphorous##############
#create data table containing means, sd, se, and N
phosphorousSE <- ddply(nutrients3, c("Inoculation", "Watering"),summarise,
                    N = length(T_P),
                    mean = mean(T_P),
                    sd = sd(T_P),
                    se = sd/ sqrt(N)
)
phosphorousSE

#plot
phos <- ggplot(data = phosphorousSE, aes(x=Inoculation, y=mean))+
  geom_line(size = 1.2, aes(group = Watering, color = Watering))+
  geom_point(size = 4, aes(color = Watering, shape = Inoculation))+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9), color = "grey") +
  ylab("Total Phosphorous")+
  theme_grey()
phos + scale_color_manual(values = c("C"="#88d8b0",
                                    "M"="#ffcc5c",
                                    "S"="#ff6f69"))


#statistical analysis ~ LM, equal variance, and ANOVA
mod15 <- lm(T_P ~ Inoculation*Watering, nutrients3)
plot(mod15)
#normal
phos_anova <- aov(mod15)
summary(phos_anova)