# Carbon content and other soil properties of near-surface peats before and after peatland restoration
# Statistics and figures
#
# last update 11.05.2023
#
# Robert Probst (HNE-Eberswalde)

#packages
install.packages("EnvStats")
install.packages("psych")
install.packages("ggpubr")
install.packages("car")
install.packages("ggplot2")
install.packages("gridExtra")
library(psych)
library(car)
library(ggplot2)
library(gridExtra)
library(EnvStats)
library(ggpubr)

#working directory and dataset / Change Working directory to your path
setwd("C:/users/robert/desktop/stechlin/stechlin/Final")
Corg<- read.table("Statistics_R_Stechlin_final.csv", header=TRUE, dec=",", sep=";", encoding="UTF-8")

#descriptive statistics ("psych")
describeBy(Corg$TOC, group = Corg$site_number_sampling)
describeBy(Corg$TN, group = Corg$site_number_sampling)
describeBy(Corg$CN, group = Corg$site_number_sampling)
describeBy(Corg$pH, group = Corg$site_number_sampling)
describeBy(Corg$DBD, group = Corg$site_number_sampling)

#boxplot TOC
plot_Corg<- ggplot(Corg, aes(x = sampling_1, y = TOC, group=site_number_sampling)) +
           stat_boxplot(geom = "errorbar", width = 0.25)+
           geom_boxplot()+
           xlab("sampling years")+
           ylab("TOC [%]")+
           scale_x_discrete()+
           scale_y_continuous(limits=c(20,60),breaks=seq(20,60,10))+
           geom_jitter(size=2,alpha = 0.6, width = 0.2, height = 0, col="grey39")+ 
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=8),
           axis.text.y=element_text(size=11),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           stat_n_text(y.pos = 20, size = 4)+
           facet_grid(cols = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_Corg, width=27.94, height=21.59, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "TOC-Boxplot-final.pdf")

#boxplot TN
plot_N<- ggplot(Corg, aes(x = sampling_1, y = TN, group=site_number_sampling)) +
           stat_boxplot(geom = "errorbar", width = 0.25)+
           geom_boxplot()+
           xlab("sampling years")+
           ylab("TN [%]")+
           scale_x_discrete()+
           scale_y_continuous(breaks=seq(0,6,0.5))+
           geom_jitter(size=2,alpha = 0.6, width = 0.2, height = 0, col="grey39")+ 
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=8),
           axis.text.y=element_text(size=11),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           stat_n_text(y.pos = 0, size = 4)+
           facet_grid(cols = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_N, width=27.94, height=21.59, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "TN-Boxplot-final.pdf")

#boxplot CN
plot_CN<- ggplot(Corg, aes(x = sampling_1, y = CN, group=site_number_sampling)) +
           stat_boxplot(geom = "errorbar", width = 0.25)+
           geom_boxplot()+
           xlab("sampling years")+
           ylab("C/N ratio")+
           scale_x_discrete()+
           scale_y_continuous(breaks=seq(0,70,10))+
           geom_jitter(size=2,alpha = 0.6, width = 0.2, height = 0, col="grey39")+ 
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=8),
           axis.text.y=element_text(size=11),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           stat_n_text(y.pos = 0, size = 4)+
           facet_grid(cols = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_CN, width=27.94, height=21.59, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "CN-Boxplot-final.pdf")

#boxplot DBD
plot_TRD<- ggplot(Corg, aes(x = sampling_1, y = DBD, group=site_number_sampling)) +
           stat_boxplot(geom = "errorbar", width = 0.25)+
           geom_boxplot()+
           xlab("sampling years")+
           ylab("dry bulk density [g/cm³]")+
           scale_x_discrete()+
           scale_y_continuous(breaks=seq(0.0,0.7,0.1))+
           geom_jitter(size=2,alpha = 0.6, width = 0.2, height = 0, col="grey39")+ 
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=8),
           axis.text.y=element_text(size=11),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           stat_n_text(y.pos = 0, size = 4)+
           facet_grid(cols = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_TRD, width=27.94, height=21.59, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "DBD-Boxplot-final.pdf")

#boxplot pH
plot_pH<- ggplot(Corg, aes(x = sampling_1, y = pH, group=site_number_sampling)) +
           stat_boxplot(geom = "errorbar", width = 0.25)+
           geom_boxplot()+
           xlab("sampling years")+
           ylab("PH-value")+
           scale_x_discrete()+
           scale_y_continuous(breaks=seq(2,8,1))+
           geom_jitter(size=2,alpha = 0.6, width = 0.2, height = 0, col="grey39")+ 
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=8),
           axis.text.y=element_text(size=11),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           stat_n_text(y.pos = 2, size = 4)+
           facet_grid(cols = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_pH, width=27.94, height=21.59, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "pH-Boxplot-final.pdf")


