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

knitr::opts_chunk$set(echo = FALSE)
library(bbstrends)
library(dplyr)
library(ggplot2)
theme_set(theme_bw())
source(here::here("analysis", "demo_fxns.R"))

Load specific route

The New Hartford route goes up and down Riverton Road and was started in 1994. It feels pretty auspicious. It is route 102, region 18.

ibd <- readRDS("~/Documents/GitHub/BBSsize/analysis/isd_data/ibd_isd_bbs_rtrg_102_18.Rds")

#ibd <- read.csv("~/Documents/GitHub/BBSsize/analysis/isd_data/isd_bbs_rtrg_7_4.csv")

ibd <- ibd %>%
  mutate(ind_size = abs(ind_size)) %>%
  mutate(log_size = log(ind_size))

Here is how species richness, abundance, biomass, and energy have changed over those years:

# 
# sv <- ibd %>%
#   group_by(year) %>%
#   summarize(abundance = dplyr::n(),
#             biomass = sum(ind_size),
#             energy = sum(ind_b)) %>%
#   ungroup() %>%
#   mutate(         mean_energy = energy/abundance)
# 

sv <- ibd %>%
  group_by(year) %>%
  summarize(abundance = dplyr::n(),
            energy = sum(ind_b),
            biomass = sum(ind_size)) %>%
  ungroup() %>%
  mutate(mean_energy = energy / abundance,
         mean_size = biomass/abundance)


sv_long <- sv %>%
  tidyr::pivot_longer(-year, names_to = "currency", values_to = "value")

ggplot(filter(sv_long), aes(x = year, y = value, color = currency)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_color_viridis_d(end = .8) +
  theme(legend.position = "none") +
  facet_wrap(vars(currency), scales = "free_y", ncol = 1)

ISD

Using the logarithm of mass.

ggplot(ibd, aes(log_size, group = year)) +
  geom_density(alpha = .1)

years <- as.list(unique(ibd$year))
min_size <- 1 * min(ibd$log_size)
max_size <- 1 * max(ibd$log_size)
npoints <- 500

annual_isds <- lapply(years, FUN = function(thisyear, ibd) return(filter(ibd, year == thisyear)), ibd = ibd)

annual_isd_smooths <- lapply(annual_isds, FUN = function(isd, min_size, max_size, npoints) return(data.frame(year = isd$year[1],
                                                                                log_size = seq(min_size, max_size, length.out = npoints),
                                                                                density = make_kde(isd$log_size, min_size = min_size, max_size = max_size, n = npoints))), min_size = min_size, max_size = max_size, npoints = npoints)


isd_smooths <- bind_rows(annual_isd_smooths) %>%
  mutate(yearf = as.factor(year),
         fiveyrs = ceiling((year - 1993.5) / 5)) %>%
  mutate(fiveyrsf = as.factor(fiveyrs))

ggplot(isd_smooths, aes(log_size, density, group = yearf, color = fiveyrsf)) +
  geom_line()


ggplot(filter(isd_smooths, fiveyrs %in% c(1, 5)), aes(log_size, density, group = yearf, color = fiveyrsf)) +
  geom_line()

library(mgcv)
library(gratia)
source(here::here("R", "gams_fxns_generalized.R"))

isd_gams_all <- gam(density ~ s(log_size, k = 50) + s(log_size, by = fiveyrsf, k = 20), data = isd_smooths, family = "tw")

#isd_gams_all <- gam(density ~ s(log_size) + s(log_size, by = fiveyrs), data = isd_smooths, family = "tw")

isd_fit <- isd_smooths %>% 
  add_fitted(isd_gams_all, exclude= c("s(yearf)"))

ggplot(isd_fit, aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3) +
  geom_line(aes(y = density, group= year), alpha = .5) +
  facet_wrap(vars(fiveyrsf))


ggplot(filter(isd_fit, fiveyrsf %in% c(1, 5)), aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3)

ggplot(filter(isd_fit, fiveyrsf %in% c(1, 5)), aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3) +
  geom_line(aes(y = density, group= year), alpha = .5) +
  facet_wrap(vars(fiveyrsf))

isd_pdat <- make_pdat(isd_smooths, np = 2000, comparison_variable = "fiveyrsf")
isd_pred <- get_predicted_vals(isd_gams_all, isd_pdat, exclude= c("s(yearf)"))
plot_fitted_pred(isd_pred, "fiveyrsf")
plot_fitted_pred(filter(isd_pred, fiveyrsf %in% c(1, 5)), "fiveyrsf")

isd_diff <- get_exclosure_diff(isd_gams_all, isd_pdat, "fiveyrsf", 1, 5, exclude = c("s(yearf)"))
plot_exclosure_diff(isd_diff)

Reshuffling years...

year_change <- data.frame(year = unlist(years),
                          newyear = sample(unlist((years)), size = length(years), replace = F))

isd_smooths_rs <- bind_rows(annual_isd_smooths) %>%
  left_join(year_change) %>%
  select(-year) %>%
  rename(year = newyear) %>%
  mutate(yearf = as.factor(year),
         fiveyrs = ceiling((year - 1993.5) / 5)) %>%
  mutate(fiveyrsf = as.factor(fiveyrs))

ggplot(isd_smooths_rs, aes(log_size, density, group = yearf, color = fiveyrsf)) +
  geom_line()


ggplot(filter(isd_smooths_rs, fiveyrs %in% c(1, 5)), aes(log_size, density, group = yearf, color = fiveyrsf)) +
  geom_line()

isd_rs_gams_all <- gam(density ~ s(log_size, k = 50) + s(log_size, by = fiveyrsf, k = 20), data = isd_smooths_rs, family = "tw")

#isd_rs_gams_all <- gam(density ~ s(log_size) + s(log_size, by = fiveyrs), data = isd_smooths_rs, family = "tw")

isd_rs_fit <- isd_smooths_rs %>% 
  add_fitted(isd_rs_gams_all, exclude= c("s(yearf)"))

ggplot(isd_rs_fit, aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3) +
  geom_line(aes(y = density, group= year), alpha = .5) +
  facet_wrap(vars(fiveyrsf))


ggplot(filter(isd_rs_fit, fiveyrsf %in% c(1, 5)), aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3)

ggplot(filter(isd_rs_fit, fiveyrsf %in% c(1, 5)), aes(log_size, .value, color = fiveyrsf)) +
  geom_line(linetype = 3) +
  geom_line(aes(y = density, group= year), alpha = .5) +
  facet_wrap(vars(fiveyrsf))

isd_rs_pdat <- make_pdat(isd_smooths_rs, np = 2000, comparison_variable = "fiveyrsf")
isd_rs_pred <- get_predicted_vals(isd_rs_gams_all, isd_rs_pdat, exclude= c("s(yearf)"))
plot_fitted_pred(isd_rs_pred, "fiveyrsf")
plot_fitted_pred(filter(isd_rs_pred, fiveyrsf %in% c(1, 5)), "fiveyrsf")

isd_rs_diff <- get_exclosure_diff(isd_rs_gams_all, isd_rs_pdat, "fiveyrsf", 1, 5, exclude = c("s(yearf)"))
plot_exclosure_diff(isd_rs_diff)