In diazrenata/isds: Analysis of Mammalian Individual Size Distributions

library(isds)
library(ggplot2)
library(dplyr)

many_sims <- read.csv(here::here("analysis", "sims_nounif.csv"))

ln_units <- .2

one_sim <- many_sims %>%
  filter(source == "constrained", time_chunk == "eighties", sim == 1) %>%
  mutate(species = as.factor(species),
         ln_mass = log(wgt), 
         size_class = ln_units * (floor(ln_mass/ln_units)),
         size_class_g = exp(size_class))

one_sim_counts <- one_sim %>%
  group_by(species, size_class, size_class_g) %>%
  summarize(nind = dplyr::n()) %>%
  ungroup() %>%
  group_by(species) %>%
  mutate(total_ind_species = sum(nind)) %>%
  ungroup()  %>%
  mutate(ind_proportional = nind / total_ind_species)

isd_plot <- ggplot(data = one_sim_counts, aes(x= size_class, y = nind)) +
  geom_col() +
  theme_bw()
isd_plot

ssd_plot <- ggplot(data = one_sim_counts, aes(x = size_class, y = nind, fill = species)) +
  geom_col(alpha = .2, position = "dodge") +
  theme_bw() +
  scale_fill_viridis_d(end = .8)
ssd_plot

two_species <- one_sim_counts %>%
  filter(species %in% c(1, 2))

two_species_plot <- ggplot(one_sim_counts, aes(x = size_class, y = ind_proportional, color = species)) +
  geom_point(data = two_species) +
  geom_line(data = two_species) +
  theme_bw() +
  scale_color_viridis_d(end = .8, option = "plasma")

two_species_plot

all_species <- expand.grid(unique(as.numeric(one_sim_counts$species)), unique(as.numeric(one_sim_counts$species))) %>%
  rename(sp1 = Var1, sp2 = Var2) %>%
  filter(sp1 < sp2) %>%
  as.matrix()

expanded <- apply(all_species, MARGIN = 1, FUN = function(sp_vect, sim_counts) return(mutate(filter(sim_counts, as.numeric(species) %in% sp_vect), sp1 = sp_vect[1], sp2 = sp_vect[2])), sim_counts = one_sim_counts)

all_species <- bind_rows(expanded)

all_comb_plots <- ggplot(all_species, aes(x = size_class, y = ind_proportional, color = species)) +
  geom_point(data = all_species) +
  geom_line(data = all_species) +
  theme_bw() +
  scale_color_viridis_d(end = .8) + 
  facet_wrap(vars(sp1, sp2), scales = "free") +
  theme(strip.text = element_blank())

all_comb_plots

library(mclust)
all_comb <- expand.grid(1:9, 1:9)  %>%
 # filter(Var1 != Var2) %>%
  mutate(p = NA)

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

reference_isd <- filter(one_sim, species == sp1)

reference_density <- densityMclust(reference_isd$wgt, modelNames = "V")

id <- data.frame(wgt = seq(0, floor(1.25 * max(one_sim$wgt)), by = .1),
                density = NA)

id$density <- predict(reference_density, newdata = id$wgt)

id$density <- id$density / sum(id$density) 

id$wgt = round(id$wgt, digits = 1)

sp2_p <- filter(one_sim, species == sp2) %>%
  mutate(wgt = round(wgt, digits = 1)) %>%
  left_join(id, by = "wgt")

illust <- ggplot(data = id, aes(x = wgt, y = density)) +
  geom_point() +
  geom_point(data = sp2_p, aes(x = wgt, y = mean(id$density)), color = "red", alpha = .1) +
  geom_point(x = mean(sp2_p$wgt), y = mean(sp2_p$density), color = "green") +
  theme_bw() +
  ggtitle(mean(sp2_p$density))

print(illust)

all_comb$p[i] <- mean(log(sp2_p$density))
}

phist <- ggplot(data = all_comb, aes(x = p)) +
  geom_histogram() +
  theme_bw()

phist