#Test on normal distribution each sampling date (Shapiro-Wilk)
Corg_11 <- subset(Corg,site_number_sampling=="11")
Corg_12 <- subset(Corg,site_number_sampling=="12")
Corg_21 <- subset(Corg,site_number_sampling=="21")
Corg_22 <- subset(Corg,site_number_sampling=="22")
Corg_31 <- subset(Corg,site_number_sampling=="31")
Corg_32 <- subset(Corg,site_number_sampling=="32")
Corg_41 <- subset(Corg,site_number_sampling=="41")
Corg_42 <- subset(Corg,site_number_sampling=="42")
Corg_51 <- subset(Corg,site_number_sampling=="51")
Corg_52 <- subset(Corg,site_number_sampling=="52")
Corg_61 <- subset(Corg,site_number_sampling=="61")
Corg_62 <- subset(Corg,site_number_sampling=="62")

shapiro.test(Corg_11$TOC)
shapiro.test(Corg_12$TOC)
shapiro.test(Corg_21$TOC)
shapiro.test(Corg_22$TOC)
shapiro.test(Corg_31$TOC)
shapiro.test(Corg_32$TOC)
shapiro.test(Corg_41$TOC)
shapiro.test(Corg_42$TOC)
shapiro.test(Corg_51$TOC)
shapiro.test(Corg_52$TOC)
shapiro.test(Corg_61$TOC)
shapiro.test(Corg_62$TOC)

shapiro.test(Corg_11$TN)
shapiro.test(Corg_12$TN)
shapiro.test(Corg_21$TN)
shapiro.test(Corg_22$TN)
shapiro.test(Corg_31$TN)
shapiro.test(Corg_32$TN)
shapiro.test(Corg_41$TN)
shapiro.test(Corg_42$TN)
shapiro.test(Corg_51$TN)
shapiro.test(Corg_52$TN)
shapiro.test(Corg_61$TN)
shapiro.test(Corg_62$TN)

shapiro.test(Corg_11$CN)
shapiro.test(Corg_12$CN)
shapiro.test(Corg_21$CN)
shapiro.test(Corg_22$CN)
shapiro.test(Corg_31$CN)
shapiro.test(Corg_32$CN)
shapiro.test(Corg_41$CN)
shapiro.test(Corg_42$CN)
shapiro.test(Corg_51$CN)
shapiro.test(Corg_52$CN)
shapiro.test(Corg_61$CN)
shapiro.test(Corg_62$CN)

shapiro.test(Corg_11$pH)
shapiro.test(Corg_12$pH)
shapiro.test(Corg_21$pH)
shapiro.test(Corg_22$pH)
shapiro.test(Corg_31$pH)
shapiro.test(Corg_32$pH)
shapiro.test(Corg_41$pH)
shapiro.test(Corg_42$pH)
shapiro.test(Corg_51$pH)
shapiro.test(Corg_52$pH)
shapiro.test(Corg_61$pH)
shapiro.test(Corg_62$pH)

shapiro.test(Corg_11$DBD)
shapiro.test(Corg_12$DBD)
shapiro.test(Corg_21$DBD)
shapiro.test(Corg_22$DBD)
shapiro.test(Corg_31$DBD)
shapiro.test(Corg_32$DBD)
shapiro.test(Corg_41$DBD)
shapiro.test(Corg_42$DBD)
shapiro.test(Corg_51$DBD)
shapiro.test(Corg_52$DBD)
shapiro.test(Corg_61$DBD)
shapiro.test(Corg_62$DBD)

#Levene Test on variance homogenity
Corg$site_number_sampling_1 <-as.factor(Corg$site_number_sampling)
Corg_1 <- subset(Corg,site_number=="1")
Corg_2 <- subset(Corg,site_number=="2")
Corg_3 <- subset(Corg,site_number=="3")
Corg_4 <- subset(Corg,site_number=="4")
Corg_5 <- subset(Corg,site_number=="5")
Corg_6 <- subset(Corg,site_number=="6")
leveneTest(Corg_1$TOC, Corg_1$site_number_sampling_1)
leveneTest(Corg_2$TOC, Corg_2$site_number_sampling_1)
leveneTest(Corg_3$TOC, Corg_3$site_number_sampling_1)
leveneTest(Corg_4$TOC, Corg_4$site_number_sampling_1)
leveneTest(Corg_5$TOC, Corg_5$site_number_sampling_1)
leveneTest(Corg_6$TOC, Corg_6$site_number_sampling_1)

