In poissonconsulting/smbr: Facilitates Bayesian Analysis using STAN

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "tools/README-"
)

smbr

Introduction

smbr (pronounced simber) is an R package to facilitate analyses using STAN. It is part of the embr family of packages.

Demonstration

library(bauw)
library(ggplot2)
library(magrittr)
library(embr)
library(smbr)

# define model in Stan language
model <- model("
  data {
      int nAnnual;
      int nObs;
      int Annual[nObs];
      int Pairs[nObs];
      real Year[nObs];
  }
  parameters {
      vector[nAnnual] bAnnual;
      real log_sAnnual;
      real alpha;
      real beta1;
      real beta2;
      real beta3;
  }
  transformed parameters {
    real sAnnual;
    sAnnual = exp(log_sAnnual);
  }
  model {
      vector[nObs] ePairs;

      log_sAnnual ~ normal(0, 10);
      bAnnual ~ normal(0, sAnnual);

      alpha ~ normal(0, 10);
      beta1 ~ normal(0, 10);
      beta2 ~ normal(0, 10);
      beta3 ~ normal(0, 10);

      for (i in 1:nObs) {
        ePairs[i] = exp(alpha + beta1 * Year[i] + beta2 * Year[i]^2 +
                      beta3 * Year[i]^3 + bAnnual[Annual[i]]);
      }
      target += poisson_lpmf(Pairs | ePairs);
  }")

# add R code to calculate derived parameters
model %<>% update_model(new_expr = "
  for (i in 1:length(Pairs)) {
    prediction[i] <- exp(alpha + beta1 * Year[i] + beta2 * Year[i]^2 +
                       beta3 * Year[i]^3 + bAnnual[Annual[i]])
  }
")

# define data types and center year
model %<>% update_model(
  select_data = list(
    "Pairs" = integer(), "Year*" = integer(),
    Annual = factor()
  ),
  derived = "sAnnual",
  random_effects = list(bAnnual = "Annual")
)

data <- bauw::peregrine
data$Annual <- factor(data$Year)

set.seed(42)

# analyse
analysis <- analyse(model, data = data, seed = 3L, glance = FALSE)

# coefficient table
coef(analysis, simplify = TRUE)

# trace plots
plot(analysis)

# make predictions by varying year with other predictors including the random effect of Annual held constant
year <- predict(analysis, new_data = "Year")

# plot those predictions
ggplot(data = year, aes(x = Year, y = estimate)) +
  geom_point(data = bauw::peregrine, aes(y = Pairs)) +
  geom_line() +
  geom_line(aes(y = lower), linetype = "dotted") +
  geom_line(aes(y = upper), linetype = "dotted") +
  expand_limits(y = 0)

Installation

# install.packages("devtools")
devtools::install_github("poissonconsulting/smbr")

Citation

citation(package = "smbr")

Contribution

Please report any issues.

Pull requests are always welcome.

Code of Conduct

Please note that the smbr project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.

poissonconsulting/smbr documentation built on Jan. 14, 2024, 5:59 a.m.