# Clean workspace
rm(list = ls(all = TRUE))

# Load required packages
library(lme4)
library(ggplot2)
library(ggeffects)
library(ggpubr)
library(cowplot)
library(lmerTest)
library(viridis) 
library(scales)
library(gridExtra)

# Read data
dat <- read.csv("Data.csv")

# Calculate saturation deficit
dat$saturation_deficit <- 6.1078 * 10^((7.5 * dat$mean_temp_May) / (237.3 + dat$mean_temp_May)) *
  (1 - dat$mean_humidity_May / 100)

# Models
m1 <- glmer(presenceAnyTypeHerbivory ~ Treatment1 * Elevation + (1 | treeID),
            data = dat, family = binomial)
summary(m1)
m2 <- lmer(leafAssymetry ~ Treatment1 * Elevation + (1 | treeID),
           data = dat)
summary(m2)

# Define elevation breaks and labels as integers (no decimals)
elev_breaks <- seq(38, 1143, by = 221)   # rounded values
elev_labels <- sprintf("%d", elev_breaks)  # format with no decimals

# Function to get ggpredict with elevation factor, and add numeric elev for colors
get_pred_with_elev_num <- function(model) {
  pred <- ggpredict(
    model,
    terms = paste0("Elevation [", min(elev_breaks), ":", max(elev_breaks),
                   " by=", elev_breaks[2] - elev_breaks[1], "]")
  )
  # convert x (factor) to numeric elevation
  pred$elev_num <- as.numeric(as.character(pred$x))
  pred
}

# Get predictions for main effect plots
pred_m1_elev <- get_pred_with_elev_num(m1)
pred_m2_elev <- get_pred_with_elev_num(m2)

# Custom color scale with integer labels on left and color swatches on right
color_scale_elev <- scale_color_viridis_c(
  option = "plasma",
  limits = range(elev_breaks),
  breaks = elev_breaks,
  labels = elev_labels,   # integers only
  name = "Elevation (m)",
  guide = guide_colorbar(
    direction = "vertical",
    label.position = "left",
    title.position = "top"
  )
)

# Calculate cap width in data units (adjust as needed)
cap_width <- 40  # roughly 5 units on the x-axis, adjust to your scale

# For p2
p2 <- ggplot(pred_m1_elev, aes(x = elev_num, y = predicted, group = 1)) +
  geom_line(size = 1) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = elev_num),
                width = 0) +  # no caps here
  geom_segment(aes(x = elev_num - cap_width / 2, xend = elev_num + cap_width / 2,
                   y = conf.low, yend = conf.low, color = elev_num), size = 0.8) +  # lower cap
  geom_segment(aes(x = elev_num - cap_width / 2, xend = elev_num + cap_width / 2,
                   y = conf.high, yend = conf.high, color = elev_num), size = 0.8) +  # upper cap
  geom_point(aes(color = elev_num), size = 3) +
  theme_classic() +
  labs(x = "Elevation (m a.s.l.)", y = "Predicted probability of folivory") +
  color_scale_elev +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = elev_breaks, labels = elev_labels)

# For p5
p5 <- ggplot(pred_m2_elev, aes(x = elev_num, y = predicted, group = 1)) +
  geom_line(size = 1) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = elev_num),
                width = 0) +
  geom_segment(aes(x = elev_num - cap_width / 2, xend = elev_num + cap_width / 2,
                   y = conf.low, yend = conf.low, color = elev_num), size = 0.8) +
  geom_segment(aes(x = elev_num - cap_width / 2, xend = elev_num + cap_width / 2,
                   y = conf.high, yend = conf.high, color = elev_num), size = 0.8) +
  geom_point(aes(color = elev_num), size = 3) +
  theme_classic() +
  labs(x = "Elevation (m a.s.l.)", y = "Predicted value of leaf asymmetry") +
  color_scale_elev +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = elev_breaks, labels = elev_labels)

# Interaction plots - modify to add the same color scale and numeric elevation values
plot_interaction_fixed <- function(model, term1, term2, xlab, ylab, col_lab) {
  pred <- ggpredict(model, terms = c(term1, term2), type = "fixed")
  # Convert factor group to numeric elevation for color mapping
  pred$elev_num <- as.numeric(as.character(pred$group))
  
  ggplot(pred, aes(x = x, y = predicted, colour = elev_num)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.05, size = 0.8) +
    geom_point(size = 3) +
    geom_line(aes(group = group), size = 1) +
    theme_classic() +
    scale_x_discrete(expand = expansion(mult = c(0.2, 0.2))) +
    color_scale_elev +
    labs(x = xlab, y = ylab, colour = col_lab)
}