leveneTest(Corg_1$TN, Corg_1$site_number_sampling_1)
leveneTest(Corg_2$TN, Corg_2$site_number_sampling_1)
leveneTest(Corg_3$TN, Corg_3$site_number_sampling_1)
leveneTest(Corg_4$TN, Corg_4$site_number_sampling_1)
leveneTest(Corg_5$TN, Corg_5$site_number_sampling_1)
leveneTest(Corg_6$TN, Corg_6$site_number_sampling_1)

leveneTest(Corg_1$CN, Corg_1$site_number_sampling_1)
leveneTest(Corg_2$CN, Corg_2$site_number_sampling_1)
leveneTest(Corg_3$CN, Corg_3$site_number_sampling_1)
leveneTest(Corg_4$CN, Corg_4$site_number_sampling_1)
leveneTest(Corg_5$CN, Corg_5$site_number_sampling_1)
leveneTest(Corg_6$CN, Corg_6$site_number_sampling_1)

leveneTest(Corg_1$pH, Corg_1$site_number_sampling_1)
leveneTest(Corg_2$pH, Corg_2$site_number_sampling_1)
leveneTest(Corg_3$pH, Corg_3$site_number_sampling_1)
leveneTest(Corg_4$pH, Corg_4$site_number_sampling_1)
leveneTest(Corg_5$pH, Corg_5$site_number_sampling_1)
leveneTest(Corg_6$pH, Corg_6$site_number_sampling_1)

leveneTest(Corg_1$DBD, Corg_1$site_number_sampling_1)
leveneTest(Corg_2$DBD, Corg_2$site_number_sampling_1)
leveneTest(Corg_3$DBD, Corg_3$site_number_sampling_1)
leveneTest(Corg_4$DBD, Corg_4$site_number_sampling_1)
leveneTest(Corg_5$DBD, Corg_5$site_number_sampling_1)
leveneTest(Corg_6$DBD, Corg_6$site_number_sampling_1)

#test normal distribution on differences between sampling dates
differences <- read.table("statistics-differences-sampling.csv", header=TRUE, dec=",", sep=";")
diff_1 <- subset(differences,site=="1")
diff_2 <- subset(differences,site=="2")
diff_3 <- subset(differences,site=="3")
diff_4 <- subset(differences,site=="4")
diff_5 <- subset(differences,site=="5")
diff_6 <- subset(differences,site=="6")

shapiro.test(diff_1$TOC)
shapiro.test(diff_2$TOC)
shapiro.test(diff_3$TOC)
shapiro.test(diff_4$TOC)
shapiro.test(diff_5$TOC)
shapiro.test(diff_6$TOC)

shapiro.test(diff_1$TN)
shapiro.test(diff_2$TN)
shapiro.test(diff_3$TN)
shapiro.test(diff_4$TN)
shapiro.test(diff_5$TN)
shapiro.test(diff_6$TN)

shapiro.test(diff_1$CN)
shapiro.test(diff_2$CN)
shapiro.test(diff_3$CN)
shapiro.test(diff_4$CN)
shapiro.test(diff_5$CN)
shapiro.test(diff_6$CN)

shapiro.test(diff_1$pH)
shapiro.test(diff_2$pH)
shapiro.test(diff_3$pH)
shapiro.test(diff_4$pH)
shapiro.test(diff_5$pH)
shapiro.test(diff_6$pH)

shapiro.test(diff_1$DBD)
shapiro.test(diff_2$DBD)
shapiro.test(diff_3$DBD)
shapiro.test(diff_4$DBD)
shapiro.test(diff_5$DBD)
shapiro.test(diff_6$DBD)

#Wilcoxon-Test (preffered because of small number of samples n=10)
wilcox.test(Corg_1$TOC~Corg_1$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_2$TOC~Corg_2$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_3$TOC~Corg_3$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_4$TOC~Corg_4$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_5$TOC~Corg_5$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_6$TOC~Corg_6$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_1$TN~Corg_1$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_2$TN~Corg_2$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_3$TN~Corg_3$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_4$TN~Corg_4$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_5$TN~Corg_5$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_6$TN~Corg_6$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_1$CN~Corg_1$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_2$CN~Corg_2$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_3$CN~Corg_3$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_4$CN~Corg_4$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_5$CN~Corg_5$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_6$CN~Corg_6$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_1$pH~Corg_1$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_2$pH~Corg_2$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_3$pH~Corg_3$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_4$pH~Corg_4$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_5$pH~Corg_5$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_6$pH~Corg_6$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_1$DBD~Corg_1$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_2$DBD~Corg_2$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_3$DBD~Corg_3$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_4$DBD~Corg_4$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_5$DBD~Corg_5$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_6$DBD~Corg_6$site_number_sampling, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)


