# ============================================================
# Ecological impacts of temperature reductions in aquatic systems
# Data analysis and visualisation
# ============================================================

# Set the working directory to the folder containing this script
# and the supplementary data file before running the analysis.
# Example:
# setwd("path/to/project/folder")


# ============================================================
# 1. Setup
# ============================================================

# 1.1 Required packages
# Install packages once if needed:
# install.packages(c(
#   "readxl", "tidyverse", "ggplot2", "viridis", "RColorBrewer",
#   "janitor", "maps", "sf", "rnaturalearth",
#   "rnaturalearthdata", "shadowtext"
# ))

library(readxl)
library(tidyverse)
library(ggplot2)
library(viridis)
library(RColorBrewer)
library(janitor)
library(maps)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(shadowtext)


# 1.2 Data import

data <- read_excel("Master Thesis Data.xlsx")


# 1.3 Common plotting theme and colour palette

theme_thesis <- function() {
  
  theme_minimal(base_size = 14) +
    
    theme(
      
      plot.title =
        element_text(
          size = 16,
          face = "bold",
          hjust = 0.5
        ),
      
      axis.title =
        element_text(
          size = 14,
          face = "bold"
        ),
      
      axis.text =
        element_text(
          size = 12,
          colour = "black"
        ),
      
      panel.grid.minor =
        element_blank(),
      
      panel.grid.major.x =
        element_blank(),
      
      legend.position =
        "right",
      
      legend.title =
        element_text(
          face = "bold"
        )
    )
}

cooling_palette <- c(
  "#1b9e77",
  "#7570b3",
  "#d95f02",
  "#e7298a"
)

# ============================================================
# 3.1 Characteristics of the evidence base
# ============================================================

# 3.1.1 Publication trend

pub_year <- data %>%
  distinct(PaperNr, Year) %>%
  count(Year) %>%
  arrange(Year) %>%
  mutate(
    Cumulative = cumsum(n)
  )

#Scale factor
scale_factor <- max(pub_year$Cumulative) / max(pub_year$n)


ggplot(
  pub_year,
  aes(
    x = Year
  )
) +
  
  # Number of studies per year 
  geom_col(
    aes(
      y = n
    ),
    fill = "#2C7FB8",
    width = 0.7
  ) +
  
  # cumulative n umber of studies
  geom_line(
    aes(
      y = Cumulative / scale_factor
    ),
    colour = "black",
    linewidth = 0.75,
    linetype = "dashed"
  ) +
  
  scale_y_continuous(
    name = "Number of Studies",
    
    sec.axis = sec_axis(
      ~ . * scale_factor,
      name = "Cumulative Number of Studies"
    )
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Publication Trend",
    x = "Year"
  )

# Additional descriptive value: year with the highest number of included studies

#year with highest nr of included studies
data %>%
  distinct(PaperNr, Year) %>%
  count(Year) %>%
  filter(n == max(n))

# 3.1.2 Geographic distribution


country_counts <- data %>%
  distinct(PaperNr, Country) %>%
  count(Country)

world <- ne_countries(
  scale = "medium",
  returnclass = "sf"
)

country_counts$Country <- recode(
  country_counts$Country,
  "USA" = "United States of America",
  "UK" = "United Kingdom"
)

world_data <- world %>%
  left_join(
    country_counts,
    by = c("name" = "Country")
  )


# Top 10 countries with most studies + color

top_countries <- world_data %>%
  filter(
    !is.na(n)
  ) %>%
  slice_max(
    n,
    n = 10
  ) %>%
  mutate(
    label_colour =
      ifelse(
        n >= median(n),
        "white",
        "black"
      )
  )


ggplot(world_data) +
  
  geom_sf(
    aes(fill = n),
    color = "white",
    linewidth = 0.1
  ) +
  
  shadowtext::geom_shadowtext(
    data = top_countries,
    aes(
      label = n,
      geometry = geometry,
      colour = label_colour
    ),
    stat = "sf_coordinates",
    size = 3.5,
    bg.colour = "grey20",
    bg.r = 0.12,
    fontface = "bold"
  ) +
  
  scale_colour_identity() +
  
  scale_fill_viridis_c(
    option = "magma",
    trans = "log10",
    na.value = "grey90"
  ) +
  
  theme_void() +
  
  labs(
    title = "Geographic Distribution of Cold-Stress Studies",
    fill = "Studies"
  )