# Plot interaction plots c and f
p3 <- plot_interaction_fixed(m1, "Treatment1",
                             paste0("Elevation [", min(elev_breaks), ":", max(elev_breaks),
                                    " by=", elev_breaks[2] - elev_breaks[1], "]"),
                             "Treatment", "Predicted probability of folivory", "Elevation")

p6 <- plot_interaction_fixed(m2, "Treatment1",
                             paste0("Elevation [", min(elev_breaks), ":", max(elev_breaks),
                                    " by=", elev_breaks[2] - elev_breaks[1], "]"),
                             "Treatment", "Predicted value of leaf asymmetry", "Elevation")

# Other simple plots for treatment only (plots a and d)
plot_pred_simple <- function(model, term, xlab, ylab, discrete_x = FALSE) {
  gp <- ggpredict(model, terms = term) %>%
    ggplot(aes(x = x, y = predicted, group = ifelse(discrete_x, 1, NA))) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                  width = ifelse(discrete_x, 0.05, 0.05), size = 0.8) +
    geom_point(size = 3) +
    geom_line(size = 1) +
    theme_classic() +
    labs(x = xlab, y = ylab)
  if (discrete_x) gp <- gp + scale_x_discrete(expand = expansion(mult = c(0.2, 0.2)))
  gp
}

p1 <- plot_pred_simple(m1, "Treatment1", "Treatment",
                       "Predicted probability of folivory", TRUE)
p4 <- plot_pred_simple(m2, "Treatment1", "Treatment",
                       "Predicted value of leaf asymmetry", TRUE)

# Remove legends from interaction plots
p3_noleg <- p3 + theme(legend.position = "none")
p6_noleg <- p6 + theme(legend.position = "none")

# Extract legend from interaction plot with integer labels
legend_elev <- cowplot::get_legend(
  p3 + labs(colour = "Elevation (m)") +
    scale_color_viridis_c(
      option = "plasma",
      limits = range(elev_breaks),
      breaks = elev_breaks,
      labels = elev_labels,   # integers
      name = "Elevation (m a.s.l.)",
      guide = guide_colorbar(
        direction = "vertical",
        label.position = "left",
        title.position = "top"
      )
    ) +
    theme(
      legend.position = "right",
      legend.title = element_text(size = 11, margin = margin(b = 11), hjust = 0.5),
      legend.text = element_text(size = 10, hjust = 1)  # right-align legend text
    )
)

# Arrange plots in grid
plots_grid <- ggarrange(
  p1, p2, p3_noleg,
  p4, p5, p6_noleg,
  ncol = 3, nrow = 2,
  labels = c("a", "b", "c", "d", "e", "f"),
  font.label = list(size = 14, face = "bold")
)

# Final combined plot with legend on the right
final_plot <- cowplot::plot_grid(plots_grid, legend_elev, rel_widths = c(1, 0.15))
final_plot

ggsave(
  filename = "Figure 3.png",
  plot = final_plot,
  width = 3500,
  height = 1900,
  units = "px"
)


##################### Mechanisms ####################################

### Q1: What affect the plants
# Fit model
m3 <- lmer(leafAssymetry ~ presenceAnyTypeHerbivory + StandDensityIndex + 
             mean_temp_May + saturation_deficit + (1|treeID), data = dat)
sum_m3 <- summary(m3)
print(sum_m3)

# Extract p-values
pvals <- coef(summary(m3))[, "Pr(>|t|)"]
sig_vec <- pvals < 0.05
names(sig_vec) <- rownames(coef(summary(m3)))
sig_vec <- sig_vec[names(sig_vec) != "(Intercept)"]

# Get marginal effects
eff_herbivory <- ggpredict(m3, terms = "presenceAnyTypeHerbivory")
eff_competition <- ggpredict(m3, terms = "StandDensityIndex")
eff_temperature <- ggpredict(m3, terms = "mean_temp_May")
eff_saturation <- ggpredict(m3, terms = "saturation_deficit")

