#Librairies####
library(openxlsx)
library(corrplot)
library(skimr)
library(tidyverse)
library(ggrepel)
library(ggtext)
library(gplots) # heatmap2
library(FactoMineR)
library(factoextra)
library(heatmaply)
library(ggsci)
library(ggpubr)
library(conover.test)
library(cluster)
library(Hmisc)
library(bestNormalize)
library(ggplot2)
library(ggpattern)
library(dplyr)
library(ggbreak)
library(gridExtra)
library(patchwork)
library(scales) 

# Change the working directory to the desired path
setwd("C:/Users/Marion/Documents/phd_lyon/articles/article1/Resultat/")

# Verify that the change has been made
getwd()

#data ####

arcachon <- read.xlsx(xlsxFile = "C:/Users/Marion/Documents/phd_lyon/articles/article1/excel_openacess.xlsx" 
                      ,sheet = "Data" 
                      ,startRow = 1 
                      ,colNames = T 
                      ,detectDates = F 
                      ,rowNames =F 
                      ,skipEmptyRows = T 
                      ,skipEmptyCols = T 
                      ,rows = NULL 
                      ,cols= NULL 
                      ,check.names = T 
                      ,namedRegion = NULL
                      ,na.strings=NA 
                      ,fillMergedCells = T 
)

summary(arcachon)
arcachon$Sample <- as.factor(arcachon$Sample)
arcachon$Site <- as.factor(arcachon$Site)
arcachon$Type <- as.factor(arcachon$Type)
arcachon$Replica <- as.factor(arcachon$Replica)
arcachon$category <- paste(arcachon$Site, arcachon$Type, sep = "-")
arcachon$category <- as.factor(arcachon$category)
arcachon$FN2O_raw  <- as.numeric(arcachon$FN2O_raw)
arcachon$FN2O_day  <- as.numeric(arcachon$FN2O_day)
arcachon$FN2O_eq <- as.numeric(arcachon$FN2O_eq)
arcachon$Offset <- as.numeric(arcachon$Offset)

summary(arcachon)


arcachon$Site <- factor(arcachon$Site, levels = c("Gail","Font", "Garr", "Jacq","Ande", "Comp","Germ")) 
arcachon$category <- factor(arcachon$category, levels = c("Gail-Bare", "Gail-Seagrass","Font-Bare", "Font-Seagrass", "Garr-Bare", "Garr-Seagrass", "Jacq-Bare", "Jacq-Seagrass","Ande-Bare", "Ande-Seagrass","Comp-Bare","Comp-Seagrass","Germ-Bare","Germ-Seagrass")) 

# Pastel color palette for the sites
couleurs_sites <- c(
  "Gail" = "#F8766D",      
  "Garr" = "#7CAE00",    
  "Font" = "#C49A00",      
  "Jacq" = "#00BFC4",      
  "Ande" = "#00A9FF",      
  "Comp" = "#C77CFF",     
  "Germ" = "#FF65AC"       
)

# Pastel color palette for categories
couleurs_categories <- c(
  "Gail-Bare" = "#F8766D",    
  "Gail-Seagrass" = "#F8766D", 
  "Garr-Bare" = "#7CAE00",    
  "Garr-Seagrass" = "#7CAE00", 
  "Font-Bare" = "#C49A00",    
  "Font-Seagrass" = "#C49A00", 
  "Jacq-Bare" = "#00BFC4",    
  "Jacq-Seagrass" = "#00BFC4", 
  "Ande-Bare" = "#00A9FF",    
  "Ande-Seagrass" = "#00A9FF", 
  "Comp-Bare" = "#C77CFF",   
  "Comp-Seagrass" = "#C77CFF", 
  "Germ-Bare" = "#FF65AC",    
  "Germ-Seagrass" = "#FF65AC"  
)


#Wilcoxon test between Bare and Seagrass for CAR ###"

# Test de Wilcoxon global (Bare vs Seagrass)
test <-print(wilcox.test(CAR_17 ~ Type, data = arcachon))
test #p-value = 1.325e-05


#Fig. 3: flux ####

