diversity.R

################################################################################

# PanTrop

# This script is a suplement to the paper: Clément et al., in prep.

# MECADEV/ECOTROP: https://mecadev.cnrs.fr/index.php?category/ECOTROP

# Ecological Synthesis Lab (SintECO): https://marcomellolab.wordpress.com

# See README for further info: https://github.com/marmello77/pantrop

################################################################################



# Setting the stage ############################################################


rm(list = ls())

packages <- c("tidyverse","here","conflicted","reshape2","vegan")
for (pkg in packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
    library(pkg, character.only = TRUE)
  }
}

conflicted::conflicts_prefer(dplyr::filter)
conflicted::conflicts_prefer(tidyr::pivot_longer)


# Setup ########################################################################


# Folders
for (folder in c("Figures", "Results")) {
  if (!dir.exists(folder)) dir.create(folder)
}

# Load helper functions used to construct networks (same as in topology.txt)
# Expecting functions like: Net_Construction_Guianas, Net_Construction_Malaysia,
# Net_Construction_Africa
source(here::here("Code","functions.R"))


# Data & networks ##############################################################


## Load raw data (same paths used in the other scripts)
Data_Guianas  <- read.delim(here::here("Data","guianas.txt"), header = TRUE)
Data_Malaysia <- read.csv(here::here("Data","Malaysia_plant.csv"), sep = ";")
Data_Africa   <- read.csv(here::here("Data","Africa.csv"), sep = ";", dec = ",")

## Build edge lists per network using the same constructors
Network_Guianas  <- Net_Construction_Guianas(Data_Guianas)
Network_Malaysia <- Net_Construction_Malaysia(Data_Malaysia)
Network_Africa   <- Net_Construction_Africa(Data_Africa)

## Bind all (columns expected: TreeSp, FrugSp, WebID, Frequency)
net_all <- dplyr::bind_rows(Network_Guianas, Network_Africa,
                            Network_Malaysia) %>%
  dplyr::mutate(
    WebID     = as.character(WebID),
    Frequency = suppressWarnings(as.numeric(Frequency))
  )


# Helpers ######################################################################


# (1) Species richness per network (plants, frugivores, total)
richness_by_network <- function(df) {
  plants <- df %>% dplyr::distinct(WebID, TreeSp) %>%
    dplyr::count(WebID, name = "n_plants")
  frugs  <- df %>% dplyr::distinct(WebID, FrugSp) %>%
    dplyr::count(WebID, name = "n_frugivores")
  dplyr::full_join(plants, frugs, by = "WebID") %>%
    dplyr::mutate(total_species = n_plants + n_frugivores) %>%
    dplyr::arrange(WebID)
}

# (2) Presence/absence matrices per guild
presence_matrix <- function(df, web_col, sp_col) {
  df %>%
    dplyr::distinct(!!rlang::sym(web_col), !!rlang::sym(sp_col)) %>%
    dplyr::mutate(value = 1L) %>%
    reshape2::dcast(stats::reformulate(sp_col, web_col),
                    value.var = "value", fill = 0)
}

# (3) Pairwise Whittaker beta-diversity between networks for a presence matrix
# βW_ij = (gamma_ij / mean(alpha_i, alpha_j)) - 1
whittaker_pairwise <- function(pa_mat) {
  pa_df <- as.data.frame(pa_mat)
  rownames(pa_df) <- pa_df[[1]]
  pa_df <- pa_df[,-1, drop = FALSE]
  # Ensure strictly numeric matrix 
  pa <- data.matrix(pa_df)
  nets <- rownames(pa); n <- length(nets)
  out  <- matrix(0, nrow = n, ncol = n, dimnames = list(nets, nets))
  if (n < 2) return(as.data.frame(out))
  for (i in seq_len(n-1)) {
    for (j in (i+1):n) {
      sp_i <- pa[i,]; sp_j <- pa[j,]
      gamma_ij   <- sum((sp_i + sp_j) > 0)
      alpha_mean <- mean(c(sum(sp_i), sum(sp_j)))
      out[i,j] <- out[j,i] <- (gamma_ij / alpha_mean) - 1
    }
  }
  as.data.frame(out)
}

# (4) Multi-site Whittaker beta across all networks for a presence matrix
# βW_multi = (gamma_all / mean(alpha across networks)) - 1
whittaker_multisite <- function(pa_mat) {
  pa_df <- as.data.frame(pa_mat)
  rownames(pa_df) <- pa_df[[1]]
  pa_df <- pa_df[,-1, drop = FALSE]
  pa <- data.matrix(pa_df)
  gamma_all   <- sum(colSums(pa) > 0)
  alpha_mean  <- mean(rowSums(pa))
  (gamma_all / alpha_mean) - 1
}


