############################################################
################## Supplementary code ########################
############## Hill-number diversity analyses #####################
############################################################

library(dplyr)
library(tidyr)
library(hillR)
library(iNEXT)
library(stringr)
library(ggplot2)

############################################################
## 1. Alpha diversity per haul
############################################################

q_values <- c(0, 1, 2)

alpha_diversity <- sapply(q_values, function(q) {
    hillR::hill_taxa(comm_matrix, q = q, MARGIN = 1)})

colnames(alpha_diversity) <- paste0("Dq_", q_values)

alpha_diversity <- as.data.frame(alpha_diversity) %>%
  mutate(haul_id = rownames(comm_matrix)) %>%
  left_join(haul_metadata, by = "haul_id") %>%
  relocate(haul_id, condition)

alpha_diversity

############################################################
## 2. Summary of alpha diversity by OMZ condition
############################################################

alpha_summary <- alpha_diversity %>%
  pivot_longer(cols = starts_with("Dq_"),
    names_to = "q",
    values_to = "diversity") %>%
  mutate(q = recode(q,
      "Dq_0" = "q = 0",
      "Dq_1" = "q = 1",
      "Dq_2" = "q = 2"),
    q = factor(q, levels = c("q = 0", "q = 1", "q = 2"))) %>%
  group_by(condition, q) %>%
  summarise(n = sum(!is.na(diversity)),
    median = median(diversity, na.rm = TRUE),
    Q1 = quantile(diversity, 0.25, na.rm = TRUE),
    Q3 = quantile(diversity, 0.75, na.rm = TRUE),
    IQR = IQR(diversity, na.rm = TRUE),
    min = min(diversity, na.rm = TRUE),
    max = max(diversity, na.rm = TRUE),
    .groups = "drop")

alpha_summary

############################################################
## 3. Multiplicative partitioning of taxonomic diversity
############################################################

partition_q0 <- hillR::hill_taxa_parti(comm_matrix, q = 0)
partition_q1 <- hillR::hill_taxa_parti(comm_matrix, q = 1)
partition_q2 <- hillR::hill_taxa_parti(comm_matrix, q = 2)

diversity_partition <- rbind(
  partition_q0,
  partition_q1,
  partition_q2)

diversity_partition

############################################################
## 4. Pairwise beta diversity among hauls
############################################################

pairwise_beta_q0 <- hillR::hill_taxa_parti_pairwise(comm_matrix, q = 0)
pairwise_beta_q1 <- hillR::hill_taxa_parti_pairwise(comm_matrix, q = 1)
pairwise_beta_q2 <- hillR::hill_taxa_parti_pairwise(comm_matrix, q = 2)

pairwise_beta <- rbind(
  pairwise_beta_q0,
  pairwise_beta_q1,
  pairwise_beta_q2
)

pairwise_beta_table <- pairwise_beta %>%
  mutate(
    q = as.integer(q),
    pairwise = paste(site1, "-", site2)
  ) %>%
  transmute(
    pairwise,
    q,
    alpha = TD_alpha,
    gamma = TD_gamma,
    beta = TD_beta,
    local_similarity = local_similarity
  )

pairwise_beta_table

############################################################
## 5. Summary of pairwise beta diversity by comparison type
############################################################

pairwise_metadata <- pairwise_beta_table %>%
  separate(
    pairwise,
    into = c("haul_1", "haul_2"),
    sep = " - ",
    remove = FALSE
  ) %>%
  left_join(
    haul_metadata %>%
      rename(haul_1 = haul_id, condition_1 = condition),
    by = "haul_1"
  ) %>%
  left_join(
    haul_metadata %>%
      rename(haul_2 = haul_id, condition_2 = condition),
    by = "haul_2"
  ) %>%
  mutate(
    comparison = case_when(
      condition_1 == "Outside" & condition_2 == "Outside" ~ "Outside-Outside",
      condition_1 == "Inside" & condition_2 == "Inside" ~ "Inside-Inside",
      condition_1 != condition_2 ~ "Inside-Outside"
    ),
    comparison = factor(
      comparison,
      levels = c("Outside-Outside", "Inside-Outside", "Inside-Inside")
    )
  )

beta_summary <- pairwise_metadata %>%
  group_by(comparison, q) %>%
  summarise(
    n = sum(!is.na(beta)),
    median = median(beta, na.rm = TRUE),
    Q1 = quantile(beta, 0.25, na.rm = TRUE),
    Q3 = quantile(beta, 0.75, na.rm = TRUE),
    IQR = IQR(beta, na.rm = TRUE),
    min = min(beta, na.rm = TRUE),
    max = max(beta, na.rm = TRUE),
    .groups = "drop"
  )

beta_summary

############################################################
## 6. Sample coverage per haul
############################################################

community_list <- lapply(
  seq_len(nrow(comm_matrix)),
  function(i) as.numeric(comm_matrix[i, ])
)

names(community_list) <- rownames(comm_matrix)

