inst/doc/trps.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

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

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

## -----------------------------------------------------------------------------
combined_iso_os <- con_os %>% 
  left_join(b1_os, by = c("id", "ecoregion"))

## -----------------------------------------------------------------------------
combined_iso_os_1 <- combined_iso_os %>% 
  group_by(ecoregion, common_name) %>% 
  mutate(
    c1 = mean(d13c_b1, na.rm = TRUE),
    n1 = mean(d15n_b1, na.rm = TRUE),
    l1 = 2
  ) %>% 
  ungroup()

## -----------------------------------------------------------------------------
combined_iso_os_1

## -----------------------------------------------------------------------------
model_output_os <- brm(
  formula = one_source_model(),
  prior = one_source_priors(),
  stanvars = one_source_priors_params(),
  data = combined_iso_os_1,
  family = gaussian(),
  chains = 2,
  iter = 4000,
  warmup = 1000,
  cores = 4,
  seed = 4,
  control = list(adapt_delta = 0.95)
)

## -----------------------------------------------------------------------------
model_output_os

## -----------------------------------------------------------------------------
plot(model_output_os)

## -----------------------------------------------------------------------------
pp_check(model_output_os)

## -----------------------------------------------------------------------------
model_output_os

## -----------------------------------------------------------------------------
get_variables(model_output_os)

## -----------------------------------------------------------------------------
post_draws <- model_output_os %>% 
  gather_draws(b_tp_Intercept) %>% 
  mutate(
    ecoregion = "Embayment",
    common_name = "Lake Trout",
    .variable = "tp"
  ) %>% 
  dplyr::select(common_name, ecoregion, .chain:.value)

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

## ----message=FALSE------------------------------------------------------------
medians_ci <- model_output_os %>%
        spread_draws(b_tp_Intercept) %>%
        median_qi() %>% 
        rename(
          tp = b_tp_Intercept
        ) %>% 
  mutate(
    ecoregion = "Embayment", 
    common_name = "Lake Trout"
  ) %>% 
  mutate_if(is.numeric, round, digits = 2) %>% 
  dplyr::select(ecoregion, common_name, tp:.interval)

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

## -----------------------------------------------------------------------------
ggplot(data = post_draws, aes(x = .value)) + 
  geom_density() + 
  theme_bw(base_size = 15) +
  theme(
    panel.grid = element_blank()
  ) + 
  labs(
    x = "P(Trophic Position | X)", 
    y = "Density"
  )

## -----------------------------------------------------------------------------
ggplot(data = post_draws, aes(y = .value, 
                              x = common_name)) + 
  stat_pointinterval() + 
  theme_bw(base_size = 15) +
  theme(
    panel.grid = element_blank()
  ) + 
  labs(
    x = "P(Trophic Position | X)", 
    y = "Density"
  )