######## Shepherd et al., 2022 Peat moss inoculation depends on water table depth: Ecosystem function analysis
rm(list=ls())
library(tidyverse); library(ggplot2); library(vegan); library(ggfortify); library(lme4); library(DHARMa)
library(ggpubr); library(FD); library(Matching); library(plotrix); library(performance); library(MuMIn)
library(visreg)

# setwd("~/Documents/Dutch_innoculation_exp/Final_Code_data") 
setwd("C:/Users/k2256287/OneDrive - King's College London/Documents/PhD_work/Inoc_study/Repository_data/PURE_round1")


###### Starting off with carbon fluxes (NEE and CH4)

read.csv("Shepherd_et_al_flux_dataR1.csv") %>% filter(Gas=="CH4") %>% filter(Flux_corrected > -5 & Flux_corrected < 5) %>%
  dplyr::mutate(Flux_corrected = (Flux_corrected)*-1) -> ch4_data
## removed outliers in the methane data: strong weight on later models 

read.csv("Shepherd_et_al_flux_dataR1.csv") %>% filter(Gas=="CO2") %>% mutate(Flux_corrected = (Flux_corrected/1000)*-1) -> co2_data
### change the sign around on the flux data: so positive value = net uptake


####### Net ecosystem exchange ##########
co2_data$Plot <- as.factor(co2_data$Plot)

## checking for outliers:
par(mfrow=c(1,2))
boxplot((co2_data$Flux_corrected))
dotchart(co2_data$Flux_corrected)

## homogeneity in variance: non-numerical variables: these are not identical but fine to continue with (difference is less than 4)
ggplot() + geom_boxplot(aes(x= Ino, y=Flux_corrected), data = co2_data) 
ggplot() + geom_boxplot(aes(x= WL, y=Flux_corrected), data = co2_data)


### testing normality (although residuals are more important but good indicator)
hist(co2_data$Flux_corrected, probability=T, main="Density plot of flux",xlab="Flux (ug CO2eq/m2/day)", breaks = 20)
lines(density(co2_data$Flux_corrected),col=2)
ggqqplot(co2_data$Flux_corrected)
shapiro.test(co2_data$Flux_corrected) # p < 0.001
# check_distribution(co2_data$Flux_corrected) ### will keep as untransformed

co2_data$Days_sc <- scale(co2_data$Days) ## scale days

### plotting fluxes against each variable
p0 <- ggplot(aes(y=Flux_corrected, x=Ino), data=co2_data) + geom_point() + geom_smooth(method = "loess")
p1 <- ggplot(aes(y=Flux_corrected, x=WL), data=co2_data) + geom_point()
p2 <- ggplot(aes(y=Flux_corrected, x=Days_sc), data=co2_data) + geom_point()
ggarrange(p0, p1, p2, ncol=2, nrow=2)

### testing random effects structures
library(nlme)
t1 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc, random = ~1|Plot, data=co2_data)
t2 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc, random = ~1|Group/Plot, data=co2_data)
t3 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc, random = ~1|Group, data=co2_data)
t4 <- lm(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc, data=co2_data)

AIC(t1) # 2264.708
AIC(t2) # 2266.708
AIC(t3) # 2330.887
AIC(t4) # 2356.476

## stick with plot as an intercept
co2.mod <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc, random = ~1|Plot, 
               data=co2_data, method = "ML")
summary(co2.mod)
anova(co2.mod)

# dredge(co2.mod) #### full model, no need to average

# check_model(co2.mod) ### looks good
as.data.frame(round(summary(co2.mod)$tTable, 3))
AIC(co2.mod) ## 2287.796
r.squaredGLMM(co2.mod)
check_collinearity(co2.mod)
car::vif(co2.mod) ## need to examine GVIFs as there are interactions

co2.mod2 <-  lmer(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Ino:WL:Days_sc +
                    (1|Plot), data=co2_data, REML=FALSE)