# Calculate the average of FCH4 for each category
mean_data <- arcachon %>%
  filter(!is.na(FCH4_day)) %>%
  group_by(category) %>%
  summarise(mean_FCH4_day = mean(FCH4_day))

# Create the boxplot with the new legends
CH4 <- ggboxplot(arcachon, x = "category", y = "FCH4_day",
               color = "Site",
               add = "jitter", shape = "Site") +
  scale_shape_manual(values = c(16, 17, 18, 19, 20, 21, 22)) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, family = "Arial", size = 12), 
    axis.text.y = element_text(family = "Arial", size = 12),                       
    axis.title.x = element_text(family = "Arial", size = 12),                       
    axis.title.y = element_text(family = "Arial", size = 12),                       
    legend.text = element_text(family = "Arial", size = 12),                        
    legend.title = element_text(family = "Arial", size = 12),                       
    plot.title = element_text(family = "Arial", size = 12, hjust = 0.5)             
  ) +
  guides(shape = "none") +  
  labs(color = "Sites", 
       x = "",        
       y = expression(CH[4] ~ "flux (" * mu * "mol CH"[4] ~ m^-2 ~ d^-1 * ")"))    

# Add the averages in black points
CH4 <- CH4 + geom_point(data = mean_data, aes(x = category, y = mean_FCH4_day), 
                    color = "black", size = 3, shape = 16)

# Show plot
CH4


#Wilcoxon test

p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Font" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Font" & arcachon$Type == "Seagrass",]$FCH4_day)
print(p_value_moy) #p-value = 0.002165


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Seagrass",]$FCH4_day) #pas de différence
print(p_value_moy) #p-value = 1


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Gail" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Gail" & arcachon$Type == "Seagrass",]$FCH4_day) #pas de différence
print(p_value_moy) #p-value = 0.002165

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Ande" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Ande" & arcachon$Type == "Seagrass",]$FCH4_day) #différence car <0.01, p-value = 0.004329 **
print(p_value_moy) #p-value = 0.004329

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Garr" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Garr" & arcachon$Type == "Seagrass",]$FCH4_day) #pas de différence p-value = 0.3095
print(p_value_moy) #p-value = 0.01515

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Comp" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Comp" & arcachon$Type == "Seagrass",]$FCH4_day) #pas de différence
print(p_value_moy) #p-value =0.1775

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Germ" & arcachon$Type == "Bare",]$FCH4_day, arcachon[arcachon$Site == "Germ" & arcachon$Type == "Seagrass",]$FCH4_day) #pas de différence
print(p_value_moy) #p-value =0.06494

# Calculate the average of FN2O_day for each category
mean_data <- arcachon %>%
  filter(!is.na(FN2O_day)) %>%
  group_by(category) %>%
  summarise(mean_FN2O_day = mean(FN2O_day))

# Create the boxplot with the new legendses
N2O <- ggboxplot(arcachon, x = "category", y = "FN2O_day",
               color = "Site",
               add = "jitter", shape = "Site") +
  scale_shape_manual(values = c(16, 17, 18, 19, 20, 21, 22)) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, family = "Arial", size = 12),  
    axis.text.y = element_text(family = "Arial", size = 12),                        
    axis.title.x = element_text(family = "Arial", size = 12),                      
    axis.title.y = element_text(family = "Arial", size = 12),                      
    legend.text = element_text(family = "Arial", size = 12),                      
    legend.title = element_text(family = "Arial", size = 12),                      
    plot.title = element_text(family = "Arial", size = 12, hjust = 0.5)             
  ) +
  guides(shape = "none") +  
  labs(color = "Sites",  
       x = "",       
       y = expression(N[2] ~ "O flux (" * mu * "mol N"[2] * "O m"^-2 ~ d^-1 * ")"))     


# Add the averages in black points
N2O <- N2O + geom_point(data = mean_data, aes(x = category, y = mean_FN2O_day), 
                    color = "black", size = 3, shape = 16)

# show plot
N2O