coverage_info <- iNEXT::DataInfo(
  community_list,
  datatype = "abundance"
)

coverage_info <- coverage_info %>%
  rename(haul_id = Assemblage) %>%
  left_join(haul_metadata, by = "haul_id")

coverage_info

############################################################
## 7. Coverage-standardized diversity estimates
############################################################

coverage_reference <- 0.80

coverage_standardized_diversity <- iNEXT::estimateD(
  community_list,
  q = c(0, 1, 2),
  datatype = "abundance",
  base = "coverage",
  level = coverage_reference
)

coverage_standardized_diversity

############################################################
########## Bayesian models and model comparison ################
############################################################

library(brms)
library(loo)
library(dplyr)

############################################################
## 1. Species richness models
############################################################

richness_null <- brm( richness ~ swept_area_km2 + year + gp(longitude, latitude),
  data = richness_data,
  family = poisson(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

richness_OMZ_condition <- brm(richness ~ swept_area_km2 + year + OMZ_condition +
    gp(longitude, latitude),
  data = richness_data,
  family = poisson(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

richness_haul_depth <- brm(richness ~ swept_area_km2 + year + haul_depth +
    gp(longitude, latitude),
  data = richness_data,
  family = poisson(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

richness_OMZ_depth <- brm(richness ~ swept_area_km2 + year + OMZ_upper_boundary_depth +
    gp(longitude, latitude),
  data = richness_data,
  family = poisson(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

richness_oxygen <- brm(richness ~ swept_area_km2 + year + dissolved_oxygen +
    gp(longitude, latitude),
  data = richness_data,
  family = poisson(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

############################################################
## 2. LOOIC comparison: species richness
############################################################

loo_richness_null <- loo(richness_null,
  moment_match = TRUE)

loo_richness_OMZ_condition <- loo(richness_OMZ_condition,
  moment_match = TRUE)

loo_richness_haul_depth <- loo(richness_haul_depth,
  moment_match = TRUE)

loo_richness_OMZ_depth <- loo(richness_OMZ_depth,
  moment_match = TRUE)

loo_richness_oxygen <- loo(richness_oxygen,
  moment_match = TRUE)

richness_loo_list <- list("Null" = loo_richness_null,
  "OMZ condition" = loo_richness_OMZ_condition,
  "Haul depth" = loo_richness_haul_depth,
  "OMZ upper-boundary depth" = loo_richness_OMZ_depth,
  "Dissolved oxygen" = loo_richness_oxygen)

richness_loo_table <- data.frame(model = names(richness_loo_list),
  LOOIC = sapply(richness_loo_list,
    function(x) x$estimates["looic", "Estimate"]),
  SE = sapply(richness_loo_list,
    function(x) x$estimates["looic", "SE"])) %>%
  mutate(delta_LOOIC = LOOIC - min(LOOIC),
    weight = exp(-0.5 * delta_LOOIC),
    weight = weight / sum(weight)) %>%
  arrange(delta_LOOIC)

richness_loo_table

############################################################
## 3. Abundance models
############################################################

abundance_null <- brm(abundance ~ swept_area_km2 + year +
    gp(longitude, latitude) +
    (1 | OTU) + (1 | haul_id),
  data = abundance_data,
  family = zero_inflated_negbinomial(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

abundance_OMZ_condition <- brm(
  bf(abundance ~ swept_area_km2 + OMZ_condition + year +
      gp(longitude, latitude) +
      (1 | OTU) + (1 | haul_id),
    zi ~ OMZ_condition),
  data = abundance_data,
  family = zero_inflated_negbinomial(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

abundance_haul_depth <- brm(
  bf(abundance ~ swept_area_km2 + haul_depth + year +
      gp(longitude, latitude) +
      (1 | OTU) + (1 | haul_id),
    zi ~ haul_depth),
  data = abundance_data,
  family = zero_inflated_negbinomial(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

abundance_OMZ_depth <- brm(
  bf(abundance ~ swept_area_km2 + OMZ_upper_boundary_depth + year +
      gp(longitude, latitude) +
      (1 | OTU) + (1 | haul_id),
    zi ~ OMZ_upper_boundary_depth),
  data = abundance_data,
  family = zero_inflated_negbinomial(),
  chains = 4,
  cores = 4,
  iter = 10000,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

abundance_oxygen <- brm(
  bf(abundance ~ swept_area_km2 + dissolved_oxygen + year +
      gp(longitude, latitude) +
      (1 | OTU) + (1 | haul_id),
    zi ~ dissolved_oxygen),
  data = abundance_data,
  family = zero_inflated_negbinomial(),
  chains = 4,
  cores = 4,
  iter = 7500,
  control = list(adapt_delta = 0.95),
  save_pars = save_pars(all = TRUE))

############################################################
## LOOIC comparison: abundance
############################################################

## The abundance candidate models were compared using the same
## PSIS-LOO workflow described above for richness models.
## The resulting LOOIC table was constructed from the abundance
## model objects using make_LOOIC_table().