simulationOutput <- simulateResiduals(fittedModel = co2.mod2, plot = F)
testDispersion(simulationOutput) ### no issue of under or over dispersion
plot(simulationOutput) ### no crossing over of quantile deviations so will stick with it: looks a bit like underdispersion but this is generally fine to continue with
outliers(simulationOutput) 

r.squaredGLMM(co2.mod) ## 25.5% of variation marginal, 47.8% conditional

par(mfrow=c(1,2))
qqnorm(residuals(co2.mod)) 
plot(resid(co2.mod), co2_data$Days)
### looks generally fine although outliers appear here as well

## examining model partial residuals: fits pretty well
par(mfrow=c(2,1))
visreg(co2.mod2, 'Days_sc', by='Ino', overlay= TRUE, jitter=TRUE, cond=list(WL="H"))
visreg(co2.mod2, 'Days_sc', by='Ino', overlay= TRUE, jitter=TRUE, cond=list(WL="L"))

### Producing time series plot 
co2_data %>% dplyr::select(Days, Plot, Ino, WL, Flux_corrected) %>%
  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("Ino", "WL", "Days"), remove=FALSE) -> co2_2

co2_2 %>% group_by(Days, Ino, WL, Treatment) %>%
  dplyr::summarise(across(c("Flux_corrected"), ~ mean(.x, na.rm = TRUE))) -> co2_means

co2_2 %>% group_by(Days, Ino, WL, Treatment) %>%
  dplyr::summarise(across(c("Flux_corrected"), ~ std.error(.x, na.rm = TRUE)))  %>%
  rename(se_co2 = Flux_corrected) %>% ungroup() %>%
  dplyr::select(Treatment, se_co2) %>% full_join(co2_means, by="Treatment") -> co2_values

co2_values$WL <- factor(co2_values$WL, levels=c("Low", "High"))

co2_2$WL <- factor(co2_2$WL, levels=c("Low", "High"))

as.data.frame(co2_2$WL) %>% mutate_if(is.factor, str_replace_all, pattern="High", replacement="High water table") %>%
  mutate_if(is.character, str_replace_all, pattern="Low", replacement="Low water table") %>%
  dplyr::rename(WTD = 'co2_2$WL')-> wtd_co2b

cbind(co2_2, wtd_co2b) -> co2_2b

as.data.frame(co2_values$WL) %>% mutate_if(is.factor, str_replace_all, pattern="High", replacement="High water table") %>%
  mutate_if(is.character, str_replace_all, pattern="Low", replacement="Low water table") %>%
  dplyr::rename(WTD = 'co2_values$WL') -> wtd_co22

cbind(co2_values, wtd_co22) %>% rename(Inoculated = "Ino") -> co2_values2

co2_values2$WTD <- factor(co2_values2$WTD, levels=c("High water table", "Low water table"))
co2_2b$WTD <- factor(co2_2b$WTD, levels=c("High water table", "Low water table"))

co2_values2$WL <- factor(co2_values2$WL, levels=c("High", "Low"))
co2_2b$WL <- factor(co2_2b$WL, levels=c("High", "Low"))


co2_values_sorted <-
  co2_values2 %>%
  arrange(Inoculated, WL, Days)