#Wilcoxon test
p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Font" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Font" & arcachon$Type == "Seagrass",]$FN2O_day)
print(p_value_moy) #p-value = 0.5887


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Seagrass",]$FN2O_day) #pas de différence
print(p_value_moy) #p-value = 0.5887


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Gail" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Gail" & arcachon$Type == "Seagrass",]$FN2O_day) #pas de différence
print(p_value_moy) #p-value = 0.004329

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Ande" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Ande" & arcachon$Type == "Seagrass",]$FN2O_day) #différence car <0.01, p-value = 0.004329 **
print(p_value_moy) #p-value = 0.5368

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Garr" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Garr" & arcachon$Type == "Seagrass",]$FN2O_day) #pas de différence p-value = 0.3095
print(p_value_moy) #p-value = 1

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Comp" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Comp" & arcachon$Type == "Seagrass",]$FN2O_day) #pas de différence
print(p_value_moy) #p-value =0.4127

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Germ" & arcachon$Type == "Bare",]$FN2O_day, arcachon[arcachon$Site == "Germ" & arcachon$Type == "Seagrass",]$FN2O_day) #pas de différence
print(p_value_moy) #p-value =0.007937


#Fig. 4: graph net climate benefit ####


# Calculate the mean, standard deviation, number of observations and standard error
net_climat_benefit_mean <- arcachon %>%
  group_by(Site, Type, category) %>%
  summarise(
    mean_net_climat_benefit = mean(net_climat_benefit, na.rm = TRUE),  
    sd_value = sd(net_climat_benefit, na.rm = TRUE),                  
    n = sum(!is.na(net_climat_benefit)),                              
    se_value = sd_value / sqrt(n)                                  
  ) %>%
  ungroup()  


net_climat_benefit_mean$Site <- factor(net_climat_benefit_mean$Site, levels = c("Gail","Font", "Garr", "Jacq","Ande", "Comp","Germ"))  
net_climat_benefit_mean$category <- factor(net_climat_benefit_mean$category, levels = c("Gail-Bare", "Gail-Seagrass","Font-Bare", "Font-Seagrass", "Garr-Bare", "Garr-Seagrass", "Jacq-Bare", "Jacq-Seagrass","Ande-Bare", "Ande-Seagrass","Comp-Bare","Comp-Seagrass","Germ-Bare","Germ-Seagrass"))  


# Create the chart with patterns for 'Type' and solid colors for 'Site'

ggplot(net_climat_benefit_mean) +
  aes(x = category, y = mean_net_climat_benefit, fill = Site, pattern = Type) +  
  geom_col_pattern(
    aes(pattern_fill = Type, pattern_color = Type),   
    pattern_density = 0.01,    
    pattern_angle = 45,       
    pattern_spacing = 0.02,   
    pattern_fill = "black",   
    pattern_color = "black"   
  ) +
  geom_errorbar(
    aes(ymin = mean_net_climat_benefit - se_value, ymax = mean_net_climat_benefit + se_value)
  ) +
  scale_pattern_manual(values = c("Bare" = "stripe", "Seagrass" = "none")) +  
  scale_y_continuous(limits = c(-20, 110)) +
  labs(
    x = NULL,  # supprime le label de l'axe X
    y = expression(gCO[2] ~ "-eq m"^-2 ~ "y"^-1), 
    fill = "Sites",     
    pattern = "Types"    
  ) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),
    axis.text.x = element_text(angle = 45, hjust = 1, family = "Arial", size = 12),  
    axis.text.y = element_text(family = "Arial", size = 12),                        
    axis.title.y = element_text(family = "Arial", size = 12),                       
    legend.text = element_text(family = "Arial", size = 12),                        
    legend.title = element_text(family = "Arial", size = 12)
  ) +
  scale_x_discrete(labels = NULL) +  # supprime les labels de l’axe X
  guides(
    fill = guide_legend(override.aes = list(pattern = "none")),  # supprime les motifs dans la légende des Sites
    pattern = guide_legend(override.aes = list(pattern_fill = c("white", "gray")))  
  )






## Wilcoxon test


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Font" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Font" & arcachon$Type == "Seagrass",]$net_climat_benefit)
print(p_value_moy) #p-value = 0.002165


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Jacq" & arcachon$Type == "Seagrass",]$net_climat_benefit) 
print(p_value_moy) #p-value = 0.002165