# 3.1.3 Taxonomic coverage

taxa <- data %>%
  distinct(PaperNr, TaxonomicGroup) %>%
  count(TaxonomicGroup)

ggplot(
  taxa,
  aes(
    x = reorder(TaxonomicGroup, n),
    y = n
  )
) +
  
  geom_col(
    fill = "#2C7FB8",
    width = 0.7
  ) +
  
  geom_text(
    aes(
      y = n + 1,
      label = paste0("n = ", n)
    ),
    size = 3,
    fontface = "bold",
    colour = "black",
    hjust = 0
  ) +
  
  coord_flip() +
  
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.15))
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Taxonomic Coverage",
    x = "",
    y = "Number of Studies"
  )

# 3.1.4 Ecosystem representation

eco <- data %>%
  distinct(PaperNr, Ecosystem) %>%
  count(Ecosystem)

ggplot(
  eco,
  aes(
    x = Ecosystem,
    y = n,
    fill = Ecosystem
  )
) +
  
  geom_col(
    width = 0.6
  ) +
  
  geom_text(
    aes(
      label = paste0("n = ", n)
    ),
    vjust = 1.5,
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  theme(
    legend.position = "none"
  ) +
  
  labs(
    title = "Ecosystem Representation",
    x = "",
    y = "Number of Studies"
  )

# 3.1.5 Life-stage representation

# Life stages  
data$LifeStageGrouped <- case_when(
  
  grepl("^Embryo$", data$LifeStage,
        ignore.case = TRUE) ~ "Embryo",
  
  grepl("Larva", data$LifeStage,
        ignore.case = TRUE) ~ "Larva",
  
  grepl("^Juvenile$", data$LifeStage,
        ignore.case = TRUE) ~ "Juvenile",
  
  grepl("^Adult$", data$LifeStage,
        ignore.case = TRUE) ~ "Adult",
  
  TRUE ~ "Others"
)


#  Count Studies per life stage and taxonomic group 
life <- data %>%
  distinct(
    PaperNr,
    LifeStageGrouped,
    TaxonomicGroup
  ) %>%
  count(
    LifeStageGrouped,
    TaxonomicGroup
  )


# Ranking of Life Stages 
life$LifeStageGrouped <- factor(
  life$LifeStageGrouped,
  levels = c(
    "Embryo",
    "Larva",
    "Juvenile",
    "Adult",
    "Others"
  )
)


# Total Number of studies per Life Stage
life_totals <- data %>%
  distinct(
    PaperNr,
    LifeStageGrouped
  ) %>%
  count(
    LifeStageGrouped
  ) 

life$TaxonomicGroup <- factor(
  life$TaxonomicGroup,
  levels = c(
    "Mollusk",
    "Crustacean",
    "Insect",
    "Algae",
    "Fish"
  )
)
# Plot
ggplot(
  life,
  aes(
    x = LifeStageGrouped,
    y = n,
    fill = TaxonomicGroup
  )
) +
  
  geom_col(
    width = 0.6,
    colour = "white",
    linewidth = 0.3
  ) +
  
  geom_text(
    data = life_totals,
    aes(
      x = LifeStageGrouped,
      y = n,
      label = paste0("n = ", n)
    ),
    inherit.aes = FALSE,
    vjust = -0.5,
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_manual(
    values = c(
      "Fish" = "#2C7FB8",
      "Mollusk" = "#F28E2B",
      "Algae" = "#B07AA1",
      "Crustacean" = "#59A14F",
      "Insect" = "#E15759"
    )
  ) +
  
  scale_y_continuous(
    expand = expansion(
      mult = c(0, 0.15)
    )
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Life Stage Representation",
    x = "Life Stage",
    y = "Number of Studies",
    fill = "Organism Group"
  )

# ============================================================
# 3.2 Characteristics of cooling experiments
# ============================================================

# 3.2.1 Distribution of cooling magnitudes

ggplot(
  data,
  aes(
    x = DeltaT
  )
) +
  
  geom_histogram(
    binwidth = 2,
    fill = "#2C7FB8",
    colour = "white",
    linewidth = 0.3
  ) +
  
  scale_x_continuous(
    breaks = seq(
      0,
      max(data$DeltaT, na.rm = TRUE),
      by = 2
    )
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Distribution of Cooling Magnitudes",
    x = expression(Delta*T~"(K)"),
    y = "Frequency"
  )

# ============================================================
# 3.3 Biological responses to cooling
# ============================================================

# 3.3.1 Biological response categories

data$ResponseCategoryShort <- gsub(
  " Response",
  "",
  data$ResponseCategory
)

responses <- data %>%
  count(ResponseCategoryShort)

ggplot(
  responses,
  aes(
    x = reorder(ResponseCategoryShort, n),
    y = n
  )
) +
  
  geom_col(
    fill = "#2C7FB8",
    width = 0.7
  ) +
  
  coord_flip() +
  
  scale_y_continuous(
    breaks = seq(
      0,
      max(responses$n),
      by = 25
    ),
    expand = expansion(
      mult = c(0, 0.15)
    )
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Biological Response Categories",
    x = "",
    y = "Number of Responses"
  )

# 3.3.2 Direction of biological responses

heat <- data %>%
  filter(
    EffectDirection != "N/A"
  ) %>%
  count(
    ResponseCategory,
    EffectDirection
  ) %>%
  mutate(
    text_colour = ifelse(n > median(n), "white", "black")
  )

# Ranking of response categories
heat$ResponseCategory <- factor(
  heat$ResponseCategory,
  levels = rev(c(
    "Survival",
    "Physiological Response",
    "Growth & Development",
    "Behavioral Response",
    "Reproductive Response",
    "Ecological / Population Response"
  ))
)

ggplot(
  heat,
  aes(
    x = EffectDirection,
    y = ResponseCategory,
    fill = n
  )
) +
  
  geom_tile(
    colour = "white"
  ) +
  
  geom_text(
    aes(label = n),
    size = 5,
    fontface = "bold",
    colour = "white"
  ) +
  
  scale_fill_viridis_c(
    option = "plasma"
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Direction of Biological Responses",
    x = "Effect Direction",
    y = "Response Category",
    fill = "Count"
  )

# ============================================================
# 3.4 Influence of cooling magnitude on biological responses
# ============================================================

# 3.4.1 Cooling-magnitude classes and biological responses

data$DeltaClass <- cut(
  data$DeltaT,
  breaks = c(0, 3, 6, 10, Inf),
  labels = c(
    "0–3]",
    "(3–6]",
    "(6–10]",
    ">10"
  ),
  include.lowest = TRUE
)

table(data$DeltaClass)


delta_effect <- data %>%
  filter(
    EffectDirection != "N/A"
  ) %>%
  count(
    DeltaClass,
    EffectDirection
  ) %>%
  group_by(
    DeltaClass
  ) %>%
  mutate(
    Percent = n / sum(n) * 100
  )


delta_effect$EffectDirection <- factor(
  delta_effect$EffectDirection,
  levels = c(
    "Mixed",
    "Negative",
    "Neutral",
    "Positive"
  )
)


ggplot(
  delta_effect,
  aes(
    x = DeltaClass,
    y = Percent,
    fill = EffectDirection
  )
) +
  
  geom_col(
    width = 0.7,
    colour = "white"
  ) +
  
  geom_text(
    aes(
      label = n
    ),
    position =
      position_stack(
        vjust = 0.5
      ),
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  labs(
    title = expression(
      Delta*T~"vs Biological Response"
    ),
    x = expression(
      Delta*T~"Class (K)"
    ),
    y = "% Responses",
    fill = "Effect Direction"
  )

# 3.4.2 Probability of negative responses across cooling magnitudes

# Does increasing cooling magnitude increase the probability 
# of negative biological responses?


# Prepare data

reg_data <- data %>%
  
  filter(
    !is.na(DeltaT),
    EffectDirection != "N/A"
  ) %>%
  
  mutate(
    NegativeResponse =
      ifelse(
        EffectDirection == "Negative",
        1,
        0
      )
  )


# Logistic regression

glm_delta <- glm(
  NegativeResponse ~ DeltaT,
  data = reg_data,
  family = binomial
)


# Model summary

summary(glm_delta)


# Extract beta, p-value and odds ratio

beta_delta <- summary(glm_delta)$coefficients[
  "DeltaT",
  "Estimate"
]


p_delta <- summary(glm_delta)$coefficients[
  "DeltaT",
  "Pr(>|z|)"
]


odds_ratio <- exp(
  beta_delta
)


# Print results

cat(
  "\nLogistic regression: Negative Response ~ DeltaT\n",
  "---------------------------------------------\n",
  "β = ",
  round(beta_delta, 4),
  "\n",
  "p-value = ",
  signif(p_delta, 4),
  "\n",
  "Odds Ratio = ",
  round(odds_ratio, 3),
  "\n"
)


# Text label

stats_label <- paste0(
  "β = ",
  round(beta_delta, 3),
  "\n",
  "p = ",
  signif(p_delta, 3)
)


# Plot

ggplot(
  reg_data,
  aes(
    x = DeltaT,
    y = NegativeResponse
  )
) +
  
  # raw data points
  
  geom_jitter(
    width = 0.15,
    height = 0.05,
    alpha = 0.35,
    size = 2.5,
    colour = "#4D4D4D"
  ) +
  
  
  # logistic regression curve
  
  geom_smooth(
    method = "glm",
    method.args =
      list(
        family = binomial
      ),
    se = TRUE,
    linewidth = 1.3,
    colour = "#D55E00"
  ) +
  
  
  # statistics
  
  annotate(
    "text",
    x = max(reg_data$DeltaT) * 0.65,
    y = 0.2,
    label = stats_label,
    size = 5,
    fontface = "bold"
  ) +
  
  
  scale_x_continuous(
    breaks =
      seq(
        0,
        max(reg_data$DeltaT),
        by = 2
      ),
    expand =
      expansion(
        mult = c(0.02,0.05)
      )
  ) +
  
  
  scale_y_continuous(
    limits = c(0,1),
    breaks = seq(
      0,
      1,
      0.25
    ),
    labels =
      scales::percent
  ) +
  
  
  theme_thesis() +
  
  
  labs(
    title = "Effect of Cooling Magnitude on Negative Responses",
    x = expression(
      Delta*T~"(K)"
    ),
    y = "Probability of Negative Response"
  )

# 3.4.3 Biological responses below and above the 3 K reference value

data$Delta3 <- ifelse(
  data$DeltaT <= 3,
  "≤3 K",
  ">3 K"
)

data$Delta3 <- factor(
  data$Delta3,
  levels = c(
    "≤3 K",
    ">3 K"
  )
)

data$NegativeResponse <- ifelse(
  data$EffectDirection == "Negative",
  1,
  0
)

threshold_model <- glm(
  NegativeResponse ~ Delta3,
  family = binomial,
  data = data
)

summary(threshold_model)

ggplot(
  data,
  aes(
    x = Delta3,
    y = NegativeResponse,
    fill = Delta3
  )
) +
  
  stat_summary(
    fun = mean,
    geom = "col",
    width = 0.6
  ) +
  
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar",
    width = 0.2
  ) +
  
  scale_y_continuous(
    labels = scales::percent
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  theme(
    legend.position = "none"
  ) +
  
  labs(
    title = expression(
      "Effect of " * Delta*T * " Threshold on Neg. Responses"
    ),
    x = expression(
      Delta*T~"Class"
    ),
    y = "Probability of Negative Response"
  )

# ============================================================
# 3.5 Ecosystem-specific biological responses
# ============================================================

# 3.5.1 Freshwater and marine biological responses

eco_response <- data %>%
  
  filter(
    EffectDirection != "N/A"
  ) %>%
  
  count(
    Ecosystem,
    EffectDirection
  ) %>%
  
  group_by(
    Ecosystem
  ) %>%
  
  mutate(
    Percent = n / sum(n) * 100
  )


eco_response$EffectDirection <- factor(
  eco_response$EffectDirection,
  levels = c(
    "Mixed",
    "Negative",
    "Neutral",
    "Positive"
  )
)


ggplot(
  eco_response,
  aes(
    x = Ecosystem,
    y = Percent,
    fill = EffectDirection
  )
) +
  
  geom_col(
    width = 0.6,
    colour = "white"
  ) +
  
  geom_text(
    aes(
      label = n
    ),
    position =
      position_stack(
        vjust = 0.5
      ),
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  labs(
    title = "Ecosystem-specific Biological Responses",
    x = "Ecosystem",
    y = "% Responses",
    fill = "Effect Direction"
  )

# ============================================================
# Appendix B. Supplementary analyses
# ============================================================

# B.1 Cooling type and biological responses

cool_response <- data %>%
  
  filter(
    CoolingType != "N/A",
    EffectDirection != "N/A"
  ) %>%
  
  count(
    CoolingType,
    EffectDirection
  ) %>%
  
  group_by(
    CoolingType
  ) %>%
  
  mutate(
    Percent = n / sum(n) * 100
  )


cool_response$CoolingType <- factor(
  cool_response$CoolingType,
  levels = c(
    "Cold shock",
    "Acute cooling",
    "Gradual cooling",
    "Seasonal cooling"
  )
)


cool_response$EffectDirection <- factor(
  cool_response$EffectDirection,
  levels = c(
    "Mixed",
    "Negative",
    "Neutral",
    "Positive"
  )
)


ggplot(
  cool_response,
  aes(
    x = CoolingType,
    y = Percent,
    fill = EffectDirection
  )
) +
  
  geom_col(
    width = 0.7,
    colour = "white"
  ) +
  
  geom_text(
    aes(
      label = n
    ),
    position = position_stack(
      vjust = 0.5
    ),
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  theme(
    axis.text.x =
      element_text(
        angle = 45,
        hjust = 1
      )
  ) +
  
  labs(
    title = "Cooling Type and Biological Responses",
    x = "Cooling Type",
    y = "% Responses",
    fill = "Effect Direction"
  )

# B.2 Cooling duration and biological responses

duration_response <- data %>%
  
  filter(
    DurationType != "N/A",
    EffectDirection != "N/A"
  ) %>%
  
  count(
    DurationType,
    EffectDirection
  ) %>%
  
  group_by(
    DurationType
  ) %>%
  
  mutate(
    Percent = n / sum(n) * 100
  )


# Ranking of Cooling Duration Categories

duration_response$DurationType <- factor(
  duration_response$DurationType,
  levels = c(
    "Acute exposure",
    "Short-/medium-term exposure",
    "Chronic exposure"
  )
)


# Ranking of Responses

duration_response$EffectDirection <- factor(
  duration_response$EffectDirection,
  levels = c(
    "Mixed",
    "Negative",
    "Neutral",
    "Positive"
  )
)


ggplot(
  duration_response,
  aes(
    x = DurationType,
    y = Percent,
    fill = EffectDirection
  )
) +
  
  geom_col(
    width = 0.6,
    colour = "white"
  ) +
  
  geom_text(
    aes(
      label = n
    ),
    position =
      position_stack(
        vjust = 0.5
      ),
    size = 3,
    fontface = "bold",
    colour = "black"
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  theme(
    axis.text.x =
      element_text(
        angle = 25,
        hjust = 1,
        size = 11
      )
  ) +
  
  labs(
    title = "Cooling Duration and Biological Responses",
    x = "Cooling Duration",
    y = "% Responses",
    fill = "Effect Direction"
  )

# B.3 Exploratory threshold analysis of biological responses at small cooling magnitudes

data$Delta3 <- ifelse(
  data$DeltaT <= 2,
  "≤2 K",
  ">2 K"
)

data$Delta3 <- factor(
  data$Delta3,
  levels = c(
    "≤2 K",
    ">2 K"
  )
)

data$NegativeResponse <- ifelse(
  data$EffectDirection == "Negative",
  1,
  0
)

threshold_model <- glm(
  NegativeResponse ~ Delta3,
  family = binomial,
  data = data
)

summary(threshold_model)

ggplot(
  data,
  aes(
    x = Delta3,
    y = NegativeResponse,
    fill = Delta3
  )
) +
  
  stat_summary(
    fun = mean,
    geom = "col",
    width = 0.6
  ) +
  
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar",
    width = 0.2
  ) +
  
  scale_y_continuous(
    labels = scales::percent
  ) +
  
  scale_fill_brewer(
    palette = "Set2"
  ) +
  
  theme_thesis() +
  
  theme(
    legend.position = "none"
  ) +
  
  labs(
    title = expression(
      "Effect of " * Delta*T * " Threshold on Neg. Responses"
    ),
    x = expression(
      Delta*T~"Class"
    ),
    y = "Probability of Negative Response"
  )