Fig4a <- ggplot() +
  geom_hline(yintercept = c(0,0), linetype="dotted") +
  geom_pointrange(aes(x= Days, y=Flux_corrected, ymax=Flux_corrected+se_co2,
                      ymin=Flux_corrected-se_co2, fill=WL, shape=Inoculated, colour=WL), position = position_jitter(width = 1.5, seed=1234), data=co2_values_sorted) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_fill_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) + 
  geom_line(aes(x=Days, y=Flux_corrected, linetype=Inoculated, colour=WL), data=co2_values_sorted, position = position_jitter(width = 1.5, seed=1234)) +
  guides(color=guide_legend("Water table"), fill = "none", linetype="none") +
  scale_linetype_manual(values = c("Yes" = "solid", "No" = "dashed")) +
  ylab(expression(CO[2]~g~m^-2~day^-1)) +
  xlab("Days since inoculation") +
  theme_bw() +
  scale_shape_manual(values=c(1, 19)) +
  theme(axis.title.x=element_blank(),
        strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_text(face="bold", size=12)) +
  ggtitle("(a) Net ecosystem exchange") +
  theme(axis.text = element_text(size=12), 
        axis.title.x=element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

#### Methane fluxes
ch4_data$Plot <- as.factor(ch4_data$Plot)

## checking for outliers
par(mfrow=c(1,2))
boxplot((ch4_data$Flux_corrected))
dotchart(ch4_data$Flux_corrected) ### dot chart doesnt look too bad in the grand sheme, but could well be underdispersion here

## homogeneity in variance: non-numerical variables: these are not identical but fine to continue with (difference is less than 4)
ggplot() + geom_boxplot(aes(x= Ino, y=Flux_corrected), data = ch4_data) 
ggplot() + geom_boxplot(aes(x= WL, y=Flux_corrected), data = ch4_data)


### testing normality (although residuals are more important)
hist(ch4_data$Flux_corrected, probability=T, main="Density plot of CH4 flux",xlab="Flux (ug CH4/m2/day)", breaks = 20)
lines(density(ch4_data$Flux_corrected),col=2)
ggqqplot(ch4_data$Flux_corrected)
shapiro.test(ch4_data$Flux_corrected) # p < 0.001
# check_distribution(ch4_data$Flux_corrected)
ch4_data$Days_sc <- scale(ch4_data$Days) ## scale days

### plotting fluxes against each variable
p0 <- ggplot(aes(y=Flux_corrected, x=Ino), data=ch4_data) + geom_point() + geom_smooth(method = "loess")
p1 <- ggplot(aes(y=Flux_corrected, x=WL), data=ch4_data) + geom_point()
p2 <- ggplot(aes(y=Flux_corrected, x=Days_sc), data=ch4_data) + geom_point()
ggarrange(p0, p1, p2, ncol=2, nrow=2)

### testing random effects structures
ch4_data$Group <- as.factor(ch4_data$Group)
t1 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:Ino:WL, random = ~1|Plot, data=ch4_data, method = "ML")
t2 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:Ino:WL, random = ~1|Group/Plot, data=ch4_data, method = "ML")
t3 <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:Ino:WL, random = ~1|Group, data=ch4_data, method = "ML")
t4 <- lm(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:Ino:WL, data=ch4_data,)

AIC(t1) # 1375.933
AIC(t2) # 1377.933
AIC(t3) # 1460.9
AIC(t4) # 1467.771

## Plot as a random intercept
CH4.mod <- lme(Flux_corrected ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:WL:Ino, random = ~1|Plot, data=ch4_data,
               method = "ML")
summary(CH4.mod)
anova(CH4.mod)

ch4_dr <- dredge(CH4.mod, rank = "AIC")
ch4_dr2 <- subset(ch4_dr, !nested(.))


CH4.mod <- lme(Flux_corrected ~ Days_sc + WL + Days_sc:WL, random = ~1|Plot, data=ch4_data,
               method = "ML")
as.data.frame(round(summary(CH4.mod)$tTable, 3))
AIC(CH4.mod) ## 763.1
r.squaredGLMM(CH4.mod)
# 21.8%, 48.8% 


# check_model(CH4.mod) ## 

CH4.mod2 <-  lmer(Flux_corrected ~ WL + Days_sc + WL:Days_sc +
                    (1|Plot), data=ch4_data, REML=FALSE)
simulationOutput <- simulateResiduals(fittedModel = CH4.mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) 
outliers(simulationOutput) 
### generally this is fine to continue, no crossing over of quantules and some underdispersion but fits ok, partial deviation from residual normality also fine
### as models tend to be able to deal with underdispersion qiute well (and are conservative in estimates) I will stick with the model as written

par(mfrow=c(1,2))
qqnorm(residuals(CH4.mod)) 
plot(resid(CH4.mod), ch4_data$Days)


