vignettes/article/estimate_trophic_position_two_source_model_ar.R

In trps: Bayesian Trophic Position Models using 'stan'

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, message=FALSE-----------------------------------------------------
{
library(bayesplot)
library(brms)
library(dplyr)
library(ggplot2)
library(ggdist)
library(grid)
library(tidybayes)
library(tidyr)
library(trps)
library(viridis)
}

## -----------------------------------------------------------------------------
consumer_iso

## -----------------------------------------------------------------------------
baseline_1_iso

## -----------------------------------------------------------------------------
baseline_2_iso

## -----------------------------------------------------------------------------
con_tsar <- consumer_iso %>% 
  arrange(ecoregion, common_name) %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    id = row_number()
  ) %>% 
  ungroup() %>% 
  dplyr::select(id, common_name:d15n)

## -----------------------------------------------------------------------------
b1_tsar <- baseline_1_iso %>% 
  arrange(ecoregion, common_name) %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    id = row_number()
  ) %>% 
  ungroup() %>% 
  dplyr::select(id, ecoregion:d15n_b1)

## -----------------------------------------------------------------------------
b2_tsar <- baseline_2_iso %>% 
  arrange(ecoregion, common_name) %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    id = row_number()
  ) %>% 
  ungroup() %>% 
  dplyr::select(id, ecoregion:d15n_b2)

## -----------------------------------------------------------------------------
combined_iso_tsar <- con_tsar %>% 
  left_join(b1_tsar, by = c("id", "ecoregion")) %>% 
  left_join(b2_tsar, by = c("id", "ecoregion")) 

## -----------------------------------------------------------------------------
combined_iso_tsar_1 <- combined_iso_tsar %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    c1 = mean(d13c_b1, na.rm = TRUE),
    n1 = mean(d15n_b1, na.rm = TRUE),
    c2 = mean(d13c_b2, na.rm = TRUE),
    n2 = mean(d15n_b2, na.rm = TRUE),
    l1 = 2
  ) %>% 
  ungroup() %>% 
  add_alpha()

## -----------------------------------------------------------------------------
combined_iso_tsar_1

## ----eval = FALSE-------------------------------------------------------------
# model_output_tsar <- brm(
#   formula = two_source_model_ar(),
#   prior = two_source_priors_ar(),
#   stanvars = two_source_priors_params_ar(),
#   data = combined_iso_tsar_1,
#   family = gaussian(),
#   chains = 2,
#   iter = 4000,
#   warmup = 1000,
#   cores = 4,
#   seed = 4,
#   control = list(adapt_delta = 0.95)
# )

## -----------------------------------------------------------------------------
model_output_tsar

## -----------------------------------------------------------------------------
plot(model_output_tsar)

## -----------------------------------------------------------------------------
pp_check(model_output_tsar, resp = "alpha")

## -----------------------------------------------------------------------------
pp_check(model_output_tsar, resp = "d15n")

## -----------------------------------------------------------------------------
model_output_tsar

## -----------------------------------------------------------------------------
get_variables(model_output_tsar)

## -----------------------------------------------------------------------------
post_draws <- model_output_tsar %>% 
  gather_draws(b_alpha_ar_Intercept, b_d15n_tp_Intercept) %>% 
  mutate(
    ecoregion = "Embayment",
    common_name = "Lake Trout",
    .variable = case_when(
      .variable %in% "b_d15n_tp_Intercept" ~ "tp",
      .variable %in% "b_alpha_ar_Intercept" ~ "ar"
    )
  ) %>% 
  dplyr::select(common_name, ecoregion, .chain:.value)

## -----------------------------------------------------------------------------
post_draws

## ----message=FALSE------------------------------------------------------------
medians_ci <- model_output_tsar %>%
        gather_draws(b_alpha_ar_Intercept, 
                     b_d15n_tp_Intercept) %>%
        median_qi() %>% 
  mutate(
    ecoregion = "Embayment", 
    common_name = "Lake Trout", 
    .variable = case_when(
      .variable %in% "b_d15n_tp_Intercept" ~ "tp",
      .variable %in% "b_alpha_ar_Intercept" ~ "ar"
    )
  ) %>% 
  mutate_if(is.numeric, round, digits = 2)

## -----------------------------------------------------------------------------
medians_ci

## -----------------------------------------------------------------------------
ggplot(data = post_draws, aes(x = .value)) + 
  geom_density() + 
  facet_wrap(~ .variable, scale = "free") +
  theme_bw(base_size = 15) +
  theme(
    panel.grid = element_blank(), 
    strip.background = element_blank(),
  ) + 
  labs(
    x = "P(Estimate | X)", 
    y = "Density"
  )

## -----------------------------------------------------------------------------
ggplot(data = post_draws, aes(y = .value, 
                              x = common_name)) + 
  stat_pointinterval() + 
   facet_wrap(~ .variable, scale = "free") +
  theme_bw(base_size = 15) +
  theme(
    panel.grid = element_blank(), 
    strip.background = element_blank(),
  ) + 
  labs(
    x = "P(Estimate | X)", 
    y = "Density"
  )