p_value_moy <- wilcox.test(arcachon[arcachon$Site == "Gail" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Gail" & arcachon$Type == "Seagrass",]$net_climat_benefit) 
print(p_value_moy) #p-value = 0.002165

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Ande" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Ande" & arcachon$Type == "Seagrass",]$net_climat_benefit) 
print(p_value_moy) #p-value = 0.002165

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Garr" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Garr" & arcachon$Type == "Seagrass",]$net_climat_benefit)
print(p_value_moy) #p-value = 0.002165

p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Comp" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Comp" & arcachon$Type == "Seagrass",]$net_climat_benefit) 
print(p_value_moy) #p-value = 0.004329
p_value_moy <-wilcox.test(arcachon[arcachon$Site == "Germ" & arcachon$Type == "Bare",]$net_climat_benefit, arcachon[arcachon$Site == "Germ" & arcachon$Type == "Seagrass",]$net_climat_benefit) 
print(p_value_moy) #p-value = 0.002165



#Fig. 5: correlation ####
# Perform the correlation test
cor_test_result <- cor.test(arcachon$net_climat_benefit, arcachon$Fbot_075, method = "spearman") 
# show results
print(cor_test_result)
#rho = -0.2893306 
#p_value = 0.007979

# on est pas en nomal donc pas en pearson donc on ne peut pas avoir d'équation à la droite
# # Ajuster un modèle linéaire
model <- lm(net_climat_benefit ~ Fbot_075, data = arcachon)
# 
# # Extraire les coefficients de l'équation
coeff <- coef(model)
equation <- paste0("y = ", round(coeff[2], 1), "x +", round(coeff[1], 1))
# Calculer le R^2
r_squared <- summary(model)$r.squared
r_squared_text <- paste0("R² = ", round(r_squared, 1))

clim_final <- ggplot(arcachon, aes(x = Fbot_075, y = net_climat_benefit, color = Site, shape = Type)) +
  geom_smooth(aes(group = 1), method = "lm", color = "black") +
  geom_point(size = 5) +
  scale_shape_manual(
    values = c("Bare" = 10, "Seagrass" = 16),
    name = "Types",
    labels = c("Bare", "Seagrass")
  ) +
  scale_color_manual(
    name = "Sites",
    values = couleurs_sites
  ) +
  labs(
    x = "Hydrodynamic forcing at the bottom",
    y = expression("Net climate benefit (gCO"[2]*"-eq m"^-2*"y"^-1*")")
  ) +
  annotate("text", x = 300000, y = 85, label = "P < 0.01", color = "black", size = 5) +
  annotate("text", x = 300000, y = 75, label = "ρ = -0.3", color = "black", size = 5) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),
    axis.text.x = element_text(family = "Arial", size = 15),
    axis.text.y = element_text(family = "Arial", size = 15),
    axis.title.x = element_text(family = "Arial", size = 15),
    axis.title.y = element_text(family = "Arial", size = 15),
    legend.text = element_text(family = "Arial", size = 15),
    legend.title = element_text(family = "Arial", size = 15)
  ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(shape = 16, linetype = 0)),  # Sites en haut
    shape = guide_legend(order = 2)  # Types en dessous
  )

clim_final




#Fig. 6 : SUM ####
# Aggregation of data to obtain the sum by category for the carbon seqne
arcachon_somme <- arcachon %>%
  group_by(Type) %>%
  dplyr::summarize(Somme = sum(CAR_17))

# Check Category Levels
unique(arcachon_somme$Type)

# Plot the barplot with the sums by category

car <- ggplot(arcachon_somme, aes(x = Type, y = Somme, fill = Type)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = format(round(Somme, 1), big.mark = ",", decimal.mark = ".")), 
    vjust = -0.5, color = "black", size = 4
  ) +  
  scale_fill_manual(values = c("Bare" = "coral4", "Seagrass" = "chartreuse3")) +
  scale_y_continuous(
    limits = y_limits,
    labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = ".")
  ) +
  labs(
    x = "",
    y = expression("Sum of CAR (gCO"[2]*"-eq m"^-2*"y"^-1*")")
  ) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),  # axes principaux noirs
    axis.text.x = element_text(family = "Arial", size = 12),      
    axis.text.y = element_text(family = "Arial", size = 12),      
    axis.title.x = element_text(family = "Arial", size = 12),     
    axis.title.y = element_text(family = "Arial", size = 12),     
    plot.title = element_text(family = "Arial", size = 12),       
    legend.position = "none"
  )
