vignettes/article/estimate_trophic_position_two_source_model.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_ts <- consumer_iso %>% 
  arrange(ecoregion, common_name) %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    id = row_number()
  ) %>% 
  ungroup() %>% 
  dplyr::select(id, common_name:d15n)

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

## -----------------------------------------------------------------------------
b2_ts <- 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_ts <- con_ts %>% 
  left_join(b1_ts, by = c("id", "ecoregion")) %>% 
  left_join(b2_ts, by = c("id", "ecoregion")) 

## -----------------------------------------------------------------------------
combined_iso_ts_1 <- combined_iso_ts %>% 
  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()

## -----------------------------------------------------------------------------
combined_iso_ts_1

## -----------------------------------------------------------------------------
model_output_ts <- brm(
  formula = two_source_model(),
  prior = two_source_priors(),
  stanvars = two_source_priors_params(),
  data = combined_iso_ts_1,
  family = gaussian(),
  chains = 2,
  iter = 4000,
  warmup = 1000,
  cores = 4,
  seed = 4,
  control = list(adapt_delta = 0.95)
)

## -----------------------------------------------------------------------------
model_output_ts

## -----------------------------------------------------------------------------
plot(model_output_ts)

## -----------------------------------------------------------------------------
pp_check(model_output_ts, resp = "d13c")

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

## -----------------------------------------------------------------------------
model_output_ts

## -----------------------------------------------------------------------------
get_variables(model_output_ts)

## -----------------------------------------------------------------------------
post_draws <- model_output_ts %>% 
  gather_draws(b_d13c_alpha_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_d13c_alpha_Intercept" ~ "alpha"
    )
  ) %>% 
  dplyr::select(common_name, ecoregion, .chain:.value)

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

## ----message=FALSE------------------------------------------------------------
medians_ci <- model_output_ts %>%
        gather_draws(b_d13c_alpha_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_d13c_alpha_Intercept" ~ "alpha"
    )
  ) %>% 
  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"
  )