In bsaul/elktoeChemistry: Data and code for Elktoe Chemistry project

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Examining Diversity in Elements within and between layers, individuals and sites

library(elktoeChemistry)
library(dplyr)
library(entropart)
library(ggplot2)
analysis_data  <- readRDS(params$inputs$analysis)

site_colors <- 
  c("LiTN 1" = "#99d8c9",
    "LiTN 2" = "#41ae76",
    "LiTN 3" = "#005824",
    "Tuck 1" = "#fdae6b",
    "Tuck 2" = "#f16913",
    "Tuck 3" = "#8c2d04",
    "Baseline" = "#525252")

dt <- 
  analysis_data(
    contrast = "all",
    test_data_FUN = identity,
    transect_opts  = list(
      .layers    = c("ipx", "ncr", "psm", "pio", "opx")
    )
  ) %>%
  select(signal, species, element, data) %>%
  tidyr::unnest(cols = data) %>%
  select(
    signal, species, final_status, river, site, 
    id, obs, transect, layer, annuli, element, value) %>%
  group_by(
    signal, species, final_status, river, site,
    id, layer, annuli, element,
  ) %>%
  summarise(
    value = sum(value)
  )

adt <- 
  dt %>% 
  ungroup() %>%
  mutate(
    layer_annuli = if_else(is.na(annuli),
                     as.character(layer),
                     paste(layer, annuli)),
    layer_annuli = factor(layer_annuli, 
                    levels = c("ipx", paste("ncr", LETTERS[1:13]),
                               "psm", "pio", "opx")),
  ) %>%
  select(
    signal, species, site, final_status, id, layer, annuli,
    layer_annuli, element, value
  )

get_diversity_data <- function(df){
  df %>%
  mutate(
    dp = purrr::map(
      .x = mc,
      .f = ~ DivProfile(
        q = seq(0, 2, by = 0.1), 
        MC = .x,
        Biased = FALSE)
    ),
    beta = purrr::map(
      dp,
      ~ tibble(
        q     = .x[["Order"]],
        value = .x[["TotalBetaDiversity"]]
      )
    )
  ) 
}

Across Nacral layers within individuals

ncrdt <-
  adt %>%
  filter(species == "Arav", layer == "ncr", annuli <= "C") %>%
  select(-final_status, -layer, -layer_annuli) %>%
  group_by(
    signal, species, site, id
  ) %>%
  filter(any(annuli > "A")) %>%
  group_nest() %>%
  mutate(
    data = purrr::map(
      .x = data,
      .f = ~ tidyr::pivot_wider(
        .x, names_from = "annuli", values_from = "value")
    ),
    mc = purrr::map(
      .x = data,
      .f = ~ MetaCommunity(.x[ , -1])
    ),
  ) %>%
  get_diversity_data()

pdt <- 
  ncrdt %>%
  select(signal, species, site, id, beta) %>%
  tidyr::unnest(beta)

ggplot(
  data = pdt,
  aes(x = q, y = value, id = id,
      color = site)
) +
  geom_line() +
  scale_color_manual(values = site_colors) +
  facet_grid(signal ~ .)

ncr_summary <- 
  ncrdt %>%
  select(
    signal, species, site, id, mc
  ) %>%
  group_by(signal, species, site) %>%
  group_nest() %>%
  mutate(
    mc = purrr::map(
      .x = data, 
      .f = ~ MergeMC(.x$mc, Weights = sapply(.x$mc, function(x) (x$N)))
    )
  ) %>%
  get_diversity_data()

pdt <- 
  ncr_summary %>%
  select(signal, species, site, beta) %>%
  tidyr::unnest(beta)

ggplot(
  data = pdt,
  aes(x = q, y = value, id = site,
      color = site)
) +
  geom_line() +
  scale_color_manual(values = site_colors) +
  facet_grid(signal ~ .)

Within Layers Between Sites

layer_dt <- 
  adt %>%
  filter(is.na(annuli) | annuli <= "C") %>%
  select(-final_status, -layer, -annuli) %>%
  group_by(signal, species, site, layer_annuli) %>%
  group_nest() %>%
  filter(species == "Arav") %>%
  mutate(
    data = purrr::map(
      .x = data,
      .f = ~ tidyr::pivot_wider(
        .x, names_from = "id", values_from = "value")
    ),
    mc = purrr::map(
      .x = data,
      .f = ~ MetaCommunity(.x[ , -1])
    )
  ) %>%
  get_diversity_data()

pdt <- 
  layer_dt %>%
  filter(signal == "base") %>%
  select(species, site, layer_annuli, beta) %>%
  tidyr::unnest(cols = beta)

ggplot(
  data = pdt,
  aes(x = q, y = value, group = site, color = site)
) +
  geom_line() + 
  scale_color_manual(
    values = site_colors
  ) + 
  scale_y_continuous(
    "Effective number of individuals"
  ) + 
  scale_x_continuous(
    "Order of diversity"
  ) + 
  facet_grid(species ~ layer_annuli, scales = "free_y") +
  theme_bw()

layer_dt_2 <- 
  layer_dt %>%
  # select(-final_status, -layer, -annuli) %>%
  group_by(signal, species, layer_annuli) %>%
  group_nest() %>%
  mutate(
    mc = purrr::map(
      .x = data, 
      .f = ~ MergeMC(.x$mc, Weights = sapply(.x$mc, function(x) (x$N)))
    )
  ) %>%
  get_diversity_data()

pdt <- 
  layer_dt_2 %>%
  filter(signal == "base") %>%
  select(species, layer_annuli, beta) %>%
  tidyr::unnest(cols = beta)

ggplot(
  data = pdt,
  aes(x = q, y = value)
) +
  geom_line() + 
  # scale_color_manual(
  #   values = site_colors
  # ) + 
  scale_y_continuous(
    "Effective number of sites"
  ) + 
  scale_x_continuous(
    "Order of diversity"
  ) + 
  facet_grid(species ~ layer_annuli, scales = "free_y") +
  theme_bw()

Within Layers by Status

layer_dt <- 
  adt %>%
  filter(is.na(annuli) | annuli <= "C") %>%
  filter(site != "Baseline", !is.na(final_status)) %>%
  select(-site, -layer, -annuli) %>%
  group_by(signal, final_status, species, layer_annuli) %>%
  group_nest() %>%
  filter(species == "Arav") %>%
  mutate(
    data = purrr::map(
      .x = data,
      .f = ~ tidyr::pivot_wider(
        .x, names_from = "id", values_from = "value")
    ),
    mc = purrr::map(
      .x = data,
      .f = ~ MetaCommunity(.x[ , -1])
    )
  ) %>%
  get_diversity_data()

pdt <- 
  layer_dt %>%
  filter(signal == "base") %>%
  select(species, final_status, layer_annuli, beta) %>%
  tidyr::unnest(cols = beta)

ggplot(
  data = pdt,
  aes(x = q, y = value, group = final_status, color = final_status)
) +
  geom_line() + 
  # scale_color_manual(
  #   values = site_colors
  # ) + 
  scale_y_continuous(
    "Effective number of individuals"
  ) + 
  scale_x_continuous(
    "Order of diversity"
  ) + 
  facet_grid(species ~ layer_annuli, scales = "free_y") +
  theme_bw()