#t-test for paired samples (optional not used because of small samples n=10)
t.test(Corg_1$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_2$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_3$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_4$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_5$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_6$TOC~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")

t.test(Corg_1$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_2$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_3$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_4$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_5$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_6$TN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")

t.test(Corg_1$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_2$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_3$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_4$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_5$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_6$CN~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")

t.test(Corg_1$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_2$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_3$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_4$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_5$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_6$pH~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")

t.test(Corg_1$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_2$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_3$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_4$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_5$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")
t.test(Corg_6$DBD~Corg_1$site_number_sampling, paired = TRUE, alternative = "two.sided")


#statistics new peat
Corg_neu<- read.table("peat_old_new.csv", header=TRUE, dec=",", sep=";", encoding="UTF-8")

describeBy(Corg_neu$TOC, group = Corg_neu$site_number_old_new)
describeBy(Corg_neu$TN, group = Corg_neu$site_number_old_new)
describeBy(Corg_neu$CN, group = Corg_neu$site_number_old_new)
describeBy(Corg_neu$pH, group = Corg_neu$site_number_old_new)
describeBy(Corg_neu$DBD, group = Corg_neu$site_number_old_new)
describeBy(Corg_neu$thickness_new, group = Corg_neu$site_number_old_new)

#tests new peat only wilcoxon because few number of samplings (e.g. Seggenkuten only n=4)
Corg_neu_4 <- subset(Corg_neu,site_number=="4")
Corg_neu_5 <- subset(Corg_neu,site_number=="5")
Corg_neu_3 <- subset(Corg_neu,site_number=="3")

wilcox.test(Corg_neu_4$TOC~Corg_neu_4$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_5$TOC~Corg_neu_5$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_3$TOC~Corg_neu_3$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_neu_4$TN~Corg_neu_4$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_5$TN~Corg_neu_5$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_3$TN~Corg_neu_3$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_neu_4$CN~Corg_neu_4$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_5$CN~Corg_neu_5$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_3$CN~Corg_neu_3$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_neu_4$pH~Corg_neu_4$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_5$pH~Corg_neu_5$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_3$pH~Corg_neu_3$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)

wilcox.test(Corg_neu_4$DBD~Corg_neu_4$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_5$DBD~Corg_neu_5$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)
wilcox.test(Corg_neu_3$DBD~Corg_neu_3$site_number_old_new, alternative = "two.sided", paired=TRUE, exact = TRUE, p.adjust="bonferroni", conf.level = 0.95)


#water level figure 
Pegel      <- read.table("Statistics_waterlevel.csv", header=TRUE, dec=",", sep=";",encoding = "UTF-8")
Pegel$Month<- as.Date(Pegel$Month, format="%d.%m.%Y")
moore <- c("Beerenwiese","Boberowseewiese","Seggenkuten","Teufelseemoor","Müllerwiese","Jägerwiese")
names(moore) <- c("1","2","3","4","5","6")

plot_Pegel <- ggplot(Pegel, aes(x = Month)) +
           geom_line(aes(y=waterlevel), color="black", linetype=1, size=1)+
           geom_hline(aes(yintercept=0), color="brown")+
           scale_x_date(date_labels = "%b %Y", date_breaks = "12 month")+
           scale_y_continuous(breaks=seq(-75,75,25))+
           xlab("date")+
           ylab("water table [cm related to surface]")+
           theme_bw()+
           theme(
           axis.title.x=element_text(size=14, vjust=-1, margin=margin(t=10,b=10)),
           axis.title.y=element_text(size=14, margin=margin(r=10,l=10)),
           axis.text.x=element_text(size=10, angle=50,hjust=1),
           axis.text.y=element_text(size=10),
		   plot.margin = unit(c(2.5, 2.5, 2.5, 2.5), "cm"))+
           facet_grid(rows = vars(site))+
           theme(strip.text.x = element_text(size = 11))
ggsave(plot = plot_Pegel, width=21.59, height=27.94, unit="cm", dpi = 300, scale=1, limitsize = FALSE, filename = "waterlevel.pdf")