# With the root.dir option below, # this vignette runs the R code in a temporary directory # so new files are written to temporary storage # and not the user's file space. knitr::opts_knit$set(root.dir = fs::dir_create(tempfile())) knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) skip <- identical(Sys.getenv("NOT_CRAN", unset = "false"), "false") || !requireNamespace("rjags", quietly = TRUE) || !requireNamespace("R2jags", quietly = TRUE) if (skip) { knitr::opts_chunk$set(eval = FALSE) } else { library(R2jags) } library(dplyr) library(targets) library(jagstargets)
The introductory vignette vignette caters to Bayesian data analysis workflows with few datasets to analyze. However, it is sometimes desirable to run one or more Bayesian models repeatedly across many simulated datasets. Examples:
This vignette focuses on (1). The goal of this particular example to simulate multiple datasets from the model below, analyze each dataset, and assess how often the estimated posterior intervals cover the true parameters from the prior predictive simulations. The quantile method by @cook2006 generalizes this concept, and simulation-based calibration [@talts2020] generalizes further. The interval-based technique featured in this vignette is not as robust as SBC, but it may be more expedient for large models because it does not require visual inspection of multiple histograms.
Consider a simple regression model with a continuous response y
with a covariate x
.
$$ \begin{aligned} y_i &\stackrel{\text{iid}}{\sim} \text{Normal}(\beta_1 + x_i \beta_2, 1) \ \beta_1, \beta_2 &\stackrel{\text{iid}}{\sim} \text{Normal}(0, 1) \end{aligned} $$
We write this model in a JAGS model file.
lines <- "model { for (i in 1:n) { y[i] ~ dnorm(beta[1] + x[i] * beta[2], 1) } for (i in 1:2) { beta[i] ~ dnorm(0, 1) } }" writeLines(lines, "model.jags")
Next, we define a pipeline to simulate multiple datasets and fit each dataset with the model. In our data-generating function, we put the true parameter values of each simulation in a special .join_data
list. jagstargets
will automatically join the elements of .join_data
to the correspondingly named variables in the summary output. This will make it super easy to check how often our posterior intervals capture the truth. As for scale, generate 20 datasets (5 batches with 4 replications each) and run the model on each of the 20 datasets.^[Internally, each batch is a dynamic branch target, and the number of replications determines the amount of work done within a branch. In the general case, batching is a way to find the right compromise between target-specific overhead and the horizontal scale of the pipeline.] By default, each of the 20 model runs computes 3 MCMC chains with 2000 MCMC iterations each (including burn-in) and you can adjust with the n.chains
and n.iter
arguments of tar_jags_rep_summary()
.
# Writes the _targets.R file shown in the next code chunk. library(targets) tar_script({ library(jagstargets) options(crayon.enabled = FALSE) tar_option_set(memory = "transient", garbage_collection = TRUE) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, "model.jags", data = generate_data(), parameters.to.save = "beta", batches = 5, # Number of branch targets. reps = 4, # Number of model reps per branch target. stdout = R.utils::nullfile(), stderr = R.utils::nullfile(), variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ) ) })
# _targets.R library(targets) library(jagstargets) options(crayon.enabled = FALSE) # Use computer memory more sparingly: tar_option_set(memory = "transient", garbage_collection = TRUE) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) # Elements of .join_data get joined on to the .join_data column # in the summary output next to the model parameters # with the same names. .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, "model.jags", data = generate_data(), parameters.to.save = "beta", batches = 5, # Number of branch targets. reps = 4, # Number of model reps per branch target. variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ) )
We now have a pipeline that runs the model 10 times: 5 batches (branch targets) with 4 replications per batch.
tar_visnetwork()
Run the computation with tar_make()
tar_make()
The result is an aggregated data frame of summary statistics, where the .rep
column distinguishes among individual replicates. We have the posterior intervals for beta
in columns q2.5
and q97.5
. And thanks to the .join_data
list we included in generate_data()
, our output has a .join_data
column with the true values of the parameters in our simulations.
tar_load(model)
model
Now, let's assess how often the estimated 95% posterior intervals capture the true values of beta
. If the model is implemented correctly, the coverage value below should be close to 95%. (Ordinarily, we would increase the number of batches and reps per batch and run batches in parallel computing.)
library(dplyr) model %>% group_by(variable) %>% dplyr::summarize(coverage = mean(q2.5 < .join_data & .join_data < q97.5))
For maximum reproducibility, we should express the coverage assessment as a custom function and a target in the pipeline.
# Writes the _targets.R file shown in the next code chunk. library(targets) tar_script({ library(jagstargets) options(crayon.enabled = FALSE) tar_option_set( packages = "dplyr", memory = "transient", garbage_collection = TRUE ) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) # Elements of .join_data get joined on to the .join_data column # in the summary output next to the model parameters # with the same names. .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, "model.jags", data = generate_data(), parameters.to.save = "beta", batches = 5, # Number of branch targets. reps = 4, # Number of model reps per branch target. stdout = R.utils::nullfile(), stderr = R.utils::nullfile(), variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ), tar_target( coverage, model %>% group_by(variable) %>% summarize( coverage = mean(q2.5 < .join_data & .join_data < q97.5), .groups = "drop" ) ) ) })
# _targets.R # packages needed to define the pipeline: library(targets) library(jagstargets) tar_option_set( packages = "dplyr", # packages needed to run the pipeline memory = "transient", # memory efficiency garbage_collection = TRUE # memory efficiency ) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) # Elements of .join_data get joined on to the .join_data column # in the summary output next to the model parameters # with the same names. .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, "model.jags", data = generate_data(), parameters.to.save = "beta", batches = 5, # Number of branch targets. reps = 4, # Number of model reps per branch target. variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ), tar_target( coverage, model %>% group_by(variable) %>% summarize( coverage = mean(q2.5 < .join_data & .join_data < q97.5), .groups = "drop" ) ) )
The new coverage
target should the only outdated target, and it should be connected to the upstream model
target.
tar_visnetwork()
When we run the pipeline, only the coverage assessment should run. That way, we skip all the expensive computation of simulating datasets and running MCMC multiple times.
tar_make()
tar_read(coverage)
tar_jags_rep_mcmc_summary()
and similar functions allow you to supply multiple jags models. If you do, each model will share the the same collection of datasets, and the .dataset_id
column of the model target output allows for custom analyses that compare different models against each other. Below, we add a new model2.jags
file to the jags_files
argument of tar_jags_rep_mcmc_summary()
. In the coverage summary below, we group by .name
to compute a coverage statistic for each model.
lines <- "model { for (i in 1:n) { y[i] ~ dnorm(beta[1] + x[i] * x[i] * beta[2], 1) # Regress on x^2, not x. } for (i in 1:2) { beta[i] ~ dnorm(0, 1) } }" writeLines(lines, "model2.jags")
# Writes the _targets.R file shown in the next code chunk. library(targets) tar_script({ library(jagstargets) options(crayon.enabled = FALSE) tar_option_set( packages = "dplyr", memory = "transient", garbage_collection = TRUE ) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) # Elements of .join_data get joined on to the .join_data column # in the summary output next to the model parameters # with the same names. .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, c("model.jags", "model2.jags"), # another model data = generate_data(), parameters.to.save = "beta", batches = 5, reps = 4, stdout = R.utils::nullfile(), stderr = R.utils::nullfile(), variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ), tar_target( coverage, model %>% group_by(.name) %>% summarize(coverage = mean(q2.5 < .join_data & .join_data < q97.5)) ) ) })
# _targets.R # packages needed to define the pipeline: library(targets) library(jagstargets) tar_option_set( packages = "dplyr", # packages needed to run the pipeline memory = "transient", # memory efficiency garbage_collection = TRUE # memory efficiency ) generate_data <- function(n = 10L) { beta <- stats::rnorm(n = 2, mean = 0, sd = 1) x <- seq(from = -1, to = 1, length.out = n) y <- stats::rnorm(n, beta[1] + x * beta[2], 1) # Elements of .join_data get joined on to the .join_data column # in the summary output next to the model parameters # with the same names. .join_data <- list(beta = beta) list(n = n, x = x, y = y, .join_data = .join_data) } list( tar_jags_rep_summary( model, c("model.jags", "model2.jags"), # another model data = generate_data(), parameters.to.save = "beta", batches = 5, reps = 4, variables = "beta", summaries = list( ~posterior::quantile2(.x, probs = c(0.025, 0.975)) ) ), tar_target( coverage, model %>% group_by(.name) %>% summarize(coverage = mean(q2.5 < .join_data & .join_data < q97.5)) ) )
In the graph below, notice how targets model_model1
and model_model2
are both connected to model_data
upstream. Downstream, model
is equivalent to dplyr::bind_rows(model_model1, model_model2)
, and it will have special columns .name
and .file
to distinguish among all the models.
tar_visnetwork()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.