######## Shepherd et al., 2022 Peat moss inoculation depends on water table depth: Belowground community analysis
rm(list=ls())
library(tidyverse); library(ggplot2); library(vegan); library(phyloseq)
library(ggpubr); library(FD); library(Matching); library(plotrix)

# setwd("~/Documents/Dutch_innoculation_exp/Final_Code_data") ### set to your working directory to data location
setwd("C:/Users/k2256287/OneDrive - King's College London/Documents/PhD_work/Inoc_study/Repository_data/PURE_round1")


mic <- read.csv("Shepherd_Inoc_Micro_data.csv") ### microbial data: produced through DADA2 pipeline

### sorting out data
mic$number <- 1:nrow(mic)
mic$asv <- "asv"
mic %>% mutate(asv_id = paste(asv, number))-> mic
ID <- mic$asv_id
rownames(mic) <- ID

mic %>% dplyr::select(X1.1A:X5.4D) -> mic_ab
names(mic_ab) = gsub(pattern = "X", replacement = "", x = names(mic_ab))
remove(ID)
mic_ab ### microbial abundance matrix

read.csv("Shepherd_Inoc_Micro_trts.csv")  %>% unite(Sample_names, c("Epi_code", "Sample_code"), 
                    remove = FALSE, sep = "") %>%
  column_to_rownames("Sample_names") %>% 
  unite(Treatment, c("Water_table", "Inoculation"), remove = FALSE) %>%
  unite(Treatment_full, c("Water_table", "Inoculation", "Days"), remove = FALSE) -> micro_trt ### treatment data


mic %>% dplyr::select(Kingdom:Species) -> mic_phylo ### 16S phylogeny

tax_mat <- as.matrix(mic_phylo)
remove(mic)

##### NMDS of the microbial recovery #######

micro_trt %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_all

trt_all$Tube_code -> Tube_code_all
mic_ab %>% dplyr::select(Tube_code_all) -> mic_ab_all
asv_mat_all <- as.matrix(mic_ab_all)

ASV_all = otu_table(asv_mat_all, taxa_are_rows = TRUE)
TAX_all = tax_table(tax_mat)
trts_all = sample_data(trt_all)

phylo_all <- phyloseq(ASV_all, TAX_all, trts_all)

total = median(sample_sums(phylo_all)) ### Normalising reads for use
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_all = transform_sample_counts(phylo_all , standf)


### ordination
phylo_all.ord <- ordinate(phylo_all, "NMDS", "bray") ## stress 0.186
#### change names on this
plot_ordination(phylo_all, phylo_all.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_all, method="bray") ~ Inoculation*Water_table*Days,
       data = trt_all, strata = trt_all$Epi_code)

### trialling different structures: all good doesn't change interpretation
adonis(distance(phylo_all, method="bray") ~ Inoculation*Days*Water_table,
       data = trt_all, strata = trt_all$Epi_code)
adonis(distance(phylo_all, method="bray") ~ Days*Water_table*Inoculation,
       data = trt_all, strata = trt_all$Epi_code)
adonis(distance(phylo_all, method="bray") ~ Days*Inoculation*Water_table,
       data = trt_all, strata = trt_all$Epi_code)
adonis(distance(phylo_all, method="bray") ~ Water_table*Days*Inoculation,
       data = trt_all, strata = trt_all$Epi_code)
adonis(distance(phylo_all, method="bray") ~ Water_table*Inoculation*Days,
       data = trt_all, strata = trt_all$Epi_code)

as.data.frame(phylo_all.ord$points) -> all_pnts
cbind(all_pnts, trt_all) %>% mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  dplyr::select(MDS1, MDS2, Inoculation, Water_table, Days, Time) %>%
  mutate_if(is.character, str_replace_all, pattern="10_days", replacement="10 days") %>%
  mutate_if(is.character, str_replace_all, pattern="1.5_months", replacement="1.5 months") %>%
  mutate_if(is.character, str_replace_all, pattern="4_months", replacement="4 months") -> pnts_all2

pnts_all2$Time <- factor(pnts_all2$Time, levels=c("10 days", "1.5 months", "4 months"))

### Splitting into individual components

pnts_all2 %>% dplyr::filter(Time == '10 days') -> all_10d
alld10_plot <- ggplot() +
  geom_point(aes(x=MDS1, y=MDS2, colour=Water_table, shape=Inoculation), data=all_10d,
             size=2) +
  scale_fill_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) + xlab("NMDS1") + ylab("NMDS2") + theme_bw() +
  xlim(-0.7, 1.3) + ylim(-0.65, 0.7) +
  guides(color=guide_legend("Water table"), fill = "none", shape = guide_legend("Inoculated")) +
  ggtitle("(a) 10 days") +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) 

pnts_all2 %>% dplyr::filter(Time == '1.5 months') -> all_m15
allm15_plot <- ggplot() +
  geom_point(aes(x=MDS1, y=MDS2, colour=Water_table, shape=Inoculation), data=all_m15,
             size=2) +
  scale_fill_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) + xlab("NMDS1") + ylab("NMDS2") + theme_bw()+
  xlim(-0.7, 1.3) + ylim(-0.65, 0.7)  +
  guides(color=guide_legend("Water table"), fill = "none", shape = guide_legend("Inoculated")) +
  ggtitle("(b) 35 days") +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

