In jasonmtroos/sbcrs: R package for implementing simulation based calibration using rank statistics

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

library(sbcrs)
library(rstan)

rstan::rstan_options("auto_write" = TRUE)
get_compiled_stan_model <- function(filename) {
  m <- NULL
  stan_file_loc <- here::here('inst', 'stan', filename)
  if (file.exists(stan_file_loc)) {
    m <- stan_model(file = stan_file_loc, save_dso = TRUE)
  }
  if (is.null(m)) {
    m <- stan_model(file = system.file('stan', filename, package = 'sbcrs'))
  }
  m
}
sbc_rstan_model <- get_compiled_stan_model('rstan_sbc_example.stan')
sbc_original_model <- get_compiled_stan_model('rstan_sbc_example_original.stan')

We will compare the ranks calculated using the SBC package against those calculated by the rstan::sbc() function.

The stan code for this model is based on the help text for rstan::sbc. The following model is given as an example:

cat(readr::read_file(system.file('stan', 'rstan_sbc_example.stan', package = 'sbcrs')))

Compile this Stan model:

sbc_rstan_model <- stan_model(file = system.file('stan', 'rstan_sbc_example.stan', package = 'sbcrs'))

Calibration involves generating data and parameters, and sampling from a Stan model many times. To speed up this process, take advantage of all of your machine's cores.

doParallel::registerDoParallel(cores = parallel::detectCores())
options(mc.cores = parallel::detectCores())

if (!identical(Sys.getenv("NOT_CRAN"), "true") && !interactive()) {
  doParallel::registerDoParallel(cores = 2)
  options(mc.cores = 2)
}

Calibrate using rstan::sbc().

calibration_data <- list(N = 10, a = 2, b = 2)
rstan_sbc <- rstan::sbc(sbc_rstan_model, data = calibration_data, 256)
plot(rstan_sbc, binwidth = 1, thin = 50)

The Stan code used in the above example has been modified from the original. In the modified version, y is generated in the transformed data block of the Stan file. The original model would have looked like this:

cat(readr::read_file(system.file('stan', 'rstan_sbc_example_original.stan', package = 'sbcrs')))

Compile this Stan model:

sbc_original_model <- stan_model(file = system.file('stan', 'rstan_sbc_example_original.stan', package = 'sbcrs'))

Create an SBC object that corresponds with the original model.

sbc <- SBC$new(
  data = function(seed) {
    calibration_data
  },
  params = function(seed, data) {
    set.seed(seed + 1e6)
    list(pi = rbeta(1, data$a, data$b))
  },
  modeled_data = function(seed, data, params) {
    set.seed(seed + 2e6)
    list(y = rbinom(1, data$N, params$pi))
  },
  sampling = function(seed, data, params, modeled_data, iters) {
    sampling(sbc_original_model, data = c(data, modeled_data), seed = seed,
             chains = 1, iter = 2 * iters, warmup = iters)
  })

sbc$calibrate(256, 50)
sbc$plot()

Assess whether the distributions of recovered ranks are similar

library(purrr)
x <- 
  map(sbc$calibrations, 'ranks') %>%
  flatten() %>%
  unlist() %>%
  unname()

y <- 
  rstan_sbc$ranks %>%
  map(~.x[seq(1, 1000, by = 20), ]) %>%
  map(~sum(.x)) %>%
  unlist()

qqplot(x, y, xlab = 'SBC', ylab = 'rstan')
abline(a = 0, b = 1, col = 'red')

wilcox.test(x, y)

jasonmtroos/sbcrs documentation built on Nov. 4, 2019, 2:20 p.m.