# Common theme
my_theme <- theme_classic(base_size = 14, base_family = "Arial") +
  theme(
    axis.line = element_line(size = 0.5),
    axis.ticks = element_line(size = 0.5),
    axis.text = element_text(family = "Arial"),
    axis.title = element_text(family = "Arial")
  )

# Calculate common y-axis range and breaks
all_y <- c(
  eff_herbivory$predicted, eff_herbivory$conf.low, eff_herbivory$conf.high,
  eff_competition$predicted, eff_competition$conf.low, eff_competition$conf.high,
  eff_temperature$predicted, eff_temperature$conf.low, eff_temperature$conf.high,
  eff_saturation$predicted, eff_saturation$conf.low, eff_saturation$conf.high
)

ymin <- min(all_y, na.rm = TRUE)
ymax <- max(all_y, na.rm = TRUE)
pad <- 0.02 * (ymax - ymin)
ymin <- ymin - pad
ymax <- ymax + pad
y_breaks <- seq(ymin, ymax, length.out = 5)
y_labels <- number_format(accuracy = 0.001)(y_breaks)

# Plot function
plot_effect <- function(eff, pred_name, x_lab, show_y = TRUE) {
  sig <- sig_vec[pred_name]
  color <- ifelse(sig, "black", "grey50")
  linetype <- ifelse(sig, "solid", "dashed")
  
  p <- ggplot(eff, aes(x = x, y = predicted)) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = color, alpha = 0.15) +
    geom_line(size = 1, color = color, linetype = linetype) +
    scale_y_continuous(
      limits = c(ymin, ymax),
      breaks = y_breaks,
      labels = y_labels
    ) +
    labs(x = x_lab, y = ifelse(show_y, "Predicted value of leaf asymmetry", "")) +
    my_theme
  
  if (!show_y) {
    p <- p + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  }
  p
}

# Create plots
p1 <- plot_effect(eff_herbivory, "presenceAnyTypeHerbivory", "Folivory probablity", show_y = TRUE)
p2 <- plot_effect(eff_competition, "StandDensityIndex", "Stand density index", show_y = FALSE)
p3 <- plot_effect(eff_temperature, "mean_temp_May", "Temperature (°C)", show_y = FALSE)
p4 <- plot_effect(eff_saturation, "saturation_deficit", "VPD (hPa)", show_y = FALSE)

# Arrange with minimal spacing
plot_grid(p1, p2, p3, p4, ncol = 4, align = "hv", axis = "tblr", rel_widths = c(1, 1, 1, 1))

# For plots without y-axis, also remove the y margin completely
remove_y_space <- theme(
  axis.title.y = element_blank(),
  axis.text.y  = element_blank(),
  axis.ticks.y = element_blank(),
  plot.margin  = margin(5, 5, 5, 0)  # no left margin
)

p2 <- p2 + remove_y_space
p3 <- p3 + remove_y_space
p4 <- p4 + remove_y_space

# Arrange in a single row with tight alignment
plot4b<-grid.arrange(p1, p2, p3, p4, ncol = 4, widths = c(1.3, 1, 1, 1))
plot4b

### Q2: What affect the herbivores

m4 <- lmer(presenceAnyTypeHerbivory ~ PredationProbablity + ParasitoidRichness + mean_temp_May + saturation_deficit + (1|treeID), data = dat)

# Print summary with p-values
sum_m4 <- summary(m4)
print(sum_m4)

# Extract p-values from summary (fixed effects)
pvals <- coef(summary(m4))[, "Pr(>|t|)"]

# Create named logical vector for significance (p < 0.05)
sig_vec <- pvals < 0.05
names(sig_vec) <- rownames(coef(summary(m4)))

# Remove intercept from vector
sig_vec <- sig_vec[names(sig_vec) != "(Intercept)"]

# Check sig_vec
print(sig_vec)

# Get marginal effects for each predictor separately
eff_Predation <- ggpredict(m4, terms = "PredationProbablity")
eff_Parasitoid <- ggpredict(m4, terms = "ParasitoidRichness")
eff_temperature <- ggpredict(m4, terms = "mean_temp_May")
eff_saturation <- ggpredict(m4, terms = "saturation_deficit")

# Custom theme with classic base and Arial font
my_theme <- theme_classic(base_size = 14, base_family = "Arial") +
  theme(
    axis.line = element_line(size = 0.5),
    axis.ticks = element_line(size = 0.5),
    axis.text = element_text(family = "Arial"),
    axis.title = element_text(family = "Arial")
  )