pnts_all2 %>% dplyr::filter(Time == '4 months') -> all_m4
allm4_plot <- ggplot() +
  geom_point(aes(x=MDS1, y=MDS2, colour=Water_table, shape=Inoculation), data=all_m4,
             size=2) +
  scale_fill_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) + xlab("NMDS1") + ylab("NMDS2") + theme_bw()+
  xlim(-0.7, 1.3) + ylim(-0.65, 0.7) + 
  guides(color=guide_legend("Water table"), fill = "none", shape = guide_legend("Inoculated")) +
  ggtitle("(c) 112 days")  +
  theme(axis.text = element_text(size=12), 
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) 


#### Part 2: community proportion of the 10 most abundant phylum

### composition over time: top 10 phylum
cbind(mic_ab, mic_phylo) %>% dplyr::select('1.1A':'5.4D', Phylum) %>% group_by(Phylum) %>%
  summarise(across(c('1.1A':'5.4D'), ~ sum(.x, na.rm = TRUE))) -> mic_df2

### find total most abundant phyla
mic_df2 %>% dplyr::select(-Phylum) -> phy_df
as.data.frame(rowSums(phy_df)) -> tot_counts_phy
mic_df2 %>% dplyr::select(Phylum) -> phy_list
cbind(phy_list, tot_counts_phy) %>% dplyr::rename(total = 'rowSums(phy_df)') -> phy_sort
phy_sort$total <- as.numeric(phy_sort$total)
phy_sort %>% arrange(desc(total)) #### order of abundance, top 10 selected from here
sum(phy_sort$total) # 5,231,567 total asv count

### top 10 phylum
tphy <- c("Proteobacteria", "Acidobacteria", "Actinobacteria", "Verrucomicrobia", "Bacteroidetes",
           "Planctomycetes", "Euryarchaeota", "Thaumarchaeota", "Cyanobacteria", "WPS-2")

### calculate the proportion in each community
micro_trt%>% dplyr::select(Inoculation:Time) %>% unite(Plot, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> mtrt2

as.data.frame(colSums(phy_df)) %>% rename(Total_asv = 'colSums(phy_df)') %>% rownames_to_column("Plot") %>%
  left_join(mtrt2, by="Plot") %>% group_by(Inoculation, Water_table, Time) %>%
  summarise(across(c('Total_asv'), ~ sum(.x, na.rm = TRUE))) %>% unite(S_code, c("Inoculation", "Water_table", "Time")) -> comm_tot

mic_df2 %>% tidyr::gather("Plot", "ASVs", 2:61) %>% filter(Phylum %in% tphy) %>%
  left_join(mtrt2, by="Plot") %>% 
  group_by(Inoculation, Water_table, Time, Phylum) %>% 
  summarise(across(c('ASVs'), ~ sum(.x, na.rm = TRUE))) %>%
  unite(S_code, c("Inoculation", "Water_table", "Time"), remove=FALSE) %>%
  left_join(comm_tot, by="S_code") %>% dplyr::select(-S_code) %>%
  mutate(prop = ASVs/Total_asv) %>% ungroup() -> lng_mic 

lng_mic %>% dplyr::select(Inoculation, Water_table) %>%  
  mutate_if(is.character, str_replace_all, pattern = "N", replacement = "Uninoculated") %>% 
  mutate_if(is.character, str_replace_all, pattern = "Y", replacement = "Inoculated") %>%
  mutate_if(is.character, str_replace_all, pattern = "H", replacement = "High WT") %>%
  mutate_if(is.character, str_replace_all, pattern = "L", replacement = "Low WT") -> lng_trts

lng_mic %>% dplyr::select(-Inoculation, -Water_table) -> lng_mic2
cbind(lng_trts, lng_mic2) %>% unite(Treatment, c("Inoculation", "Water_table"), remove=FALSE, sep=" + ") %>%
  mutate_if(is.character, str_replace_all, pattern = "1.5_months", replacement = "35 days") %>%
  mutate_if(is.character, str_replace_all, pattern = "10_days", replacement = "10 days") %>%
  mutate_if(is.character, str_replace_all, pattern = "4_months", replacement = "112 days") -> lng_mic3

lng_mic3$Time <- factor(lng_mic3$Time, levels=c("10 days", "35 days", "112 days"))
lng_mic3$Phylum <- factor(lng_mic3$Phylum, levels=c("WPS-2", "Cyanobacteria", "Thaumarchaeota", "Euryarchaeota", "Planctomycetes", 
                                                    "Bacteroidetes", "Verrucomicrobia", "Actinobacteria", "Acidobacteria", "Proteobacteria"))
#### Split phylum plot

lng_mic3 %>% dplyr::filter(Inoculation == "Uninoculated" & Water_table == "High WT") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop)-> Un_Hg


Un_Hg %>% group_by(Time, Treatment, Inoculation, Water_table) %>% summarise(total_prop = sum(prop)) %>% dplyr::mutate(prop = 1-total_prop) %>%
  dplyr::select(-total_prop) %>% dplyr::mutate(Phylum = "Other") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> other_Un_Hg

rbind(Un_Hg, other_Un_Hg)  %>%
  mutate_if(is.factor, str_replace_all, pattern = "35 days", replacement = "35") %>%
  mutate_if(is.character, str_replace_all, pattern = "10 days", replacement = "10") %>%
  mutate_if(is.character, str_replace_all, pattern = "112 days", replacement = "112") -> Un_Hg2

Un_Hg2$Time <- factor(Un_Hg2$Time, levels = c("10", "35", "112"))

