In diazrenata/isds: Analysis of Mammalian Individual Size Distributions

knitr::opts_chunk$set(echo = TRUE)
library(isds)
library(dplyr)
library(ggplot2)

sample_community <- function(means, sds, abunds) {

  comm <- list() 

  for(i in 1:length(means)) {
    comm[[i]] <- data.frame(
      species = i,
      wgt = rnorm(mean = means[i], sd = sds[i], n = abunds[i])
    )
  }

  comm <- bind_rows(comm)

  return(comm)
}


get_overlap <- function(community_df) {
  nspp <- length(unique(community_df$species))
  all_pairs <- expand.grid(1:nspp, 1:nspp) %>%
    as.data.frame() %>%
    filter(Var1 < Var2) %>%
    mutate(overlap = 0)


  for(i in 1:nrow(all_pairs)) {
    sp1 = all_pairs[i, 1]
    sp2 = all_pairs[i, 2]

    sp1_density <- density(filter(community_df, species == sp1)$wgt, from = 0, to = 600, n = 1024)$y

    sp1_density <- sp1_density / sum(sp1_density)

    sp2_density <- density(filter(community_df, species == sp2)$wgt, from = 0, to = 600, n = 1024)$y
    sp2_density <- sp2_density / sum(sp2_density)

    min_density <- data.frame(
      sp1 = sp1_density,
      sp2 = sp2_density
    )

    min_density <- min_density %>%
      mutate(index = row_number()) %>%
      group_by(index) %>%
      mutate(min_den = min(sp1, sp2)) %>%
      ungroup()

    overlap <- sum(min_density$min_den)


    all_pairs$overlap[i] <- overlap

  }

  abund_df <- community_df %>%
    group_by(species) %>%
    summarize(abund = dplyr::n()) %>%
    ungroup()

  all_pairs_w <- all_pairs %>%
    left_join(abund_df, by = c("Var1" = "species")) %>%
    rename(abund1 = abund) %>%
    left_join(abund_df, by = c("Var2" = "species")) %>%
    rename(abund2 = abund) %>%
    mutate(index = row_number()) %>%
    group_by(index) %>%
    mutate(hm_abund = psych::harmonic.mean(c(abund1, abund2)),
           total_abund = sum(abund1, abund2),
           prop_abund = sum(abund1, abund2) / sum(abund_df$abund)) %>%
    ungroup()
  return(all_pairs_w)
}


scenario_names <- c("six_one", "six_ov", "six_two", "six_six", "two_two_two")

summ_plots <- list()
overlaps <- list()

Scenarios:

Six species, two overlapped groups

set.seed(2)
means_62 <- c(20, 21, 22, 80, 81, 82) 
sds_62 <- .2 * means_62
an_sad <- c(50, 50, 50, 50, 50, 50)

an_sad[ which(an_sad == 1)] <- 2

six_two <- sample_community(means_62, sds_62, an_sad)


six_two_plot_all <- ggplot(data = six_two, aes(x = wgt)) +
  geom_density() +
  theme_bw() +
  xlim(0, 200)

six_two_plot_all

six_two_plot <- ggplot(data = six_two, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) +  theme_bw() +
   theme(legend.position = "none") +
  xlim(0, 200)

six_two_plot

six_two_overlap <- get_overlap(six_two)

