In diazrenata/isds: Analysis of Mammalian Individual Size Distributions

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

Toy communities

set.seed(2019)
means <- c(40, 41, 100)
sds <- .1 * means

abund <- unlist(scads::sample_METE(s = 3, n = 200, nsamples = 1))

#abund <- sample(abund, size = 3, replace = F)

abund <- sort(abund, decreasing = T)

comm <- list() 

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

comm <- bind_rows(comm)

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

comm_plot_all

comm_plot <- ggplot(data = comm, aes(x = wgt, fill = species, group = species)) +
  geom_density(alpha =  .4) +
  theme_bw() +
  xlim(0, 200)

comm_plot

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)
}

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 <- unlist(scads::sample_METE(s = 6, n = 200, nsamples = 1))

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() +
  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

six_two_o_e <- list()
for(i in 1:nrow(six_two_overlap)) {
  six_two_o_e[[i]] <- rep(six_two_overlap$overlap[i], times = round(six_two_overlap$hm_abund[i]))
}
six_two_o_e <- data.frame(
  overlap = unlist(six_two_o_e))
six_two_hm <- ggplot(data = six_two_o_e, aes(x= overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()
six_two_hm

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() +
  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

six_six_o_e <- list()
for(i in 1:nrow(six_six_overlap)) {
  six_six_o_e[[i]] <- rep(six_six_overlap$overlap[i], times = round(six_six_overlap$hm_abund[i]))
}
six_six_o_e <- data.frame(
  overlap = unlist(six_six_o_e))
six_six_hm <- ggplot(data = six_six_o_e, aes(x= overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()
six_six_hm

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() +
  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
two_two_two_o_e <- list()
for(i in 1:nrow(two_two_two_overlap)) {
  two_two_two_o_e[[i]] <- rep(two_two_two_overlap$overlap[i], times = round(two_two_two_overlap$hm_abund[i]))
}
two_two_two_o_e <- data.frame(
  overlap = unlist(two_two_two_o_e))
two_two_hm <- ggplot(data = two_two_two_o_e, aes(x= overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()
two_two_hm

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() +
  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

six_one_o_e <- list()
for(i in 1:nrow(six_one_overlap)) {
  six_one_o_e[[i]] <- rep(six_one_overlap$overlap[i], times = round(six_one_overlap$hm_abund[i]))
}
six_one_o_e <- data.frame(
  overlap = unlist(six_one_o_e))
six_one_hm <- ggplot(data = six_one_o_e, aes(x= overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()
six_one_hm

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() +
  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

six_ov_o_e <- list()
for(i in 1:nrow(six_ov_overlap)) {
  six_ov_o_e[[i]] <- rep(six_ov_overlap$overlap[i], times = round(six_ov_overlap$hm_abund[i]))
}
six_ov_o_e <- data.frame(
  overlap = unlist(six_ov_o_e))
six_ov_hm <- ggplot(data = six_ov_o_e, aes(x= overlap)) +
  geom_histogram(boundary = 0) + geom_vline(xintercept = c(0, 1), color = "red") + xlim(-.15, 1.15) +
  theme_bw()
six_ov_hm