# Plot function with significance based line style and color
plot_effect <- function(eff, pred_name, x_lab) {
  sig <- sig_vec[pred_name]
  color <- ifelse(sig, "black", "grey50")
  linetype <- ifelse(sig, "solid", "dashed")
  
  ggplot(eff, aes(x = x, y = predicted)) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = color, alpha = 0.15) +
    geom_line(size = 1, color = color, linetype = linetype) +
    labs(x = x_lab, y = "Herbivory") +
    my_theme
}


# Calculate common y-axis range and breaks 
all_y <- c(
  eff_Predation$predicted, eff_Predation$conf.low, eff_Predation$conf.high,
  eff_Parasitoid$predicted, eff_Parasitoid$conf.low, eff_Parasitoid$conf.high,
  eff_temperature$predicted, eff_temperature$conf.low, eff_temperature$conf.high,
  eff_saturation$predicted, eff_saturation$conf.low, eff_saturation$conf.high
)

ymin <- min(all_y, na.rm = TRUE)
ymax <- max(all_y, na.rm = TRUE)
pad <- 0.02 * (ymax - ymin)
ymin <- ymin - pad
ymax <- ymax + pad
y_breaks <- seq(ymin, ymax, length.out = 5)
y_labels <- number_format(accuracy = 0.001)(y_breaks)

# Plot function with optional y-axis display (like m3)
plot_effect <- function(eff, pred_name, x_lab, show_y = TRUE) {
  sig <- sig_vec[pred_name]
  color <- ifelse(sig, "black", "grey50")
  linetype <- ifelse(sig, "solid", "dashed")
  
  p <- ggplot(eff, aes(x = x, y = predicted)) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = color, alpha = 0.15) +
    geom_line(size = 1, color = color, linetype = linetype) +
    scale_y_continuous(
      limits = c(ymin, ymax),
      breaks = y_breaks,
      labels = y_labels
    ) +
    labs(x = x_lab, y = ifelse(show_y, "Predicted probablity of folivory", "")) +
    my_theme
  
  if (!show_y) {
    p <- p + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  }
  p
}

# Create aligned plots
p1 <- plot_effect(eff_Predation, "PredationProbablity", "Predation probability", show_y = TRUE)
p2 <- plot_effect(eff_Parasitoid, "ParasitoidRichness", "Parasitoid richness", show_y = FALSE)
p3 <- plot_effect(eff_temperature, "mean_temp_May", "Temperature (°C)", show_y = FALSE)
p4 <- plot_effect(eff_saturation, "saturation_deficit", "VPD (hPa)", show_y = FALSE)

# Remove y margins for non-y-axis plots
remove_y_space <- theme(
  axis.title.y = element_blank(),
  axis.text.y  = element_blank(),
  axis.ticks.y = element_blank(),
  plot.margin  = margin(5, 5, 5, 0)
)

p2 <- p2 + remove_y_space
p3 <- p3 + remove_y_space
p4 <- p4 + remove_y_space

# Arrange in a single row, same as m3
plot4a<-grid.arrange(p1, p2, p3, p4, ncol = 4, widths = c(1.3, 1, 1, 1))
plot4a

# Arrange plots in grid
plot4 <- ggarrange(
  plot4a, plot4b,
  ncol = 1, nrow = 2,
  labels = c("a", "b"),
  font.label = list(size = 14, face = "bold")
)

plot4

ggsave(
  filename = "Figure4.png",
  plot = plot4,
  width = 3600,
  height = 2300,
  units = "px"      # size given in pixels
)


############################################################################
#                         Folivory overview plots                          #
############################################################################

# Load packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

# Define elevation bands explicitly
elev_bands <- list(
  b = c(38, 259),
  c = c(260, 480),
  d = c(481, 701),
  e = c(702, 922),
  f = c(923, 1143)
)

# Fixed guild order
guild_levels <- c("Suckers", "Chewers", "Gallers", "Miners", "Skeletonizers", "Rollers")