ch4_data %>% dplyr::select(Days, Plot, Ino, WL, Flux_corrected) %>%
  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("Ino", "WL", "Days"), remove=FALSE) -> ch4_2

ch4_2 %>% group_by(Days, Ino, WL, Treatment) %>%
  summarise(across(c("Flux_corrected"), ~ mean(.x, na.rm = TRUE))) -> ch4_means

ch4_2 %>% group_by(Days, Ino, WL, Treatment) %>%
  summarise(across(c("Flux_corrected"), ~ std.error(.x, na.rm = TRUE)))  %>%
  rename(se_ch4 = "Flux_corrected") %>% ungroup() %>%
  dplyr::select(Treatment, se_ch4) %>% full_join(ch4_means, by="Treatment") %>% rename(Inoculated = "Ino") -> ch4_values

ch4_values$WL <- factor(ch4_values$WL, levels=c("High", "Low"))
ch4_2$WL <- factor(ch4_2$WL, levels=c("High", "Low"))

ch4_values_sorted <-
  ch4_values %>%
  arrange(Inoculated, WL, Days)

Fig4b <- ggplot() +
  geom_hline(yintercept = c(0,0), linetype="dotted") +
  geom_pointrange(aes(x= Days, y=Flux_corrected, ymax=Flux_corrected+se_ch4,
                      ymin=Flux_corrected-se_ch4, fill=WL, shape=Inoculated, colour=WL), position = position_jitter(width = 2.5, seed=1234), data=ch4_values_sorted) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  guides(color=guide_legend("Water table"), fill = "none", linetype="none") +
  scale_linetype_manual(values = c("Yes" = "solid", "No" = "dashed")) +
  geom_line(aes(x=Days, y=Flux_corrected, linetype=Inoculated, colour=WL), data=ch4_values_sorted, position = position_jitter(width = 2.5, seed=1234)) +
  ylab(expression(CH[4]~mg~m^-2~day^-1)) +
  xlab("Days since inoculation") +
  theme_bw() +
  scale_shape_manual(values=c(1, 19)) +
  # facet_wrap(~WL) +
  theme(axis.title.x=element_blank(),
        strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_blank()) +
  ggtitle("(b) Methane exchange") +
  theme(axis.text = element_text(size=12), 
        axis.title.x=element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

######### Pore water composition

pw <- read.csv("Shepherd_Inoc_PWcomp.csv", row.names = 1)
pw_pca <- prcomp(pw[,c(13:19)], center=TRUE, scale=TRUE)
summary(pw_pca)
as.data.frame(pw_pca$x) -> pc_coords
pw %>% dplyr::select(Treatment:Rep) -> pw_trt

cbind(pw_trt, pc_coords) %>% dplyr::select(Treatment:PC2) -> pw_pca_res

## checking for outliers
par(mfrow=c(1,2))
boxplot(pw_pca_res$PC1)
dotchart(pw_pca_res$PC1)

## homogeneity in variance: non-numerical variables: these are not identical but fine to continue with (difference is less than 4)
ggplot() + geom_boxplot(aes(x= Ino, y=PC1), data = pw_pca_res) 
ggplot() + geom_boxplot(aes(x= WL, y=PC1), data = pw_pca_res)


### testing normality (although residuals are more important)
hist(pw_pca_res$PC1, probability=T, main="Density plot of flux",xlab="Flux (ug CO2eq/m2/day)", breaks = 20)
lines(density(pw_pca_res$PC1),col=2)
ggqqplot(pw_pca_res$PC1)
shapiro.test(pw_pca_res$PC1) # p = 0.00958
# check_distribution(pw_pca_res$PC1) ### will keep as untransformed

pw_pca_res$Days_sc <- scale(pw_pca_res$Date) ## scale days

### plotting fluxes against each variable
p0 <- ggplot(aes(y=PC1, x=Ino), data=pw_pca_res) + geom_point() + geom_smooth(method = "loess")
p1 <- ggplot(aes(y=PC1, x=WL), data=pw_pca_res) + geom_point()
p2 <- ggplot(aes(y=PC1, x=Days_sc), data=pw_pca_res) + geom_point()
ggarrange(p0, p1, p2, ncol=2, nrow=2)

### testing random effects structures
pw_pca_res$meso <- as.factor(pw_pca_res$meso)
pw_pca_res$Rep <- as.factor(pw_pca_res$Rep)
t1 <- lme(PC1 ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:WL:Ino, random = ~1|meso, data=pw_pca_res, method="ML")
t2 <- lme(PC1 ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:WL:Ino, random = ~1|Rep/meso, data=pw_pca_res, method="ML")
t3 <- lme(PC1 ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:WL:Ino, random = ~1|Rep, data=pw_pca_res, method="ML")
t4 <- lm(PC1 ~ Days_sc + Ino + WL + Days_sc:Ino + Ino:WL + Days_sc:WL + Days_sc:WL:Ino, data=pw_pca_res)

AICc(t1) # 235.45
AICc(t2) # 238.21
AICc(t3) # 236.19
AICc(t4) # 233.53

## Stick with mesocosm as a random intercept: even though no RE performs better: stay true to structure
pw_pca_mod <- lme(PC1 ~ Ino + WL + Days_sc + Ino:WL + Ino:Days_sc + WL:Days_sc +  Days_sc:WL:Ino, 
                  random = ~1|meso, data=pw_pca_res, method="ML")
summary(pw_pca_mod)
as.data.frame(round(summary(pw_pca_mod)$tTable, 3))
anova(pw_pca_mod)

# check_model(pw_pca_mod) 

car::vif(pw_pca_mod) ## need to examine GVIFs as there are interactions


pw_pca_mod2 <-  lmer(PC1 ~ Ino + WL + Days_sc + Ino:WL + Ino:Days_sc + WL:Days_sc + Days_sc:Ino:WL +
                       (1|meso), data=pw_pca_res, REML=FALSE)

simulationOutput <- simulateResiduals(fittedModel = pw_pca_mod2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) ### 
outliers(simulationOutput) ##

r.squaredGLMM(pw_pca_mod) ## 73.6% marginal, 75.9% conditional

pw_pca_dr <- dredge(pw_pca_mod, rank="AIC") ### best performing model has inoculation and date in it, and their interaction
subset(pw_pca_dr, !nested(.)) ### just the single model then, so continue with that


pca_bpm <- lme(PC1 ~ Ino + Days_sc + Ino:Days_sc, 
               random = ~1|meso, data=pw_pca_res, method="ML")
as.data.frame(round(summary(pca_bpm)$tTable, 3))
AICc(pca_bpm)
r.squaredGLMM(pca_bpm)
# check_model(pca_bpm)
pca.bpm2 <-  lmer(PC1 ~ Ino + Days_sc + Ino:Days_sc +
                    (1|meso), data=pw_pca_res, REML=FALSE)

simulationOutput <- simulateResiduals(fittedModel = pca.bpm2, plot = F)
testDispersion(simulationOutput) 
plot(simulationOutput) #
outliers(simulationOutput)

par(mfrow=c(1,1))
visreg(pca.bpm2, "Days_sc", by="Ino", overlay=TRUE, jitter= TRUE, ylab="Pore water PCA")

pw_pca_res %>%
  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(Treatment2, c("Treatment", "Date"), remove = FALSE) -> pw_pca_res2

pw_pca_res2 %>% group_by(Date, Ino, WL, Treatment2) %>%
  dplyr::summarise(across(c("PC1"), ~ mean(.x, na.rm = TRUE))) -> pw_means

pw_pca_res2 %>% group_by(Date, Ino, WL, Treatment2) %>%
  dplyr::summarise(across(c("PC1"), ~ std.error(.x, na.rm = TRUE))) %>%
  rename(se_pw = PC1) %>% ungroup() %>%
  dplyr::select(Treatment2, se_pw) %>% full_join(pw_means, by="Treatment2") -> pw_values

pw_values_sorted <-
  pw_values %>%
  arrange(Ino, WL, Date) %>% rename(Inoculated = Ino)

Fig4c <- ggplot() +
  geom_hline(yintercept = c(0,0), linetype="dotted") +
  geom_pointrange(aes(x=Date, y=PC1, ymax=PC1+se_pw,
                      ymin=PC1-se_pw, fill=WL, shape=Inoculated, colour=WL), position = position_jitter(width=2, seed=1234), data=pw_values_sorted) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  geom_line(aes(x=Date, y=PC1, linetype=Inoculated, colour=WL, fill=WL), data=pw_values_sorted, position = position_jitter(width=2, seed=1234)) +
  guides(color=guide_legend("Water table"), fill = "none", linetype="none") +
  scale_linetype_manual(values = c("Yes" = "solid", "No" = "dashed")) +
  ylab("PCA1") +
  xlab("Days since inoculation") +
  theme_bw() +
  scale_shape_manual(values=c(1, 19)) +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_blank()) +
  ggtitle("(c) Porewater composition") +
  theme(axis.text = element_text(size=12), 
        axis.title.x=element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  xlim(0, 100)

####### Soil CN ratio

read.csv("Shepherd_Inoc_CN_R1.csv") %>% filter(SoilC != "NA") -> cn

par(mfrow=c(1,2))
boxplot((cn$SoilCN))
dotchart(cn$SoilCN) ### seems fine

## homogeneity in variance: non-numerical variables: these are not identical but fine to continue with (difference is less than 4)
ggplot() + geom_boxplot(aes(x= Ino, y=SoilCN), data = cn) 
ggplot() + geom_boxplot(aes(x= WL, y=SoilCN), data = cn)


### testing normality (although residuals are more important)
hist(log(cn$SoilCN), probability=T, main="",xlab="CN ratio", breaks = 20)
lines(density(log(cn$SoilCN)),col=2)
ggqqplot(log(cn$SoilCN))
shapiro.test(log(cn$SoilCN)) 
# check_distribution(log(cn$SoilCN)) ## gonna log transform


### checking colinearity of numerical parameters
num_eff <- dplyr::select(cn, Ino, WL) 
num_eff$Ino <- as.numeric(num_eff$Ino)
num_eff$WL <- as.numeric(num_eff$WL)
cor(num_eff) ### all fine

### plotting fluxes against each variable
p0 <- ggplot(aes(y=SoilCN, x=Ino), data=cn) + geom_point() + geom_smooth(method = "loess")
p1 <- ggplot(aes(y=SoilCN, x=WL), data=cn) + geom_point()
ggarrange(p0, p1, ncol=2, nrow=1)


cnmod <- lm(log(SoilCN) ~ Ino + WL + Ino:WL, data=cn)
summary(cnmod)
anova(cnmod)
cn_aov <- aov(log(SoilCN) ~ Ino + WL + Ino:WL, data=cn)
summary(cn_aov)
# leveneTest(residuals(cnmod)~ Ino*WL, data=cn) 
ks.test(residuals(cnmod), "pnorm")

cn_aov <- aov(log(SoilCN) ~ Ino + WL + Ino:WL, data=cn)
qqnorm(residuals(cn_aov))
# check_model(cnmod) ### will stick with as is

cn %>% mutate(SoilCN = log(SoilCN)) %>% 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") -> cn2

cn2 %>% group_by(Ino, WL, Treatment) %>% dplyr::select(Code:VegCN) %>%
  dplyr::summarise(across(c("SoilC":"VegCN"), ~ mean(.x, na.rm = TRUE))) -> mn_cn

cn2 %>% group_by(Ino, WL, Treatment) %>% dplyr::select(Code:VegCN) %>%
  dplyr::summarise(across(c("SoilC":"VegCN"), ~ std.error(.x, na.rm = TRUE))) %>%
  rename(se_soilcn = SoilCN) %>% ungroup() %>%
  dplyr::select(Treatment, se_soilcn) %>% full_join(mn_cn, by="Treatment") -> cn_vals


Fig4d <- ggplot(aes(x=Ino, y=SoilCN, fill=WL, colour=WL, shape=Ino), data=cn_vals) +
  geom_pointrange(aes(ymax=SoilCN+se_soilcn,
                      ymin=SoilCN-se_soilcn), position = position_dodge(0.4)) +
  geom_jitter(aes(x=Ino, y=SoilCN, fill=WL), data=cn2,
              alpha=0.6, position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.4),
              size=1.4) +
  guides(color=guide_legend("Water table"), fill = "none") +
  scale_shape_manual(values=c(1, 19)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  ylab("C:N (log)") +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_text(face="bold", size=12)) +
  xlab("Inoculated") +
  ggtitle("(d) Peat C:N ratio") +
  theme(axis.text = element_text(size=12), 
        axis.title.x=element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  scale_y_continuous(breaks = c(3.2, 3.6, 4))

###### Dry weight aboveground biomass

read.csv("Shepherd_Inoc_VegDW.csv") %>% mutate(NetDW = NetDW*0.16)-> dwveg
### converts to M squared for easy comparison
check_distribution((dwveg$NetDW))
hist(dwveg$NetDW)
shapiro.test(dwveg$NetDW) ## fine

dwmod <- lm(NetDW ~ Ino*WL, data=dwveg)
summary(dwmod)
anova(dwmod)
summary(aov((NetDW) ~ Ino*WL, data=dwveg)) ### use result here
# check_model(dwmod)

options(na.action = na.fail)
dredge(dwmod, rank="AICc") ## water table in best performing model

plot(resid(dwmod), (dwveg$NetDW))
qqnorm(residuals(dwmod)) ## seems alright
ks.test(residuals(dwmod), "pnorm") ## violates however wuth small sample size this is difficult to pass

# leveneTest(residuals(dwmod)~ Ino*WL, data=dwveg) # all fine

dwveg %>% 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") -> dw2

dw2 %>% group_by(Ino, WL, Treatment) %>%
  dplyr::summarise(across(c("DW":"NetDW"), ~ mean(.x, na.rm = TRUE))) -> mn_dw

dw2 %>% group_by(Ino, WL, Treatment) %>%
  dplyr::summarise(across(c("DW":"NetDW"), ~ std.error(.x, na.rm = TRUE))) %>%
  rename(se_dwnet = NetDW) %>% ungroup() %>%
  dplyr::select(Treatment, se_dwnet) %>% full_join(mn_dw, by="Treatment") -> dw_vals


Fig4e <- ggplot(aes(x=Ino, y=NetDW, fill=WL, colour=WL, shape=Ino), data=dw_vals) +
  geom_pointrange(aes(ymax=NetDW+se_dwnet,
                      ymin=NetDW-se_dwnet), position = position_dodge(0.4)) +
  geom_jitter(aes(x=Ino, y=NetDW, fill=WL, shape = Ino), data=dw2,
              alpha=0.6, position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.4), size=1.4) +
  guides(color=guide_legend("Water table"), fill = "none") +
  ylab(expression(Dry~weight~(g~m^-2))) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  scale_shape_manual(values=c(1, 19)) +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text = element_text(face="bold", size=12)) +
  xlab("Inoculated") +
  ggtitle("(e) Aboveground biomass") +
  theme(axis.text = element_text(size=12), 
        axis.title.x=element_text(size=14),
        axis.title.y = element_text(size=14),
        plot.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14)) +
  ylim(0,11)