# Calculations #################################################################


## Richness
rich_tbl <- richness_by_network(net_all)

## Presence matrices
pa_plants <- presence_matrix(net_all, "WebID", "TreeSp")
pa_frugs  <- presence_matrix(net_all, "WebID", "FrugSp")

## Whittaker beta (pairwise) for plants and frugivores
betaW_plants_pair <- whittaker_pairwise(pa_plants)
betaW_frugs_pair  <- whittaker_pairwise(pa_frugs)

## Whittaker beta (multi-site, all three networks)
betaW_plants_multi <- whittaker_multisite(pa_plants)
betaW_frugs_multi  <- whittaker_multisite(pa_frugs)


# Outputs ######################################################################


## Save tables
readr::write_csv(rich_tbl, here::here("Results",
                                      "diversity_richness_summary.csv"))
readr::write_csv(betaW_plants_pair %>% tibble::rownames_to_column("Network"),
                 here::here("Results","betaW_pairwise_plants.csv"))
readr::write_csv(betaW_frugs_pair %>% tibble::rownames_to_column("Network"),
                 here::here("Results","betaW_pairwise_frugivores.csv"))

## Save a human-readable summary
summary_lines <- c(
  "================ Diversity summary ================",
  "",
  "Richness by network (plants, frugivores, total):",
  capture.output(print(rich_tbl)),
  "",
  "Pairwise Whittaker beta-diversity (plants):",
  capture.output(print(betaW_plants_pair)),
  "",
  "Pairwise Whittaker beta-diversity (frugivores):",
  capture.output(print(betaW_frugs_pair)),
  "",
  sprintf("Multi-site Whittaker beta (plants, 3 networks): %.4f",
          betaW_plants_multi),
  sprintf("Multi-site Whittaker beta (frugivores, 3 networks): %.4f",
          betaW_frugs_multi),
  "==================================================="
)
writeLines(summary_lines, con = here::here("Results","diversity_summary.txt"))

## Also save the workspace for reproducibility
save.image(here::here("Results","Diversity.RData"))

## Console summary
cat(paste(summary_lines, collapse = "\n"))


## Save a CSV summary

# 1) Pairwise βW in long format (one row per network pair)
pair_plants_long <- betaW_plants_pair %>%
  tibble::rownames_to_column("Network1") %>%
  tidyr::pivot_longer(-Network1, names_to = "Network2", values_to = "betaW") %>%
  dplyr::filter(Network1 < Network2) %>%
  dplyr::mutate(guild = "plants")

pair_frugs_long <- betaW_frugs_pair %>%
  tibble::rownames_to_column("Network1") %>%
  tidyr::pivot_longer(-Network1, names_to = "Network2", values_to = "betaW") %>%
  dplyr::filter(Network1 < Network2) %>%
  dplyr::mutate(guild = "frugivores")

beta_pairwise_long <- dplyr::bind_rows(pair_plants_long, pair_frugs_long) %>%
  dplyr::mutate(betaW = round(betaW, 4))

readr::write_csv(beta_pairwise_long, here::here("Results","betaW_pairwise_long.csv"))

# 2) Multi-site βW (one row per guild)
beta_multisite_tbl <- tibble::tibble(
  guild = c("plants", "frugivores"),
  betaW_multisite = round(c(betaW_plants_multi, betaW_frugs_multi), 4)
)
readr::write_csv(beta_multisite_tbl, here::here("Results","betaW_multisite.csv"))

# 3) Single unified summary CSV (richness, pairwise βW, multi-site βW)
summary_csv <- list(
  # Richness by network
  rich_tbl %>%
    dplyr::mutate(section = "richness_network") %>%
    dplyr::rename(Network = WebID),
  # Pairwise βW
  beta_pairwise_long %>%
    dplyr::mutate(section = "beta_pairwise") %>%
    dplyr::rename(value = betaW),
  # Multi-site βW
  beta_multisite_tbl %>%
    dplyr::mutate(section = "beta_multisite")
) |>
  (\(lst) { suppressWarnings(dplyr::bind_rows(lst, .id = NULL)) })()

readr::write_csv(summary_csv, here::here("Results","diversity_summary_table.csv"))


################################################################################