In karinorman/biodivTS: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE)

library(dplyr)
library(tidyr)
library(ggplot2)

load(here::here("data/model_output.rda"))
load(here::here("data/model_data.rda"))

Types of Biodiversity Change

wide_table <- model_output %>%
  filter(!is.na(study_id)) %>%
  filter(p.adjust < 0.05, term == "year_scaled", model_id == "base") %>% 
  dplyr::select(study_id, metric, estimate) %>% 
  pivot_wider(names_from = "metric", values_from = "estimate") %>% 
  dplyr::select(., -starts_with("CWM_")) %>% 
  dplyr::select(study_id, S, Jaccard_base, everything()) %>%
  left_join(model_data %>% 
              dplyr::select(study_id, duration) %>%
              distinct() %>%
              group_by(study_id) %>% 
              mutate(duration = max(duration)) %>% distinct() %>% 
              mutate(study_id = as.character(study_id)))

Add decriptive columns to table

fun_change_table <- wide_table %>%
  mutate(rich_change = ifelse(!is.na(S), "Richness Change", "No Richness Change"),
         turnover = ifelse(!is.na(Jaccard_base), "Turnover", "No Turnover"),
         fun_change = case_when(
           if_all(c(FRic, SES_FRic, FDiv, SES_FDiv, FEve, SES_FEve), ~is.na(.)) ~ "No Functional Change",
           .default = "Functional Change"
         ))

labels = fun_change_table %>%
  dplyr::select(rich_change, turnover, fun_change) %>%
  distinct() %>%
  mutate(group = LETTERS[1:8])

fun_change_table <- fun_change_table %>%
  left_join(labels)

fun_change_table %>% 
  count(rich_change, turnover, fun_change, group) 

#usethis::use_data(fun_change_table)

Let's split up the corrected and uncorrected functional metrics

ses_wide_table <-  wide_table %>%
  dplyr::select(-FRic, -FDiv, -FEve) %>%
  mutate(rich_change = ifelse(!is.na(S), TRUE, FALSE),
         turnover = ifelse(!is.na(Jaccard_base), TRUE, FALSE),
         fun_change = case_when(
           if_all(c(SES_FRic, SES_FDiv, SES_FEve), ~is.na(.)) ~ FALSE,
           .default = TRUE
         ))

ses_wide_table %>% 
  count(rich_change, turnover, fun_change)

uncor_wide_table <-  wide_table %>%
  dplyr::select(-SES_FRic, -SES_FDiv, -SES_FEve) %>%
  mutate(rich_change = ifelse(!is.na(S), TRUE, FALSE),
         turnover = ifelse(!is.na(Jaccard_base), TRUE, FALSE),
         fun_change = case_when(
           if_all(c(FRic, FDiv, FEve), ~is.na(.)) ~ FALSE,
           .default = TRUE
         ))

uncor_wide_table %>% 
  count(rich_change, turnover, fun_change)

Look at distribution of durations

duration_table <- fun_change_table %>%
  dplyr::select(study_id, rich_change, turnover, fun_change, group) %>%
  left_join(model_data %>% 
              dplyr::select(study_id, rarefyID, duration) %>% 
              distinct())

plot_table <- duration_table %>%
  filter(group %in% c("A", "B", "C", "D"))

group_types <- as_labeller(
  c("A" = "No richness change or turnover, functional change", "B" = "No turnover, richness and functional change", "C" = "No change", "D" = "Turnover, richness and functional change")
)

my_pal <- c("#d9be3b", "#77976e", "#086788", "#cc8214", "#bfac88", "#99adbf")
group_colors <- c("A" = "#d9be3b", "B" = "#77976e", "C" = "#086788", "D" = "#cc8214")

ggplot(plot_table, aes(x = duration, fill = group)) +
  geom_histogram(data = dplyr::select(plot_table, -group), fill = "grey", alpha = .5) +
  geom_histogram(colour = "black") +
  facet_wrap(~ group, labeller = group_types) +
  ylab("count of time series") +
  xlab("time series duration") +
  guides(fill = FALSE) +  # to remove the legend
  scale_fill_manual(values = group_colors) +
  theme_classic() 

ggsave(here::here("figures/duration_hist.jpeg"), bg = "white")

library(kableExtra)

knitr::include_graphics(here::here("figures/duration_hist.jpeg"))

