rm(list=ls())
# Load required packages
library(ggplot2)
library(readxl)
library(gridExtra)

# Load data for the first plot
ruta_archivo <- "cz_values_interaccionesparana.xlsx"
datos <- read_excel(ruta_archivo)

# Convert variables
datos$STATUS <- as.factor(datos$STATUS)
datos$C <- as.numeric(datos$C)
datos$Z <- as.numeric(datos$Z)

# Kruskal-Wallis for participation coefficient (c) (Olesen et al. 2007)
kruskal.test(C ~ STATUS, data = datos)

# Kruskal-Wallis for within-module degree (z) (Olesen et al. 2007)
kruskal.test(Z ~ STATUS, data = datos)

# First plot
p1 <- ggplot(datos, aes(x = STATUS, y = Z, fill = STATUS)) +
  geom_boxplot(
    width = 0.4, 
    outlier.shape = 21, 
    outlier.fill = "black", 
    outlier.size = 1.8,
    outlier.stroke = 0.5,
    notch = FALSE
  ) +
  scale_fill_manual(values = c("R" = "#B0BEC5", "M" = "#90A4AE")) +
  labs(
    title = "Within-Module Degree (Z) Comparison",
    x = "Species Status (R = Resident, M = Migratory)",
    y = "Within-Module Degree (Z)"
  ) +
  theme_classic(base_size = 15, base_family = "Times New Roman") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 17),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.text = element_text(color = "black"),
    legend.position = "none"
  )

# Load data for the second plot
datos$sform <- as.factor(datos$sform)
# 📌 2. Convert variables to the correct type
datos$sform <- as.factor(datos$sform)  # Convert status to factor
datos$C <- as.numeric(datos$C)  # Ensure numeric values
datos$Z <- as.numeric(datos$Z)

# Kruskal-Wallis C
kruskal.test(C ~ sform, data = datos)

# Kruskal-Wallis for Z
kruskal.test(Z ~ sform, data = datos)

# Second plot
p2 <- ggplot(datos, aes(x = sform, y = Z, fill = sform)) +
  geom_boxplot(
    width = 0.4, 
    outlier.shape = 21, 
    outlier.fill = "black", 
    outlier.size = 1.8,
    outlier.stroke = 0.5,
    notch = FALSE
  ) +
  scale_fill_manual(values = c("ANNUAL HERB" = "#a6bddb", 
                               "CLIMBING PLANT" = "#d0d1e6", 
                               "NN" = "#b8cde6", 
                               "PARASITE BUSH" = "#9ebcda",
                               "SHRUB" = "#8c96c6",
                               "SHRUB/TREE" = "#6e91c6",
                               "TREE" = "#4f79b6")) +
  labs(
    title = "Within-Module Degree (Z) Comparison",
    x = "Life form",
    y = "Within-Module Degree (Z)"
  ) +
  theme_classic(base_size = 15, base_family = "Times New Roman") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 17),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.text = element_text(color = "black"),
    legend.position = "none"
  )

# Combine both plots into one figure with 2 columns
grid.arrange(p1, p2, ncol = 2)

# Save the combined figure
ggsave("combined_graph.png", plot = grid.arrange(p1, p2, ncol = 2),
       width = 16, height = 6, dpi = 600, units = "in")