# Function to summarise guild composition within an elevation band
guild_summary_band <- function(data, elev_range) {
  data %>%
    filter(Elevation >= elev_range[1], Elevation < elev_range[2]) %>%
    summarise(
      Chewers = mean(presenceChewers, na.rm = TRUE) * 100,
      Suckers = mean(presenceSuckers, na.rm = TRUE) * 100,
      Miners = mean(presenceMiners, na.rm = TRUE) * 100,
      Gallers = mean(presenceGallers, na.rm = TRUE) * 100,
      Skeletonizers = mean(presenceSkeletonizers, na.rm = TRUE) * 100,
      Rollers = mean(presenceRollers, na.rm = TRUE) * 100
    ) %>%
    pivot_longer(everything(), names_to = "Guild", values_to = "freq") %>%
    mutate(Guild = factor(Guild, levels = guild_levels))
}

# Base plotting function (with fixed ylim expansion)
plot_guild_bar <- function(summary_data, elev_range) {
  ggplot(summary_data, aes(x = Guild, y = freq, fill = Guild)) +
    geom_col() +
    geom_text(aes(label = sprintf("%.1f%%", freq)), 
              vjust = -0.5, size = 3.5, family = "Arial") +
    labs(x = paste0("Folivore guild (", elev_range[1], "–", elev_range[2], " m asl)"), 
         y = "Folivory occurance (%)") +
    theme_classic() +
    theme(
      legend.position = "none",
      axis.text.x = element_text(angle = 45, hjust = 1, family = "Arial"),
      axis.title = element_text(family = "Arial", size = 12)
    ) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.10)))  # add 10% headroom
}


# Panel a (all data, full elevational range)
guild_summary_all <- guild_summary_band(dat, c(38, 1143))
p3a <- plot_guild_bar(guild_summary_all, c(38, 1143))

# Panels b–f (elevation bands)
p3b <- plot_guild_bar(guild_summary_band(dat, elev_bands$b), elev_bands$b)
p3c <- plot_guild_bar(guild_summary_band(dat, elev_bands$c), elev_bands$c)
p3d <- plot_guild_bar(guild_summary_band(dat, elev_bands$d), elev_bands$d)
p3e <- plot_guild_bar(guild_summary_band(dat, elev_bands$e), elev_bands$e)
p3f <- plot_guild_bar(guild_summary_band(dat, elev_bands$f), elev_bands$f)

# Combine into one figure with labels
fig3 <- ggarrange(
  p3a, p3b, p3c,
  p3d, p3e, p3f,
  ncol = 3, nrow = 2,
  labels = c("a", "b", "c", "d", "e", "f"),
  font.label = list(size = 14, face = "bold")
)

fig3

# Save
ggsave("FigureS2.png", fig3, width = 10, height = 7, dpi = 300)

######################
#Elevation histogram
##################

# Define elevation bands
elev_bands <- list(
  b = c(38.14, 259.21),
  c = c(259.22, 480.29),
  d = c(480.30, 701.37),
  e = c(701.38, 922.45),
  f = c(922.46, 1143.54)
)

# Create a factor variable assigning each patch to an elevational band
patch_elev <- dat %>%
  dplyr::select(patchLabelBETAFOR, Elevation) %>%
  dplyr::distinct() %>%
  dplyr::mutate(
    elev_band = cut(
      Elevation,
      breaks = c(38.14, 259.21, 480.29, 701.37, 922.45, 1143.54),
      labels = c("38–259", "260–480", "481–701", "702–922", "923–1143"),
      include.lowest = TRUE,
      right = TRUE
    )
  )

# Plot as bar chart (counts per band)
ggplot(patch_elev, aes(x = elev_band)) +
  geom_bar(color = "black", fill = "steelblue") +
  theme_classic() +
  labs(x = "Elevation band (m asl)", y = "Number of patches") +
  theme(
    axis.title = element_text(family = "Arial", size = 12),
    axis.text = element_text(family = "Arial", size = 10, angle = 45, hjust = 1)
  )

library(dplyr)
library(ggplot2)

# Count number of patches per band
patch_counts <- patch_elev %>%
  count(elev_band)

# Plot with labels on top of bars
ggplot(patch_counts, aes(x = elev_band, y = n)) +
  geom_col(color = "black", fill = "steelblue") +
  geom_text(aes(label = n), vjust = -0.5, size = 4, family = "Arial") +
  theme_classic() +
  labs(x = "Elevation band (m)", y = "Number of patches") +
  theme(
    axis.title = element_text(family = "Arial", size = 12),
    axis.text = element_text(family = "Arial", size = 10, angle = 45, hjust = 1)
  ) +  scale_y_continuous(expand = expansion(mult = c(0, 0.10)))