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

knitr::opts_chunk$set(echo = FALSE)
#knitr::opts_chunk$set(fig.dim = c(5,3))

library(dplyr)
library(gratia)
library(ggplot2)
load_mgcv()

ts <- read.csv(here::here("analysis", "from_stories", "results", "ts_w_rescaled_e_100bbs.csv"))

set.seed(1977)


max_to_plot = 16

ones_to_plot <- unique(ts$site_name)[sample.int(n = length(unique(ts$site_name)), size = max_to_plot, replace = F)]

ts_long <- ts %>%
  select(year, abundance, scaled_energy, site_name) %>%
  tidyr::pivot_longer(c(abundance, scaled_energy), names_to = "currency", values_to = "value")

ggplot(filter(ts_long, site_name %in% ones_to_plot), aes(year, value, color= currency)) +
  geom_line() +
  theme_bw() +
  ggtitle("Actual values")+
  scale_color_viridis_d(end = .8) +
  facet_wrap(vars(site_name), scales = "free")

joint_samples <- read.csv(here::here("analysis", "from_stories",  "results", "joint_samples_5_100bbs.csv"))

ggplot(filter(joint_samples, identifier %in% ones_to_plot), aes(year, mean, color = currency, fill = currency)) +
  geom_line(aes(year, mean)) +
  geom_ribbon(aes(year, ymin = lower, ymax = upper),  alpha = .25) +
  ggtitle(paste0("Fits for energy and abundance, k = ", joint_samples$k[1])) +
  theme_bw() +
  scale_color_viridis_d(end = .8) +
  scale_fill_viridis_d(end = .8) + 
  facet_wrap(vars(identifier), scales = "free")

both_change <- read.csv(here::here("analysis", "from_stories",  "results", "both_change_5_100bbs.csv"))

both_summary <- read.csv(here::here("analysis", "from_stories",  "results", "both_summary_5_100bbs.csv"))

ggplot(both_summary, aes(x = mean_net_proportional)) +
  geom_histogram(binwidth = .1) +
  theme_bw() +
  facet_wrap(vars(currency), scales = "free") +
  geom_vline(xintercept = 0)

both_summary <- both_summary %>%
  group_by_all() %>%
  mutate(over_zero = ((lower_net_proportional * upper_net_proportional) < 0),
         increasing = all((lower_net_proportional > 0), upper_net_proportional > 0),
         decreasing = all((lower_net_proportional < 0), upper_net_proportional < 0),
         doubling = all(lower_net_proportional > 0, upper_net_proportional > 0, mean_net_proportional >= 1),
         halving = all(lower_net_proportional < 0, upper_net_proportional < 0, mean_net_proportional <= -.5)) %>%
  ungroup()

ggplot(filter(both_summary, !over_zero), aes(mean_net_proportional)) +
  geom_histogram(binwidth = .1) +
  facet_wrap(vars(currency)) +
  theme_bw()

both_summary %>%
  group_by(currency) %>%
  summarize(ntotal = dplyr::n(),
            noverzero = sum(over_zero),
            nincreasing = sum(increasing),
            ndecreasing = sum(decreasing),
            ndouble = sum(doubling),
            nhalf = sum(halving))

both_summary %>%
  select(identifier, currency, mean_net_proportional, over_zero, increasing, decreasing) %>%
  tidyr::pivot_longer(c(over_zero, increasing, decreasing), names_to = "description", values_to = "tf") %>%
  filter(tf) %>%
  select(-tf) %>%
  group_by(currency, description) %>%
  summarize(mean_shift = mean(mean_net_proportional))

both_summary_wide <- both_summary %>%
  tidyr::pivot_wider(id_cols = c(identifier, k), names_from = currency, values_from = c(mean_net_proportional, lower_net_proportional, upper_net_proportional))

ggplot(both_summary_wide, aes(x = mean_net_proportional_abundance, y = mean_net_proportional_scaled_energy, color = identifier)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_net_proportional_scaled_energy,
                    ymax = upper_net_proportional_scaled_energy)) +
  geom_errorbarh(aes(xmin = lower_net_proportional_abundance,
                     xmax = upper_net_proportional_abundance)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) + 
  geom_abline(slope = 1, intercept = 0) +
  theme(legend.position="none")



ggplot(both_summary_wide, aes(x = mean_net_proportional_abundance, y = mean_net_proportional_scaled_energy, color = identifier)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_net_proportional_scaled_energy,
                    ymax = upper_net_proportional_scaled_energy)) +
  geom_errorbarh(aes(xmin = lower_net_proportional_abundance,
                     xmax = upper_net_proportional_abundance)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) + 
  geom_abline(slope = 1, intercept = 0) +
  theme(legend.position="none") +
  xlim(-1, 1.5) +
  ylim(-1,1.5)

e_v_n_resid <- both_summary_wide %>%
  mutate(resid = mean_net_proportional_abundance - mean_net_proportional_scaled_energy)

ggplot(e_v_n_resid, aes(resid)) +
  geom_histogram(binwidth = .1) +
  xlab("abundance minus energy")


ggplot(e_v_n_resid, aes(mean_net_proportional_abundance, resid)) +
  geom_point() +
  ylab("abundance minus energy")