centrality.R

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

# PanTrop

# This script is a supplement 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

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


# WARNING: this script cannot be used independently. Follow these instructions:
# 1. Run the script topology.R, so you can produce all required results and
# store them as .RData files. The analyses proposed in that script are very
# resource-consuming, as they involve a lot of permutations of large matrices.
# Therefore, it is better to run that script only once and save the results.
# 2. Run the script reload.R to load all .RData files produced by the script
# topology.R. 
# 3. Run this script to analyze centrality.


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



packages <- c("bipartite", "ggmap", "ggplot2", "ggrepel", "ggspatial",
              "here", "igraph", "questionr", "rnaturalearth",
              "rnaturalearthdata", "sf", "taxize", "tidyverse","conflicted",
              "purrr", "broom")

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(bipartite::nullmodel)
conflicted::conflicts_prefer(tidyr::pivot_longer)


if (!requireNamespace("rnaturalearthhires", quietly = TRUE)) {
  message("Package 'rnaturalearthhires' is not installed. Installing from the specified source...")
  
  install.packages(
    "rnaturalearthhires",
    repos = "https://ropensci.r-universe.dev",
    type = "source"
  )
  
  if (requireNamespace("rnaturalearthhires", quietly = TRUE)) {
    message("Package 'rnaturalearthhires' successfully installed.")
    library(rnaturalearthhires)
  } else {
    warning("Failed to install the package 'rnaturalearthhires'.")
  }
  
} else {
  message("Package 'rnaturalearthhires' is already installed.")
  library(rnaturalearthhires)
}


source(here('Code',"functions.R"))


for (folder in c("Figures", "Results")) {
  if (!dir.exists(folder)) {
    dir.create(folder)
    cat(sprintf("Folder '%s' created in the working directory!\n", folder))
  } else {
    cat(sprintf("Folder '%s' already exists.\n", folder))
  }
}



## Guiana Network ####


### Frugivore body Mass dataframe ####

Guiana_data <- Process_Traits_Guiana(
  frug_csv = "Guiana_frug.csv",
  plant_csv = "Guiana_plant.csv",
  sorted_web = Web_Guiana_Sort,
  region_name = "Guiana_shield",
  plant_region_col = "Tree_Class"
)


### Frugivore centrality metrics ####

# Each centrality metric is added to the previously created data frame

Guiana_Frug_Traits <- Guiana_data$Frugivore
Guiana_Frug_Traits <- compute_centralities(Guiana_igraph,
                                           Guiana_Frug_Traits)

write_csv2(Guiana_Frug_Traits,
           file =here("Results",
                      "Guiana_Frug_Traits_Centrality_metrics.csv"))


### Plant centrality metrics ####

Guiana_Plant_Traits <- Guiana_data$Plant


### Network level degree ####

Guiana_Plant_Traits <-
  compute_plant_network_traits(Guiana_igraph, 
                               Guiana_Plant_Traits,Web_Guiana_Sort)

write_csv2(Guiana_Plant_Traits,
           file =here("Results",
                      "Guiana_Plant_Traits_Degree.csv"))


#### Saving the section's workspace ####

save.image(here::here("Results/Guiana_Centrality.RData"))



## Central Africa Network ####



### Frugivore body Mass dataframe ####

Central_Africa_results <- Process_Traits_Central_Africa(
  data = Data_Africa,
  region_name = "Central_Africa",
  study_zones = c("Dja", "Lope", "Makokou",
                  "Belinga", "Loango", "Lekedi", "Moukalaba-Doudou"),
  sorted_web = Web_Central_Africa_Sort
)


### Frugivore centrality metrics ####

Central_Africa_Frug_Traits <- Central_Africa_results$Frugivore

# Each centrality metric is added to the data frame previously created
Central_Africa_Frug_Traits <- 
  compute_centralities(Central_Africa_igraph,
                       Central_Africa_Frug_Traits)

write_csv2(Central_Africa_Frug_Traits,
           file =here("Results",
                      "Central_Africa_Frug_Traits_Centrality_metrics.csv"))


### Plant centrality metrics ####

Central_Africa_Plant_Traits <- Central_Africa_results$Plant


### Network level degree ####

Central_Africa_Plant_Traits <-
  compute_plant_network_traits(Central_Africa_igraph,
                               Central_Africa_Plant_Traits,
                               Web_Central_Africa_Sort)

write_csv2(Central_Africa_Plant_Traits,
           file =here("Results",
                      "Central_Africa_Plant_Traits_Degree.csv"))


#### Saving the section's workspace ####

save.image(here::here("Results/Central_Africa_Centrality.RData"))



## Malaysia Network ####



### Frugivore Body Mass dataframe ####

malaysia_traits <- process_Traits_Malaysia(
  plant_csv = "Malaysia_plant.csv",
  frug_csv = "Malaysia_frug.csv",
  sorted_web = Web_Malaysia_Sort,
  region_name = "Malaysia"
)