# show graph
car

# Aggregating data to obtain the sum by category for flows 
arcachon_somme <- arcachon %>%
  group_by(Type) %>%
  dplyr::summarize(Somme = sum(TotalCO2esource))

# Check Category Levels
unique(arcachon_somme$Type)

# Plot the barplot with the sums by category

# Creating the graph without legend
flux <- ggplot(arcachon_somme, aes(x = Type, y = Somme, fill = Type)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = format(round(Somme, 1), big.mark = ",", decimal.mark = ".")), 
    vjust = -0.5, color = "black", size = 4
  ) +  
  scale_fill_manual(values = c("Bare" = "coral4", "Seagrass" = "chartreuse3")) +
  scale_y_continuous(
    limits = y_limits,
    labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = ".")
  ) +
  labs(
    x = "",
    y = expression("Sum of GHG fluxes (gCO"[2]*"-eq m"^-2*"y"^-1*")")
  ) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),  # axes principaux noirs
    axis.text.x = element_text(family = "Arial", size = 12),      
    axis.text.y = element_text(family = "Arial", size = 12),      
    axis.title.x = element_text(family = "Arial", size = 12),     
    axis.title.y = element_text(family = "Arial", size = 12),     
    plot.title = element_text(family = "Arial", size = 12),       
    legend.position = "none"
  )

# show graph
flux




# Aggregation of data to obtain the sum by category for the net climate benefit balance
arcachon_somme <- arcachon %>%
  group_by(Type) %>%
  dplyr::summarize(Somme = sum(net_climat_benefit))

# Check Category Levels
unique(arcachon_somme$Type)

# Plot the barplot with the sums by category
bilan <- ggplot(arcachon_somme, aes(x = Type, y = Somme, fill = Type)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = format(round(Somme, 1), big.mark = ",", decimal.mark = ".")), 
    vjust = -0.5, color = "black", size = 4
  ) +  
  scale_fill_manual(values = c("Bare" = "coral4", "Seagrass" = "chartreuse3")) +
  scale_y_continuous(
    limits = y_limits,
    labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = ".")
  ) +
  labs(
    x = "",
    y = expression("Sum of net climate benefit (gCO"[2]*"-eq m"^-2*"y"^-1*")")
  ) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),  # axes principaux noirs
    axis.text.x = element_text(family = "Arial", size = 12),      
    axis.text.y = element_text(family = "Arial", size = 12),      
    axis.title.x = element_text(family = "Arial", size = 12),     
    axis.title.y = element_text(family = "Arial", size = 12),     
    plot.title = element_text(family = "Arial", size = 12),       
    legend.position = "none"
  )

# show graph
bilan



# Set Y-axis limits
y_limits <- c(0, 3400)  # Adjust these values according to the data

# Apply Y scale with formatting to all three charts
car <- car + 
  scale_y_continuous(limits = y_limits, labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = "."))

flux <- flux + 
  scale_y_continuous(limits = y_limits, labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = "."))

bilan <- bilan + 
  scale_y_continuous(limits = y_limits, labels = function(x) format(round(x, 1), big.mark = ",", decimal.mark = "."))

# Combines graphics with patchwork (line layout)
combined_plot <- car + flux + bilan + plot_layout(ncol = 3)  

# show graph
combined_plot




#supplementary data ####

arcachon_supp <- read.xlsx(xlsxFile = "C:/Users/Marion/Documents/phd_lyon/articles/article1/excel_openacess.xlsx" 
                      ,sheet = "Fig.sup.hydro" 
                      ,startRow = 1 
                      ,colNames = T 
                      ,detectDates = F 
                      ,rowNames =F 
                      ,skipEmptyRows = T 
                      ,skipEmptyCols = T 
                      ,rows = NULL 
                      ,cols= NULL 
                      ,check.names = T 
                      ,namedRegion = NULL
                      ,na.strings=NA 
                      ,fillMergedCells = T 
)