#### Figure 4: compile from
tseries_gg  <- ggarrange(Fig4a,
                         Fig4b,
                         nrow=2, common.legend = TRUE, legend = "right", align = "hv")
ggarrange(Fig4c, Fig4d, Fig4e, ncol=3, legend = "none")

fig_top <- tseries_gg
fig_bot <- ggarrange(Fig4c, Fig4d, Fig4e, ncol=3, nrow=1, common.legend = TRUE, legend=FALSE, align = "hv") 


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_top, vp = define_region(row = 1:65, col = 1:100))   # Span over two columns
print(fig_bot, vp = define_region(row = 68:100, col = 1:91))
# export 1200 x 750


###### Supplementary figure: plotting individual C and N values

### Plotting relationship between CN and N

SCplot <- ggplot() +
  geom_point(aes(x=SoilCN, y=SoilC,  fill=WL, colour=WL, shape=Ino), data=cn2) +
  guides(color=guide_legend("Water table"), fill = "none", shape = guide_legend("Inoculated")) +
  scale_shape_manual(values=c(1, 19)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  theme_bw() +
  ylab("Peat C (%)") +
  xlab("Peat C:N")

SNplot <- ggplot() +
  geom_point(aes(x=SoilCN, y=SoilN,  fill=WL, colour=WL, shape=Ino), data=cn2) +
  guides(color=guide_legend("Water table"), fill = "none", shape = guide_legend("Inoculated")) +
  scale_shape_manual(values=c(1, 19)) +
  scale_colour_manual(values=c("High" = "#0072B2", "Low" = "#D55E00")) +
  theme_bw() +
  stat_smooth(aes(x=SoilCN, y=SoilN), data=cn2, method="lm", col="black") +
  ylab("Peat N (%)") +
  xlab("Peat C:N")

ggarrange(SCplot, SNplot, ncol=2, nrow=1, common.legend = TRUE, legend="right") 

##### Porewater PCA

pw_trt %>% 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::rename(Inoculated = Ino) %>% dplyr::rename(Day = Date) %>% unite(Trt, c("Inoculated", "WL")) -> nut_trt2

nut_trt2$Day <- as.factor(nut_trt2$Day)

autoplot(pw_pca, data=nut_trt2, colour = "Trt", shape = "Day", loadings=TRUE, loadings.colour = 'gray4',
         loadings.label = TRUE, loadings.label.size = 3, label.size=5, loadings.label.vjust = 1.3, loadings.label.hjust = 0.8) +  
  guides(color=guide_legend("Inoculated_Water table")) + 
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "inside",
        strip.text=element_text(face="bold", size=12)) +
  ylab("PC2 (15%)") +
  xlab("PC1 (54%)")