Un_Hg2$Phylum <- factor(Un_Hg2$Phylum, levels=c("Acidobacteria", "Actinobacteria", "Bacteroidetes", "Cyanobacteria",
                                                "Euryarchaeota", "Planctomycetes", "Proteobacteria", "Thaumarchaeota",
                                                "Verrucomicrobia", "WPS-2", "Other"))

prop_Un_hg <- ggplot() +
  geom_bar(aes(x=Time, y=prop, fill=Phylum), data=Un_Hg2, position = "stack",
           stat="identity") +
  ylim(0, 1) +
  theme_bw() +
  ylab("Proportion of all ASVs") +
  xlab("Days since inoculation") +
  scale_fill_brewer("Phylum", palette = "Paired") +
  theme(strip.background = element_blank(),
        strip.placement = "inside") +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  ggtitle("(e) Uninoculated + \n     High WT")+
  scale_y_continuous(breaks = c(0, 0.5, 1))

lng_mic3 %>% dplyr::filter(Inoculation == "Uninoculated" & Water_table == "Low WT") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop)-> Un_Lw

Un_Lw %>% group_by(Time, Treatment, Inoculation, Water_table) %>% summarise(total_prop = sum(prop)) %>% dplyr::mutate(prop = 1-total_prop) %>%
  dplyr::select(-total_prop) %>% dplyr::mutate(Phylum = "Other") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> other_Un_Lw

rbind(Un_Lw, other_Un_Lw)  %>%
  mutate_if(is.factor, str_replace_all, pattern = "35 days", replacement = "35") %>%
  mutate_if(is.character, str_replace_all, pattern = "10 days", replacement = "10") %>%
  mutate_if(is.character, str_replace_all, pattern = "112 days", replacement = "112") -> Un_Lw2

Un_Lw2$Time <- factor(Un_Lw2$Time, levels = c("10", "35", "112"))

Un_Lw2$Phylum <- factor(Un_Lw2$Phylum, levels=c("Acidobacteria", "Actinobacteria", "Bacteroidetes", "Cyanobacteria",
                                                "Euryarchaeota", "Planctomycetes", "Proteobacteria", "Thaumarchaeota",
                                                "Verrucomicrobia", "WPS-2", "Other"))

prop_Un_lw <- ggplot() +
  geom_bar(aes(x=Time, y=prop, fill=Phylum), data=Un_Lw2, position = "stack",
           stat="identity") +
  ylim(0, 1) +
  theme_bw() +
  ylab("Proportion of all ASVs") +
  xlab("Days since inoculation") +
  scale_fill_brewer("Phylum", palette = "Paired") +
  theme(strip.background = element_blank(),
        strip.placement = "inside") +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  ggtitle("(d) Uninoculated + \n     Low WT")+
  scale_y_continuous(breaks = c(0, 0.5, 1))

lng_mic3 %>% dplyr::filter(Inoculation == "Inoculated" & Water_table == "Low WT")  %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> In_Lw

In_Lw %>% group_by(Time, Treatment, Inoculation, Water_table) %>% summarise(total_prop = sum(prop)) %>% dplyr::mutate(prop = 1-total_prop) %>%
  dplyr::select(-total_prop) %>% dplyr::mutate(Phylum = "Other") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> other_In_Lw

rbind(In_Lw, other_In_Lw)  %>%
  mutate_if(is.factor, str_replace_all, pattern = "35 days", replacement = "35") %>%
  mutate_if(is.character, str_replace_all, pattern = "10 days", replacement = "10") %>%
  mutate_if(is.character, str_replace_all, pattern = "112 days", replacement = "112") -> In_Lw2

In_Lw2$Time <- factor(In_Lw2$Time, levels = c("10", "35", "112"))

In_Lw2$Phylum <- factor(In_Lw2$Phylum, levels=c("Acidobacteria", "Actinobacteria", "Bacteroidetes", "Cyanobacteria",
                                                "Euryarchaeota", "Planctomycetes", "Proteobacteria", "Thaumarchaeota",
                                                "Verrucomicrobia", "WPS-2", "Other"))

prop_In_lw <- ggplot() +
  geom_bar(aes(x=Time, y=prop, fill=Phylum), data=In_Lw2, position = "stack",
           stat="identity") +
  ylim(0, 1) +
  theme_bw() +
  ylab("Proportion of all ASVs") +
  xlab("Days since inoculation") +
  scale_fill_brewer("Phylum", palette = "Paired") +
  theme(strip.background = element_blank(),
        strip.placement = "inside") +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  ggtitle("(f) Inoculated + \n    Low WT")+
  scale_y_continuous(breaks = c(0, 0.5, 1))

lng_mic3 %>% dplyr::filter(Inoculation == "Inoculated" & Water_table == "High WT")  %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> In_Hg

In_Hg %>% group_by(Time, Treatment, Inoculation, Water_table) %>% summarise(total_prop = sum(prop)) %>% dplyr::mutate(prop = 1-total_prop) %>%
  dplyr::select(-total_prop) %>% dplyr::mutate(Phylum = "Other") %>%
  dplyr::select(Treatment, Inoculation, Water_table, Time, Phylum, prop) -> other_In_Hg

rbind(In_Hg, other_In_Hg) %>%
  mutate_if(is.factor, str_replace_all, pattern = "35 days", replacement = "35") %>%
  mutate_if(is.character, str_replace_all, pattern = "10 days", replacement = "10") %>%
  mutate_if(is.character, str_replace_all, pattern = "112 days", replacement = "112") -> In_Hg2

In_Hg2$Time <- factor(In_Hg2$Time, levels = c("10", "35", "112"))