### Frugivores centrality metrics ####
# Each centrality metric is added to the data frame previously created

Malaysia_Frug_Traits <- malaysia_traits$Frug_Traits
Malaysia_Frug_Traits <- compute_centralities(Malaysia_igraph,
                                             Malaysia_Frug_Traits)

write_csv2(Malaysia_Frug_Traits,
           file =here("Results",
                      "Malaysia_Frug_Traits_Centrality_metrics.csv"))


### Plants centrality metrics ####

Malaysia_Plant_Traits <- malaysia_traits$Plant_Traits


### Network level degree  ####
Malaysia_Plant_Traits <-
  compute_plant_network_traits(Malaysia_igraph,
                               Malaysia_Plant_Traits,
                               Web_Malaysia_Sort)

write_csv2(Malaysia_Plant_Traits,
           file =here("Results",
                      "Malaysia_Plant_Traits_Degree.csv"))


#### Saving the section's workspace ####


save.image(here::here("Results/Malaysia_Centrality.RData"))



# Centrality summaries #########################################################



## Recovering previous results #################################################


plant_files <- c(
  here::here("Results", "Central_Africa_Plant_Traits_Degree.csv"),
  here::here("Results", "Guiana_Plant_Traits_Degree.csv"),
  here::here("Results", "Malaysia_Plant_Traits_Degree.csv")
)

frug_files <- c(
  here::here("Results", "Central_Africa_Frug_Traits_Centrality_metrics.csv"),
  here::here("Results", "Guiana_Frug_Traits_Centrality_metrics.csv"),
  here::here("Results", "Malaysia_Frug_Traits_Centrality_metrics.csv")
)



### Helper function to read files ###############################################


# with ; separator and comma as decimal mark
read_semicolon <- function(path) {
  readr::read_delim(
    file = path,
    delim = ";",
    locale = readr::locale(decimal_mark = ","),
    show_col_types = FALSE,
    trim_ws = TRUE
  )
}


## Plants ######################################################################


# Functional roles already given in the files
# (Module_hub, Network_hub, Module_connector, Network_connector)

plants_raw <- purrr::map_dfr(
  plant_files,
  ~ read_semicolon(.x) %>%
    mutate(
      region  = dplyr::coalesce(.data$Biogeographical_region, "Unknown"),
      species = .data$TreeSp,
      guild   = "plant"
    ) %>%
    select(region, species, guild, Species_role)
)