pw %>% dplyr::group_by(Ino, WL, Treatment, Date) %>%
  dplyr::summarise(across(c("pH":"Cl"), ~ mean(.x, na.rm = TRUE))) %>%
  tidyr::gather(nut, conc, 5:11) %>%
  unite(Tr2, c(Treatment, Date, nut), remove = FALSE) -> nut_mns

pw %>%dplyr::group_by(Ino, WL, Treatment, Date) %>%
  dplyr::summarise(across(c("pH":"Cl"), ~ std.error(.x, na.rm = TRUE))) %>%
  ungroup() %>%
  tidyr::gather(nut, se, 5:11) %>%
  unite(Tr2, c(Treatment, Date, nut)) %>% 
  dplyr::select(Tr2, se) %>%
  full_join(nut_mns, by="Tr2") -> nut_vals

ggplot() +
  geom_pointrange(aes(x=Date, y=conc, ymax=conc+se,
                      ymin=conc-se, fill=WL, shape=Ino, colour=WL), position = position_dodge(3), data=nut_vals) +
  scale_colour_manual(values=c("H" = "#0072B2", "L" = "#D55E00")) +
  scale_linetype_manual(values = c("N" = "dashed", "Y" = "solid")) +
  guides(color=guide_legend("Water table"), fill = "none", shape=guide_legend("Inoculated"), linetype =guide_legend("Inoculated")) +
  geom_line(aes(x=Date, y=conc, linetype=Ino, colour=WL), data=nut_vals, position = position_dodge(3)) +
  xlab("Days since inoculation") +
  ylab("") +
  theme_bw() +
  scale_shape_manual(values=c(1, 19)) +
  facet_wrap(~nut, scales = "free_y", nrow=3, ncol=3) +
  theme(strip.background = element_blank(),
        strip.placement = "inside")