six_two_o <- ggplot(data = six_two_overlap, aes(x = overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()

six_two_o

summ_plots$six_two <- list(six_two_plot_all, six_two_plot, six_two_o)
overlaps$six_two <- six_two_overlap

Six nonoverlapping species

means_66 <- c(20,80, 130, 170, 250, 350) 
sds_66 <- .05 * means_66
six_six <- sample_community(means_66, sds_66, an_sad)


six_six_plot_all <- ggplot(data = six_six, aes(x = wgt)) +
  geom_density() +
  theme_bw() +
  xlim(0, 500)

six_six_plot_all

six_six_plot <- ggplot(data = six_six, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) +  theme_bw() +
   theme(legend.position = "none") +
  xlim(0, 500)

six_six_plot


six_six_overlap <- get_overlap(six_six)

six_six_o <- ggplot(data = six_six_overlap, aes(x = overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()

six_six_o
summ_plots$six_six <- list(six_six_plot_all, six_six_plot, six_six_o)

overlaps$six_six <- six_six_overlap

Two groups of two overlapping with two that span the overlap

means_222 <- c(31, 32, 45, 60, 85, 86)
sds_222 <- .15 * means_222

two_two_two <- sample_community(means_222, sds_222, an_sad)

two_two_two_plot_all <- ggplot(data = two_two_two, aes(x = wgt)) +
  geom_density() +
  theme_bw() +
  xlim(0, 200)

two_two_two_plot_all

two_two_two_plot <- ggplot(data = two_two_two, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) + theme_bw() +
   theme(legend.position = "none") +
  xlim(0, 200)

two_two_two_plot


two_two_two_overlap <- get_overlap(two_two_two)

two_two_two_o <- ggplot(data = two_two_two_overlap, aes(x = overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()

two_two_two_o
summ_plots$two_two_two <- list(two_two_two_plot_all, two_two_two_plot, two_two_two_o)
overlaps$two_two_two <- two_two_two_overlap

Six species, one group

means_61 <- c(45, 46, 47, 48, 49, 50)
sds_61 <- .15 * means_61

six_one <- sample_community(means_61, sds_61, an_sad)

six_one_plot_all <- ggplot(data = six_one, aes(x = wgt)) +
  geom_density() +
  theme_bw() +
  xlim(0, 200)

six_one_plot_all

six_one_plot <- ggplot(data = six_one, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) + theme_bw() +
   theme(legend.position = "none") +
  xlim(0, 200)

six_one_plot


six_one_overlap <- get_overlap(six_one)

six_one_o <- ggplot(data = six_one_overlap, aes(x = overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()

six_one_o
summ_plots$six_one <- list(six_one_plot_all, six_one_plot, six_one_o)
overlaps$six_one <- six_one_overlap

Six overlapping

means_6 <- c(30, 40, 50, 60, 70, 80)
sds_6 <- .15 * means_6

six_ov <- sample_community(means_6, sds_6, an_sad)

six_ov_plot_all <- ggplot(data = six_ov, aes(x = wgt)) +
  geom_density() +
  theme_bw() +
  xlim(0, 200)

six_ov_plot_all

six_ov_plot <- ggplot(data = six_ov, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) +
  theme_bw() +
   theme(legend.position = "none") +
  xlim(0, 200)

six_ov_plot


six_ov_overlap <- get_overlap(six_ov)

six_ov_o <- ggplot(data = six_ov_overlap, aes(x = overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()

six_ov_o
summ_plots$six_ov <- list(six_ov_plot_all, six_ov_plot, six_ov_o)
overlaps$six_ov <- six_ov_overlap

for(i in 1:length(summ_plots)) {
  gridExtra::grid.arrange(grobs = summ_plots[[i]], nrow = 1)
}

overlaps <- bind_rows(overlaps, .id = "source")

overlaps <- overlaps %>%
  mutate(abs_overlap = abs(overlap - .5))

overlap_hists <- ggplot(data = overlaps, aes(x = overlap)) +
  geom_histogram() +
  facet_wrap(vars(source), scales = "free_y") +
  theme_bw()
overlap_hists

overlaps <- group_by(overlaps, source) %>%
  mutate(med_overlap = median(abs_overlap),
         quarter_overlap = quantile(abs_overlap, probs = .25)) %>%
  ungroup()

overlap_hists1 <- ggplot(data = overlaps, aes(x = abs_overlap)) +
  geom_histogram()  + 
  geom_point(aes(x =quarter_overlap, y = 5), color = "pink") +
  facet_wrap(vars(source), scales = "free_y") +
  theme_bw()
overlap_hists1