summary(arcachon_supp)

arcachon_supp$Sample <- as.factor(arcachon_supp$Sample)
arcachon_supp$Site <- as.factor(arcachon_supp$Site)
arcachon_supp$Type <- as.factor(arcachon_supp$Type)
arcachon_supp$Replica <- as.factor(arcachon_supp$Replica)
summary(arcachon_supp)


# Perform the correlation test
cor_test_result <- cor.test(arcachon_supp$FCH4_dayT0, arcachon_supp$FCH4_dayT1, method = "spearman") 
# show results
print(cor_test_result)
#rho = 0.7378741 
#p_value = 0.000702 ***


# Fitting a linear model
model <- lm(FCH4_dayT0 ~ FCH4_dayT1, data = arcachon_supp)

# Extract the coefficients from the equation
coeff <- coef(model)
equation <- paste0("y = ", round(coeff[2], 1), "x +", round(coeff[1], 1))
equation <- expression("T"[0] == 0.9 %.% "T"[1] + 1.5)
# Calcul R^2
r_squared <- summary(model)$r.squared
r_squared_text <- paste0("R² = ", round(r_squared, 1))
time <- ggplot(arcachon_supp, aes(x = FCH4_dayT1, y = FCH4_dayT0)) +
  geom_smooth(method = "lm", color = "black") +
  geom_point(color = "black", size = 2) +
  labs(
    x = expression("CH"[4]*" fluxes T"[1]*" ("*mu*"mol CH"[4]*" m"^-2*" d"^-1*")"), 
    y = expression("CH"[4]*" fluxes T"[0]*" ("*mu*"mol CH"[4]*" m"^-2*" d"^-1*")")
  ) +
  annotate("text", x = 20, y = max(65), label = equation, hjust = 0, vjust = 1, size = 6, family = "Arial", color = "black") +
  annotate("text", x = 20, y = max(60), label = r_squared_text, hjust = 0, vjust = 1.5, size = 6, family = "Arial", color = "black") +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black", size = 0.5),  # axes principaux noirs
    axis.text.x = element_text(family = "Arial", size = 15),
    axis.text.y = element_text(family = "Arial", size = 15),
    axis.title.x = element_text(family = "Arial", size = 15),
    axis.title.y = element_text(family = "Arial", size = 15)
  )
time




#correlation flux and production N-1
arcachon_supp_prod <- read.xlsx(xlsxFile = "C:/Users/Marion/Documents/phd_lyon/articles/article1/excel_openacess.xlsx" 
                           ,sheet = "Sup.prod.n-1" 
                           ,startRow = 1 
                           ,colNames = T 
                           ,detectDates = F 
                           ,rowNames =F 
                           ,skipEmptyRows = T 
                           ,skipEmptyCols = T 
                           ,rows = NULL 
                           ,cols= NULL 
                           ,check.names = T 
                           ,namedRegion = NULL
                           ,na.strings=NA 
                           ,fillMergedCells = T 
)

summary(arcachon_supp_prod)

arcachon_supp_prod$Sample <- as.factor(arcachon_supp_prod$Sample)
arcachon_supp_prod$Site <- as.factor(arcachon_supp_prod$Site)
arcachon_supp_prod$Type <- as.factor(arcachon_supp_prod$Type)
arcachon_supp_prod$Replica <- as.factor(arcachon_supp_prod$Replica)
arcachon_supp_prod$CH4_N_1 <- as.numeric(arcachon_supp_prod$CH4_N_1)

summary(arcachon_supp_prod)

arcachon_supp_prod <- arcachon_supp_prod %>%
  filter(!is.na(CH4_N_1), !is.na(FCH4_day))


cor_test_result_prod <- cor.test(arcachon_supp_prod$FCH4_day, arcachon_supp_prod$CH4_N_1, method = "spearman") 
# show results
print(cor_test_result_prod)
#rho = 0.4063807     
#p_value =  0.0005306***