In_Hg2$Phylum <- factor(In_Hg2$Phylum, levels=c("Acidobacteria", "Actinobacteria", "Bacteroidetes", "Cyanobacteria",
                                                    "Euryarchaeota", "Planctomycetes", "Proteobacteria", "Thaumarchaeota",
                                                    "Verrucomicrobia", "WPS-2", "Other"))

prop_In_Hg <- ggplot() +
  geom_bar(aes(x=Time, y=prop, fill=Phylum), data=In_Hg2, position = "stack",
           stat="identity") +
  ylim(0, 1) +
  theme_bw() +
  ylab("") +
  xlab("Days since inoculation") +
  scale_fill_brewer("Phylum", palette = "Paired") +
  theme(strip.background = element_blank(),
        strip.placement = "inside") +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  ggtitle("(g) Inoculated + \n     High WT") +
  scale_y_continuous(breaks = c(0, 0.5, 1))

library(grid)

phylo_figure <- (ggarrange(prop_Un_lw + rremove("ylab") + rremove("xlab"),
                           prop_Un_hg + rremove("ylab") + rremove("xlab"), 
                           prop_In_lw + rremove("ylab") + rremove("xlab"),
                           prop_In_Hg + rremove("ylab") + rremove("xlab"), 
                           labels = NULL, ncol=4, common.legend = TRUE, legend="right",
                           align="hv", 
                           font.label = list(size = 14, color = "black", face = "bold", family = NULL, position = "top")))

remove(phy_sort, phy_list, comm_tot, tot_counts_phy, lng_mic2, lng_mic, lng_trts)

### Main figure (Figure 3)
fig_top2 <- ggarrange(alld10_plot, allm15_plot, allm4_plot, ncol=3, common.legend=TRUE, legend = "right")
fig_bot2 <- annotate_figure(phylo_figure, left = textGrob("Proportion of all ASVs", rot = 90, 
                                                          vjust = 1, gp = gpar(cex = 1.4)),
                           bottom = textGrob("Days since inoculation", gp = gpar(cex = 1.4)))
# ggarrange(fig_top, fig_bot, nrow=2, ncol=1)

library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow = 100, ncol = 100)))
define_region <- function(row, col){
  viewport(layout.pos.row = row, layout.pos.col = col)
} 

print(fig_top2, vp = define_region(row = 1:33, col = 1:95))   # Span over two columns
print(fig_bot2, vp = define_region(row = 37:100, col = 1:100))
# export 1200 x 800

##### Producing supplementary figures 


######## Individual daily NMDS plots with water table split
#### 10 days: HWT 
micro_trt %>% filter(Time == "10_days" & Water_table == "H") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_day10_hwt

trt_day10_hwt$Tube_code -> Tube_code_day10_hwt
mic_ab %>% dplyr::select(Tube_code_day10_hwt) -> mic_ab_day10_hwt
asv_mat_day10_hwt <- as.matrix(mic_ab_day10_hwt)

ASV_day10_hwt = otu_table(asv_mat_day10_hwt, taxa_are_rows = TRUE)
TAX_day10_hwt = tax_table(tax_mat)
trts_day10_hwt = sample_data(trt_day10_hwt)

phylo_day10_hwt <- phyloseq(ASV_day10_hwt, TAX_day10_hwt, trts_day10_hwt)