## Normalize role labels to "hub" vs "connector"
plants_roles <- plants_raw %>%
  mutate(
    role_simple = case_when(
      Species_role %in% c("Module_hub", "Network_hub") ~ "hub",
      Species_role %in% c("Module_connector", "Network_connector") ~ "connector",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(role_simple)) %>%
  transmute(region, guild, species, role = role_simple, role_method = "given")


## Frugivores ##################################################################


# Infer "hub" and "connector" roles based on centrality metrics

# Default criteria: top 10% in eigenvector = hub; top 10% in betweenness =
# connector. You can adjust the threshold in 'q_thr' (e.g., 0.95 for top 5%).

q_thr <- 0.90


### Helper: robust coercion to numeric #########################################


coerce_num <- function(x) {
  if (is.character(x) || is.factor(x)) {
    readr::parse_number(as.character(x), locale = readr::locale(decimal_mark = ","))
  } else {
    suppressWarnings(as.numeric(x))
  }
}


frug_raw <- purrr::map_dfr(
  frug_files,
  ~ read_semicolon(.x) %>%
    mutate(
      region  = dplyr::coalesce(.data$Biogeographical_region, "Unknown"),
      species = .data$FrugSp,
      guild   = "frugivore"
    ) %>%
    ## Coerce relevant columns to numeric
    mutate(
      Body.mass         = coerce_num(Body.mass),
      degree            = coerce_num(degree),
      normalized_degree = coerce_num(normalized_degree),
      betweeness        = coerce_num(betweeness),
      closeness         = coerce_num(closeness),
      eigen_vector      = coerce_num(eigen_vector)
    ) %>%
    select(region, species, guild, Body.mass, degree, normalized_degree, betweeness, closeness, eigen_vector)
)


# Apply quantile thresholds by region to assign roles

frug_roles <- frug_raw %>%
  group_by(region) %>%
  mutate(
    thr_hub = quantile(eigen_vector, probs = q_thr, na.rm = TRUE, type = 7),
    thr_con = quantile(betweeness,   probs = q_thr, na.rm = TRUE, type = 7),
    is_hub       = ifelse(!is.na(eigen_vector) & eigen_vector >= thr_hub, TRUE, FALSE),
    is_connector = ifelse(!is.na(betweeness)   & betweeness   >= thr_con, TRUE, FALSE),
    role_inferred = case_when(
      is_hub &  is_connector ~ "hub+connector",
      is_hub & !is_connector ~ "hub",
      !is_hub & is_connector ~ "connector",
      TRUE ~ NA_character_
    )
  ) %>%
  ungroup() %>%
  filter(!is.na(role_inferred)) %>%
  transmute(region, guild, species, role = role_inferred, role_method = "inferred")


## Summary 1: by network #######################################################


summary_by_network <- bind_rows(plants_roles, frug_roles) %>%
  arrange(region, guild, role, species)

## Save:

readr::write_csv(
  summary_by_network,
  here::here("Results", "centrality_roles_by_network_summary.csv")
)


## Summary 2: species in multiple networks #####################################


# Expand "hub+connector" into separate rows for clarity by region

summary_multi_network <- summary_by_network %>%
  tidyr::separate_rows(role, sep = "\\+") %>%
  group_by(guild, species) %>%
  filter(dplyr::n_distinct(region) > 1) %>%
  summarise(
    n_regions = dplyr::n_distinct(region),
    roles_by_region = paste0(region, ":", role, collapse = " | "),
    .groups = "drop"
  ) %>%
  arrange(guild, desc(n_regions), species)


# Save:

readr::write_csv(
  summary_multi_network,
  here::here("Results", "centrality_species_multi_network_summary.csv")
)


## Summary 3: body mass vs centrality (frugivores) ##############################



### 1) Ensure frug_raw and frug_for_tests exist ################################


if (!exists("frug_raw")) {
  stop("Object 'frug_raw' not found. It must be created earlier from the frugivore CSVs.")
}

frug_for_tests <- frug_raw %>%
  mutate(
    log_mass = suppressWarnings(log10(as.numeric(Body.mass)))
  )

### 2) Identify which centrality metrics are actually available ################


wanted_metrics <- c("degree","normalized_degree","betweeness","closeness","eigen_vector")
metrics <- intersect(wanted_metrics, names(frug_for_tests))
if (length(metrics) == 0) {
  warning("No centrality metrics found among columns: ", paste(wanted_metrics, collapse=", "))
  ## create an empty tibble but ensure the object exists
  mass_centrality_tests <- tibble(
    level=character(), region=character(), metric=character(),
    rho=numeric(), p.value=numeric(), n=integer(), interpretation=character()
  )
} else {
  
  
### 3) Safe Spearman correlation function ######################################

  
  cor_spearman_safe <- function(x, y) {
    df <- tibble(x=x, y=y) %>%
      filter(is.finite(x), is.finite(y))
    if (nrow(df) < 3) {
      return(tibble(rho = NA_real_, p.value = NA_real_, n = nrow(df)))
    }
    ct <- suppressWarnings(cor.test(df$x, df$y, method = "spearman", exact = FALSE))
    tibble(rho = unname(ct$estimate), p.value = unname(ct$p.value), n = nrow(df))
  }
  
  
### 4) Compute by region (avoiding reframe/cur_data for compatibility) #########

  
  regions <- sort(unique(frug_for_tests$region))
  by_region_list <- list()
  
  for (rg in regions) {
    df_rg <- frug_for_tests %>% filter(region == rg)
    for (met in metrics) {
      res <- cor_spearman_safe(df_rg$log_mass, df_rg[[met]])
      by_region_list[[length(by_region_list) + 1]] <- tibble(
        level = "by_region",
        region = rg,
        metric = met,
        rho = res$rho,
        p.value = res$p.value,
        n = res$n
      )
    }
  }
  by_region <- dplyr::bind_rows(by_region_list)
  
  
### 5) Compute overall correlations (all networks combined) ####################

  
  overall_list <- lapply(metrics, function(met) {
    res <- cor_spearman_safe(frug_for_tests$log_mass, frug_for_tests[[met]])
    tibble(
      level = "overall",
      region = "ALL",
      metric = met,
      rho = res$rho,
      p.value = res$p.value,
      n = res$n
    )
  })
  overall <- dplyr::bind_rows(overall_list)
  
  
### 6) Merge results and add interpretation ####################################

  
  mass_centrality_tests <- bind_rows(by_region, overall) %>%
    mutate(
      interpretation = dplyr::case_when(
        is.na(rho) ~ "insufficient sample",
        p.value < 0.05 & rho > 0 ~ "positive (significant)",
        p.value < 0.05 & rho < 0 ~ "negative (significant)",
        TRUE ~ "not significant"
      )
    ) %>%
    arrange(level, region, metric)
}


### 7) Export (keeping your preferred style) ###################################


readr::write_csv(
  mass_centrality_tests,
  here::here("Results", "centrality_mass_vs_centrality_spearman_summary.csv")
)
message("✔ Saved: Results/centrality_mass_vs_centrality_spearman_summary.csv")



# Optional: friendly console messages ##########################################



message("# Summary 1 stored in object: summary_by_network (roles by network)")
message("# Summary 2 stored in object: summary_multi_network (multi-network species and roles)")
message("# Summary 3 stored in object: mass_centrality_tests (Spearman log-mass vs centrality)")



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