In diazrenata/neonbecs: What the Package Does (One Line, Title Case)

library(neonbecs)
# replicatebecs::download_raw_paper_data()
# replicatebecs::process_raw_data()
dat <- replicatebecs::load_paper_data()[[1]]

generate an ISD

dat_isd <- dat %>%
  replicatebecs::add_energy_sizeclass() %>%
  neonbecs::make_isd()

fit a gmm

dat_gmm <- fit_gmm(dat_isd)

dat_pdf <- get_pdf(dat_gmm)

integrate density and find peaks and pits

source(here::here("functions", "integrated_density.R"))

interval_size <- 0.001
min_size <- 0
max_size <- 8

integrated_density <- get_integrated_density(dat_gmm = dat_gmm, 
                                             interval_size = interval_size,
                                             min_size = min_size,
                                             max_size = max_size) %>%
  add_all_pit_boundaries(threshold = .05) 
integrated_density_plot_wb <- plot_integrated_density(integrated_density = integrated_density,
                                                      threshold_lines = FALSE,
                                                      pit_boundaries = TRUE)

print(integrated_density_plot_wb)

modes <- dplyr::filter(integrated_density, start_is_peak) %>%
  dplyr::select(start, by_max)
pits <- dplyr::filter(integrated_density, start_is_pit) %>%
  dplyr::select(start, by_max)

modes
pits

Look at species positions along the PDF...

intd_tojoin <- integrated_density %>%
  dplyr::rename(ln_size2 = start) %>%
  dplyr::mutate(ln_size2 = as.character(ln_size2))

dat_p <- dat_isd %>%
  dplyr::mutate(ln_size2 = as.character(signif(ln_size, digits = 3))) %>%
  dplyr::left_join(intd_tojoin, by = "ln_size2") %>%
  dplyr::select(-ln_size2)

species_summaries <- dat_p %>%
  dplyr::group_by(individual_species_ids) %>%
  dplyr::summarise(mean_size = mean(ln_size),
                   mean_p = mean(by_max),
                   sd_p = sd(by_max),
                   prop_in_trough = mean(is_in_trough),
                   n_ind = dplyr::n()) %>%
  dplyr::ungroup()
species_summaries <- species_summaries %>%
  dplyr::mutate(sd_p = tidyr::replace_na(data = species_summaries$sd_p, replace = 0))


individuals_plot <- ggplot2::ggplot(data = dat_p, ggplot2::aes(x = ln_size, y = by_max, color = individual_species_ids)) + 
  ggplot2::geom_jitter(size = .5) +
  ggplot2::geom_point(data = species_summaries, ggplot2::aes(x = mean_size, y = mean_p, color = individual_species_ids), size = 3) +
  ggplot2::geom_hline(yintercept = 0.05) +

  ggplot2::theme_bw()

individuals_plot

Look at the composition of the troughs

trough_sizes <- integrated_density %>%
  dplyr::filter(is_in_trough) %>%
  dplyr::select(start) %>%
  dplyr::mutate(ln_size = signif(start, digits = 4)) %>%
  dplyr::select(ln_size) %>%
  dplyr::distinct()

trough_ind <- dat_isd %>%
  dplyr::mutate(ln_size = signif(ln_size, digits = 4)) %>%
  dplyr::filter(ln_size %in% trough_sizes$ln_size) %>%
  dplyr::left_join(species_summaries, by = "individual_species_ids") %>%
  dplyr::select(-individual_sizes, -individual_energy, -size_class, -size_class_g, -ln_energy) %>%
  dplyr::arrange(dplyr::desc(individual_species_ids))

trough_plot <- ggplot2::ggplot(data = trough_ind, ggplot2::aes(x = prop_in_trough, y = ln_size/mean_size, color = individual_species_ids, size = n_ind)) +
  ggplot2::xlim(0, 1) + 
  ggplot2::ylim(.5, 1.5) +
  ggplot2::geom_hline(yintercept = 1) +
  ggplot2::geom_jitter() +
  ggplot2::theme_bw()

trough_plot

Mode structure

get_mode_dom <- function(threshold, dat_p) {

  N <- as.numeric(nrow(dat_p))
  S <- as.numeric(length(unique(dat_p$individual_species_ids)))

  modes <- dat_p %>%
    dplyr::filter(by_max > threshold)

  mode_richness <- as.numeric(length(unique(modes$individual_species_ids)))
  mode_abund <- as.numeric(nrow(modes))

  mode_dom <- (mode_richness / S) / (mode_abund / N)

  return(mode_dom)
}
mode_thresholds <- data.frame(
  threshold = seq(.01, 1, by = 0.01),
  mode_dom = vapply(seq(0.01, 1, by = 0.01), FUN = get_mode_dom, 
                    dat_p = dat_p, FUN.VALUE = .75)
)

thresh_plot <- ggplot2::ggplot(data = mode_thresholds, 
                               ggplot2::aes(x = threshold, y = mode_dom)) +
  ggplot2::geom_point() + 
  ggplot2::geom_hline(yintercept = 1) +
  ggplot2::theme_bw()

thresh_plot