fun_change_table %>% 
  dplyr::select(-group) %>%
  count(rich_change, turnover, fun_change) %>%
  pivot_wider(names_from = rich_change, values_from = n) %>%
  select(fun_change, everything()) %>%
  arrange(fun_change) %>%
  mutate(`No Richness Change` = cell_spec(`No Richness Change`,
                                           color = case_when(
                                             `No Richness Change` == 12 ~ "#d9be3b",
                                             `No Richness Change` == 11 ~ "#086788",
                                             .default = "black"))) %>%
  mutate(`Richness Change` = cell_spec(`Richness Change`,
                                           color = case_when(
                                             `Richness Change` == 6 ~ "#77976e",
                                             `Richness Change` == 7 ~ "#cc8214",
                                             .default = "black"))
         ) %>%
  kable(format = "latex", col.names = c("", "", "No Richness Change", "Richness Change"), 
        booktabs = TRUE, align = "c", escape = FALSE) %>%
  kable_styling("striped") %>%
  # row_spec(0, bold = TRUE) %>%
  # column_spec(1:2, bold = TRUE) %>%
  collapse_rows(columns = 1, valign = "middle") #%>%
  #save_kable(here::here("figures/change_table.jpeg"), latex_header_includes = c("\\\\usepackage{booktabs}"), keep_tex = TRUE)

Let's look at more specific combinations

load(here::here("data/meta_clean.rda"))

study_meta <- meta_clean %>% 
  select(-rarefyid) %>%
  distinct()

#remove NA entries if another country was assigned
na_studies <- study_meta %>% 
  group_by(study_id) %>% 
  filter(n() > 1, is.na(country)) 

study_meta <- setdiff(study_meta, na_studies)

# change_groups <- wide_table %>% 
#   select(study_id, S, Jaccard_base, starts_with("SES")) %>%
#   mutate(across(!study_id, ~case_when(
#     . > 0 ~ "pos",
#     . < 0 ~ "neg",
#     .default = "0"
#   ))) 

all_metrics <- c("S", "Jaccard_base", "FRic", "FEve", "FDiv", "SES_FRic", "SES_FEve", "SES_FDiv")

change_groups <- wide_table %>% 
  select(study_id,all_of(all_metrics)) %>%
  mutate(across(!study_id, ~case_when(
    . > 0 ~ "pos",
    . < 0 ~ "neg",
    .default = "0"
  ))) 

change_group_count <- change_groups %>%
  select(-study_id, -starts_with("SES_")) %>%
  group_by(across(everything())) %>%
  count() 

change_group_count_ses <- change_groups %>%
  select(-c(study_id, FRic, FEve, FDiv)) %>%
  group_by(across(everything())) %>%
  count() 

change_group_count_complete <- change_groups %>%
    select(-c(study_id)) %>%
    group_by(across(everything())) %>%
    count() %>%
  mutate(change_description = case_when(
    S == "pos" ~ "increase in richness, not beyond random expectation", 
    S == "neg" ~ "loss of redundancy", 
    S == "0" & Jaccard_base == "neg" ~ "changes in functional space due to turnover",
    n == 7 ~ "no change", 
    .default = NA
  ))

change_group_ses_desc <- change_group_count_ses %>%
   mutate(change_description = case_when(
    S == "pos" & n != "7" ~ "increase in species richness", 
    S == "neg" ~ "loss of redundancy", 
    S == "0" ~ "functional change",
    n == "7" ~ "no change", 
    .default = NA
  ))

usethis::use_data(change_group_ses_desc)

Look at CWM's

cwm_counts <- model_output %>% 
  filter(!is.na(study_id), model_id == "base", p.adjust < 0.05, 
         !metric %in% c("S", "FRic", "FEve", "FDiv", "SES_FRic", "SES_FEve", "SES_FDiv", "Jaccard_base")) %>% 
  select(metric, study_id) %>% 
  count(metric)

# get study-level significant CWM counts with their other change characteristics
cwm_counts %>% 
  full_join(change_groups) %>% 
  select(-c("FRic", "FEve", "FDiv")) %>% 
  View()

# change <-  change_groups %>% 
#   mutate(study_id = as.integer(study_id)) %>%
#   filter(if_any(all_of(c("S", "Jaccard_base", "SES_FRic", "SES_FEve", "SES_FDiv")), ~. != 0))
# 
# no_change <- study_meta %>% filter(!study_id %in% change$study_id, study_id %in% model_output$study_id)
# 
# redundancy <- change_groups %>% 
#   mutate(study_id = as.integer(study_id)) %>%
#   filter(S == "neg", SES_FRic == "0", SES_FDiv == "0", SES_FEve == "0") %>%
#   left_join(study_meta)
# 
# redundancy_gain <- change_groups %>% 
#   mutate(study_id = as.integer(study_id)) %>%
#   filter(S == "pos", SES_FRic == "0", SES_FDiv == "0", SES_FEve == "0") %>%
#   left_join(study_meta)
# 
#   
# evenness <- change_groups %>% 
#   mutate(study_id = as.integer(study_id)) %>%
#   filter(across(!study_id & !SES_FEve, ~. == 0)) %>%
#   left_join(study_meta)
# 
# increase <- change_groups %>% 
#   mutate(study_id = as.integer(study_id)) %>%
#   filter(Jaccard_base == "neg", S == "pos", SES_FRic == "pos", SES_FDiv == "0", SES_FEve == "0") %>%
#   left_join(study_meta)