### normalise number of reads
total = median(sample_sums(phylo_day10_hwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_day10_hwt  = transform_sample_counts(phylo_day10_hwt , standf)

### ordination
phylo_day10_hwt.ord <- ordinate(phylo_day10_hwt, "NMDS", "bray") ## stress 0.07
#### change names on this
plot_ordination(phylo_day10_hwt, phylo_day10_hwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_day10_hwt, method="bray") ~ Inoculation,
       data = trt_day10_hwt) ### no differences here

#### 10 days: LWT 
micro_trt %>% filter(Time == "10_days" & Water_table == "L") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_day10_lwt

trt_day10_lwt$Tube_code -> Tube_code_day10_lwt
mic_ab %>% dplyr::select(Tube_code_day10_lwt) -> mic_ab_day10_lwt
asv_mat_day10_lwt <- as.matrix(mic_ab_day10_lwt)

ASV_day10_lwt = otu_table(asv_mat_day10_lwt, taxa_are_rows = TRUE)
TAX_day10_lwt = tax_table(tax_mat)
trts_day10_lwt = sample_data(trt_day10_lwt)

phylo_day10_lwt <- phyloseq(ASV_day10_lwt, TAX_day10_lwt, trts_day10_lwt)

### normalise number of reads
total = median(sample_sums(phylo_day10_lwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_day10_lwt  = transform_sample_counts(phylo_day10_lwt , standf)

### ordination
phylo_day10_lwt.ord <- ordinate(phylo_day10_lwt, "NMDS", "bray") ## stress 0.06
#### change names on this
plot_ordination(phylo_day10_lwt, phylo_day10_lwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_day10_lwt, method="bray") ~ Inoculation,
       data = trt_day10_lwt) ### no differences here


#### 1.5 months 
micro_trt %>% filter(Time == "1.5_months" & Water_table == "H") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_mnth1.5_hwt

trt_mnth1.5_hwt$Tube_code -> Tube_code_mnth1.5_hwt
mic_ab %>% dplyr::select(Tube_code_mnth1.5_hwt) -> mic_ab_mnth1.5_hwt
asv_mat_mnth1.5_hwt <- as.matrix(mic_ab_mnth1.5_hwt)

ASV_mnth1.5_hwt = otu_table(asv_mat_mnth1.5_hwt, taxa_are_rows = TRUE)
TAX_mnth1.5_hwt = tax_table(tax_mat)
trts_mnth1.5_hwt = sample_data(trt_mnth1.5_hwt)

phylo_mnth1.5_hwt <- phyloseq(ASV_mnth1.5_hwt, TAX_mnth1.5_hwt, trts_mnth1.5_hwt)

### normalise number of reads
total = median(sample_sums(phylo_mnth1.5_hwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_mnth1.5_hwt  = transform_sample_counts(phylo_mnth1.5_hwt , standf)

### ordination
phylo_mnth1.5_hwt.ord <- ordinate(phylo_mnth1.5_hwt, "NMDS", "bray") ## stress 0.07
#### change names on this
plot_ordination(phylo_mnth1.5_hwt, phylo_mnth1.5_hwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_mnth1.5_hwt, method="bray") ~ Inoculation,
       data = trt_mnth1.5_hwt) ### no differences here

## 1.5 months low water table
micro_trt %>% filter(Time == "1.5_months" & Water_table == "L") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_mnth1.5_lwt

trt_mnth1.5_lwt$Tube_code -> Tube_code_mnth1.5_lwt
mic_ab %>% dplyr::select(Tube_code_mnth1.5_lwt) -> mic_ab_mnth1.5_lwt
asv_mat_mnth1.5_lwt <- as.matrix(mic_ab_mnth1.5_lwt)

ASV_mnth1.5_lwt = otu_table(asv_mat_mnth1.5_lwt, taxa_are_rows = TRUE)
TAX_mnth1.5_lwt = tax_table(tax_mat)
trts_mnth1.5_lwt = sample_data(trt_mnth1.5_lwt)

phylo_mnth1.5_lwt <- phyloseq(ASV_mnth1.5_lwt, TAX_mnth1.5_lwt, trts_mnth1.5_lwt)

### normalise number of reads
total = median(sample_sums(phylo_mnth1.5_lwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_mnth1.5_lwt  = transform_sample_counts(phylo_mnth1.5_lwt , standf)

### ordination
phylo_mnth1.5_lwt.ord <- ordinate(phylo_mnth1.5_lwt, "NMDS", "bray") ## stress 0.07
#### change names on this
plot_ordination(phylo_mnth1.5_lwt, phylo_mnth1.5_lwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_mnth1.5_lwt, method="bray") ~ Inoculation,
       data = trt_mnth1.5_lwt) ### no differences here


#### 4 months 
micro_trt %>% filter(Time == "4_months" & Water_table == "H") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_mnth4_hwt

trt_mnth4_hwt$Tube_code -> Tube_code_mnth4_hwt
mic_ab %>% dplyr::select(Tube_code_mnth4_hwt) -> mic_ab_mnth4_hwt
asv_mat_mnth4_hwt <- as.matrix(mic_ab_mnth4_hwt)

ASV_mnth4_hwt = otu_table(asv_mat_mnth4_hwt, taxa_are_rows = TRUE)
TAX_mnth4_hwt = tax_table(tax_mat)
trts_mnth4_hwt = sample_data(trt_mnth4_hwt)

phylo_mnth4_hwt <- phyloseq(ASV_mnth4_hwt, TAX_mnth4_hwt, trts_mnth4_hwt)

### normalise number of reads
total = median(sample_sums(phylo_mnth4_hwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_mnth4_hwt  = transform_sample_counts(phylo_mnth4_hwt , standf)

### ordination
phylo_mnth4_hwt.ord <- ordinate(phylo_mnth4_hwt, "NMDS", "bray") ## stress 0.07
#### change names on this
plot_ordination(phylo_mnth4_hwt, phylo_mnth4_hwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_mnth4_hwt, method="bray") ~ Inoculation,
       data = trt_mnth4_hwt) ### weak effect here

micro_trt %>% filter(Time == "4_months" & Water_table == "L") %>% 
  unite(Tube_code, c("Epi_code", "Sample_code"), sep="", remove=FALSE) -> trt_mnth4_lwt

trt_mnth4_lwt$Tube_code -> Tube_code_mnth4_lwt
mic_ab %>% dplyr::select(Tube_code_mnth4_lwt) -> mic_ab_mnth4_lwt
asv_mat_mnth4_lwt <- as.matrix(mic_ab_mnth4_lwt)

ASV_mnth4_lwt = otu_table(asv_mat_mnth4_lwt, taxa_are_rows = TRUE)
TAX_mnth4_lwt = tax_table(tax_mat)
trts_mnth4_lwt = sample_data(trt_mnth4_lwt)

phylo_mnth4_lwt <- phyloseq(ASV_mnth4_lwt, TAX_mnth4_lwt, trts_mnth4_lwt)

### normalise number of reads
total = median(sample_sums(phylo_mnth4_lwt ))
standf = function(x, t=total) round(t * (x / sum(x)))
phylo_mnth4_lwt  = transform_sample_counts(phylo_mnth4_lwt , standf)

### ordination
phylo_mnth4_lwt.ord <- ordinate(phylo_mnth4_lwt, "NMDS", "bray") ## stress 0.07
#### change names on this
plot_ordination(phylo_mnth4_lwt, phylo_mnth4_lwt.ord, type="samples", 
                shape="Inoculation") + geom_point(size=2) 
adonis(distance(phylo_mnth4_lwt, method="bray") ~ Inoculation,
       data = trt_mnth4_lwt) ### no differences here

###### Plotting microbial NMDS

as.data.frame(phylo_day10_hwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "10 days") -> micr0_d10_pnts_hwt
as.data.frame(phylo_mnth1.5_hwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "1.5 months") -> micr0_mnth1.5_pnts_hwt
as.data.frame(phylo_mnth4_hwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "4 months") -> micr0_mnth4_pnts_hwt

as.data.frame(phylo_day10_lwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "10 days") -> micr0_d10_pnts_lwt
as.data.frame(phylo_mnth1.5_lwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "1.5 months") -> micr0_mnth1.5_pnts_lwt
as.data.frame(phylo_mnth4_lwt.ord$points) %>% rownames_to_column("Mesocosm") %>%
  dplyr::mutate(Time = "4 months") -> micr0_mnth4_pnts_lwt

rbind(micr0_d10_pnts_hwt, micr0_mnth1.5_pnts_hwt, micr0_mnth4_pnts_hwt,
      micr0_d10_pnts_lwt, micr0_mnth1.5_pnts_lwt, micr0_mnth4_pnts_lwt) -> micro_pnts_wt


micro_trt %>% dplyr::select(Inoculation, Water_table, Epi_code, Sample_code) %>%
  unite(Mesocosm, c("Epi_code", "Sample_code"), sep="") %>%
  right_join(micro_pnts_wt, by="Mesocosm") -> micro_nmds2_wt

micro_nmds2_wt %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") -> micro_nmds3_wt

micro_nmds3_wt$Time <- factor(micro_nmds3_wt$Time, levels=c("10 days", "1.5 months", "4 months"))

ggplot() +
  geom_point(aes(x=MDS1, y=MDS2, colour=Water_table, shape=Inoculation), data=micro_nmds3_wt,
             size=2) +
  scale_fill_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) +
  facet_wrap(~Water_table + Time, ncol=3, strip.position = "top") +
  xlab("NMDS1") +
  ylab("NMDS2") +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text=element_text(face="bold", size=12)) 


###### Plotting individual phylum plus linear models of changing abundance
as.data.frame(colSums(phy_df)) %>% rename(Total_asv = 'colSums(phy_df)') %>% rownames_to_column("Plot") -> plot_asvs

mic_df2 %>% tidyr::gather("Plot", "ASVs", 2:61) %>% filter(Phylum %in% tphy) %>%
  left_join(mtrt2, by="Plot") %>% left_join(plot_asvs, by="Plot") %>%
  mutate(prop = ASVs/Total_asv) %>% dplyr::select(Phylum, prop, Plot, Inoculation, Water_table, Days, Time) %>%
  mutate_if(is.character, str_replace_all, pattern="N", replacement="No") %>%
  mutate_if(is.character, str_replace_all, pattern="Y", replacement="Yes") %>%
  mutate_if(is.character, str_replace_all, pattern="H", replacement="High") %>%
  mutate_if(is.character, str_replace_all, pattern="L", replacement="Low") %>%
  unite(Treatment, c("Inoculation", "Water_table"), remove=FALSE) -> plot_prop


plot_prop %>% group_by(Inoculation, Water_table, Treatment, Phylum, Time, Days) %>%
  summarise(across(c("prop"), ~ mean(.x, na.rm = TRUE))) -> prop_means

plot_prop %>% group_by(Inoculation, Water_table, Treatment, Phylum, Time, Days) %>%
  summarise(across(c("prop"), ~ std.error(.x, na.rm = TRUE))) %>%
  rename(ci_prop = prop) %>% ungroup() %>%
  dplyr::select(Treatment, ci_prop) %>% full_join(prop_means, by="Treatment") -> prop_values

prop_values$Time <- as.character(prop_values$Time)
prop_values %>%
  mutate_if(is.character, str_replace_all, pattern="10_days", replacement="10 days") %>%
  mutate_if(is.character, str_replace_all, pattern="1.5_months", replacement="35 days") %>%
  mutate_if(is.character, str_replace_all, pattern="4_months", replacement="112 days") -> prop_values2
prop_values2$Time <- factor(prop_values2$Time, levels=c("10 days", "35 days", "112 days"))

ggplot() +
  geom_pointrange(aes(x=Time, y=prop, ymax=prop+ci_prop,
                      ymin=prop-ci_prop, fill=Water_table, colour=Water_table, shape=Inoculation),
                  position = position_dodge(0.6), data=prop_values2) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  guides(color=guide_legend("Water table"), fill = "none") +
  ylab("Proportion of total ASVs") + xlab("Days since inoculation") +
  theme_bw() +
  scale_shape_manual(values=c(1, 19)) +
  facet_wrap(~Phylum, scales="free_y", ncol=3) +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_text(size=12))

plot_prop$Days_sc <- scale(plot_prop$Days)

mtrt2 %>% dplyr::select(Plot, Epi_code) %>% right_join(plot_prop, by="Plot") -> plot_prop2


### selection process: dredge top 6 AICc, simplify down removing all those that are nested, then average rest

## Acidobacteria
plot_prop2 %>% filter(Phylum == "Acidobacteria") -> phy1
acido_mod <- lme(log(prop) ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                   Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                 random = ~1|Epi_code, data=phy1, method="ML")
summary(acido_mod)
dredge(acido_mod, rank = "AICc")


acido_mod_bp1 <- lme(prop ~ Days_sc, random = ~1|Epi_code, data=phy1, method="ML")
summary(acido_mod_bp1)
r.squaredGLMM(acido_mod_bp1)
AICc(acido_mod_bp1)
check_model(acido_mod_bp1)


library(lme4)
acido.mod2 <-  lmer(prop ~  Days_sc + (1|Epi_code), data=phy1, REML=FALSE)
library(DHARMa)
simulationOutput <- simulateResiduals(fittedModel = acido.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 
### this is fine

## Actinobacteria
plot_prop2 %>% filter(Phylum == "Actinobacteria") -> phy2
actino_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                    Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                  random = ~1|Epi_code, data=phy2, method="ML")
summary(actino_mod)
dredge(actino_mod)
actino_dr <- dredge(actino_mod)

actino_mod1 <- lme(prop ~ Inoculation+ Days_sc, random = ~1|Epi_code, data=phy2, method="ML")
actino_mod2 <- lme(prop ~ Days_sc, random = ~1|Epi_code, data=phy2, method="ML")

AICc(actino_mod1)
r.squaredGLMM(actino_mod1)
r.squaredGLMM(actino_mod2)

actino_dr2 <- subset(actino_dr1, !nested(.))
summary(model.avg(actino_dr2, subset = delta <= 6))

actino_mod_bp <- lme(prop ~ Inoculation + Days_sc, 
                     random = ~1|Epi_code, data=phy2, method="ML")
check_model(actino_mod_bp)

actino.mod2 <-  lmer(prop ~  Days_sc + Inoculation + (1|Epi_code), data=phy2, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = actino.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

## Bacteroidetes
plot_prop2 %>% filter(Phylum == "Bacteroidetes") -> phy3
bacter_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                    Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                  random = ~1|Epi_code, data=phy3, method="ML")
summary(bacter_mod)
dredge(bacter_mod)
bacter_dr <- dredge(bacter_mod)
bacter_dr2 <- subset(bacter_dr, !nested(.))
summary(model.avg(bacter_dr2, subset = delta <= 6))

bacter_mod1 <- lme(prop ~ Inoculation+ Days_sc + Water_table + Inoculation:Water_table, random = ~1|Epi_code, data=phy3, method="ML")
bacter_mod2 <- lme(prop ~ Inoculation+ Days_sc + Water_table, random = ~1|Epi_code, data=phy3, method="ML")
bacter_mod3 <- lme(prop ~ Days_sc + Water_table, random = ~1|Epi_code, data=phy3, method="ML")
bacter_mod4 <- lme(prop ~ Inoculation+ Days_sc, random = ~1|Epi_code, data=phy3, method="ML")

r.squaredGLMM(bacter_mod1)
r.squaredGLMM(bacter_mod2)
r.squaredGLMM(bacter_mod3)
r.squaredGLMM(bacter_mod4)
AICc(bacter_mod1)

bacter_bpm <- lme((prop) ~ Inoculation + Days_sc + Water_table + Inoculation:Water_table, 
                  random = ~1|Epi_code, data=phy3, method="ML")
summary(bacter_bpm)

check_model(bacter_bpm)

bacter.mod2 <-  lmer(prop ~  Inoculation + Days_sc + Water_table + Inoculation:Water_table + (1|Epi_code), data=phy3, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = bacter.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 


## Cyanobacteria
plot_prop2 %>% filter(Phylum == "Cyanobacteria") -> phy4
cyano_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                   Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                 random = ~1|Epi_code, data=phy4, method="ML")
summary(cyano_mod)
dredge(cyano_mod)
cyano_dr <- dredge(cyano_mod)
cyano_dr2 <- subset(cyano_dr, !nested(.))
summary(model.avg(cyano_dr2, subset = delta <= 6))

cyano_mod_ss1 <- lme(prop ~ Inoculation + Days_sc + Water_table + Inoculation:Water_table + Inoculation:Days_sc, 
                     random = ~1|Plot, data=phy4, method="ML")
cyano_mod_ss2 <- lme(prop ~ Inoculation + Days_sc + Water_table + Inoculation:Water_table, 
                     random = ~1|Plot, data=phy4, method="ML")
cyano_mod_ss3 <- lme(prop ~ Inoculation + Days_sc + Water_table + Inoculation:Days_sc, 
                     random = ~1|Plot, data=phy4, method="ML")
cyano_mod_ss4 <- lme(prop ~ Inoculation + Days_sc + Water_table, 
                     random = ~1|Plot, data=phy4, method="ML")
cyano_mod_ss5 <- lme(prop ~ Inoculation + Days_sc + Inoculation:Days_sc, 
                     random = ~1|Plot, data=phy4, method="ML")
AICc(cyano_mod_ss1)
r.squaredGLMM(cyano_mod_ss1)
r.squaredGLMM(cyano_mod_ss2)
r.squaredGLMM(cyano_mod_ss3)
r.squaredGLMM(cyano_mod_ss4)
r.squaredGLMM(cyano_mod_ss5)

cyano_bpm <- lme((prop) ~ Inoculation + Days_sc + Water_table + Inoculation:Water_table + Inoculation:Days_sc, 
                 random = ~1|Epi_code, data=phy4, method="ML")
check_model(cyano_bpm)

cyano.mod2 <-  lmer(prop ~  Inoculation + Days_sc + Water_table + Inoculation:Water_table + Inoculation:Days_sc + (1|Epi_code), data=phy4, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = cyano.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

## Euryarchaeota
plot_prop2 %>% filter(Phylum == "Euryarchaeota") -> phy5
eury_mod <- lme(sqrt(prop) ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                  Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                random = ~1|Epi_code, data=phy5, method="ML")
summary(eury_mod)
dredge(eury_mod)
eury_mod_bp1 <- lme(prop ~ Days_sc, 
                    random = ~1|Epi_code, data=phy5, method="ML")
eury_mod_bp2 <- lme(prop ~ 1, 
                    random = ~1|Epi_code, data=phy5, method="ML")

eury_bpm <- lme((prop) ~ Days_sc, random = ~1|Epi_code, data=phy5, method="ML")
check_model(eury_bpm)
summary(eury_bpm)
AICc(eury_bpm)
r.squaredGLMM(eury_bpm)

eury.mod2 <-  lmer((prop) ~  Days_sc + (1|Epi_code), data=phy5, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = eury.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

## Planctomycetes
plot_prop2 %>% filter(Phylum == "Planctomycetes") -> phy6
plan_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                  Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                random = ~1|Epi_code, data=phy6, method="ML")
summary(plan_mod)
dredge(plan_mod)
plan_mod_bp1 <- lme(prop ~ Days_sc, 
                    random = ~1|Epi_code, data=phy6, method="ML")
summary(plan_mod_bp1)
r.squaredGLMM(plan_mod_bp1)
AICc(plan_mod_bp1)
check_model(plan_mod_bp1)

plan.mod2 <-  lmer(prop ~  Days_sc + (1|Epi_code), data=phy6, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = plan.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 


### Proteobacteria
plot_prop2 %>% filter(Phylum == "Proteobacteria") -> phy7
prot_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                  Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                random = ~1|Epi_code, data=phy7, method="ML")
summary(prot_mod)
dredge(prot_mod)
prot_mod_bp1 <- lme(prop ~ Days_sc, 
                    random = ~1|Epi_code, data=phy7, method="ML")
summary(prot_mod_bp1)
AICc(prot_mod_bp1)
r.squaredGLMM(prot_mod_bp1)

check_model(prot_mod_bp1)

prot.mod2 <-  lmer(prop ~  Days_sc + (1|Epi_code), data=phy7, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = prot.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

### Thaumarchaeota
plot_prop2 %>% filter(Phylum == "Thaumarchaeota") -> phy8
tham_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                  Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                random = ~1|Epi_code, data=phy8, method="ML")
summary(tham_mod)
dredge(tham_mod)
tham_mod_bp1 <- lme(prop ~ Days_sc, 
                    random = ~1|Epi_code, data=phy8, method="ML")
summary(tham_mod_bp1)
AICc(tham_mod_bp1)
r.squaredGLMM(tham_mod_bp1)
check_model(tham_mod_bp1)

tham.mod2 <-  lmer(prop ~  Days_sc + (1|Epi_code), data=phy8, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = tham.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

## Verrucomicrobia
plot_prop2 %>% filter(Phylum == "Verrucomicrobia") -> phy9
verr_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                  Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
                random = ~1|Epi_code, data=phy9, method="ML")
summary(verr_mod)
verr_dr <- dredge(verr_mod)
verr_dr2 <- subset(verr_dr, !nested(.))
summary(model.avg(verr_dr2, subset = delta <= 6))

vmod1 <- lme(prop ~  Days_sc + Water_table + Water_table:Days_sc, random = ~1|Epi_code, data=phy9, method="ML")
vmod2 <- lme(prop ~  Days_sc + Water_table + Inoculation, random = ~1|Epi_code, data=phy9, method="ML")
vmod3 <- lme(prop ~  Water_table + Days_sc, random = ~1|Epi_code, data=phy9, method="ML")

AICc(vmod1)
r.squaredGLMM(vmod1)
r.squaredGLMM(vmod2)
r.squaredGLMM(vmod3)

check_model(vmod1)

vmod2 <-  lmer(prop ~  Days_sc + Water_table + Water_table:Days_sc + (1|Epi_code), data=phy8, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = vmod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 

## WPS-2
plot_prop2 %>% filter(Phylum == "WPS-2") -> phy10
wps_mod <- lme(prop ~ Inoculation + Water_table + Days_sc + Inoculation:Water_table + Inoculation:Days_sc +
                 Water_table:Days_sc +  Days_sc:Water_table:Inoculation, 
               random = ~1|Epi_code, data=phy10, method="ML")
summary(wps_mod)
wps_dr <- dredge(wps_mod)
wps_dr2 <- subset(wps_dr, !nested(.))
summary(model.avg(wps_dr2, subset = delta <= 6))


wps_mod_bp1 <- lme(prop ~ Water_table + Days_sc, 
                   random = ~1|Epi_code, data=phy10, method="ML")
wps_mod_bp2 <- lme(prop ~ Days_sc, 
                   random = ~1|Epi_code, data=phy10, method="ML")
wps_mod_bp3 <- lme(prop ~ Water_table, 
                   random = ~1|Epi_code, data=phy10, method="ML")
wps_mod_bp4 <- lme(prop ~ 1, 
                   random = ~1|Epi_code, data=phy10, method="ML")

AICc(wps_mod_bp1)
r.squaredGLMM(wps_mod_bp1)
r.squaredGLMM(wps_mod_bp2)
r.squaredGLMM(wps_mod_bp3)
r.squaredGLMM(wps_mod_bp4)

check_model(wps_mod_bp1)
wps.mod2 <-  lmer(prop ~  Days_sc + Water_table + (1|Epi_code), data=phy10, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = wps.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput)