Introduction to jagstargets

The jagstargets package makes it easy to run a single jags model and keep track of the results. R2jags fits the models, and targets manages the workflow and helps avoid unnecessary computation.

Consider the simple regression model below with response variable y and covariate x.

\[ \begin{aligned} y_i &\stackrel{\text{iid}}{\sim} \text{Normal}(x_i \beta, 1) \\ \beta &\sim \text{Normal}(0, 1) \end{aligned} \]

We write this model in the JAGS model file below.

lines <- "model {
  for (i in 1:n) {
    y[i] ~ dnorm(x[i] * beta, 1)
  }
  beta ~ dnorm(0, 1)
}"
writeLines(lines, "x.jags")

A typical workflow proceeds as follows:

  1. Prepare a list of input data to JAGS, including vector elements x and y.
  2. Fit the JAGS model using the list of input data.
  3. Use the fitted model object to compute posterior summaries and convergence diagnostics.
  4. Use the fitted model object to extract posterior draws of parameters and store them in a tidy data frame.
  5. If there are other models to compare, use the fitted model object to compute the deviance information criterion (DIC).

jagstargets encapsulates this workflow with the tar_jags() function. To use it in a targets pipeline, invoke it from the _targets.R script of the project.

# _targets.R
library(targets)
library(jagstargets)

generate_data <- function(n = 10) {
  true_beta <- stats::rnorm(n = 1, mean = 0, sd = 1)
  x <- seq(from = -1, to = 1, length.out = n)
  y <- stats::rnorm(n, x * true_beta, 1)
  out <- list(n = n, x = x, y = y)
}

# The _targets.R file ends with a list of target objects
# produced by jagstargets::tar_jags(), targets::tar_target(), or similar.
list(
  tar_jags(
    example,
    jags_files = "x.jags",
    parameters.to.save = "beta",
    data = generate_data()
  )
)

tar_jags() only defines the pipeline. It does not actually run JAGS, it declares the targets that will eventually run JAGS. The specific targets are as follows. Run tar_manifest() to show specific details about the targets declared.

tar_manifest()
#> # A tibble: 7 × 2
#>   name              command                                                     
#>   <chr>             <chr>                                                       
#> 1 example_data      "tar_jags_example_data()"                                   
#> 2 example_file_x    "\"x.jags\""                                                
#> 3 example_lines_x   "readLines(con = example_file_x)"                           
#> 4 example_mcmc_x    "jagstargets::tar_jags_run(jags_lines = example_lines_x, \n…
#> 5 example_summary_x "jagstargets::tar_jags_df(example_mcmc_x, data = example_da…
#> 6 example_dic_x     "jagstargets::tar_jags_df(fit = example_mcmc_x, data = exam…
#> 7 example_draws_x   "jagstargets::tar_jags_df(fit = example_mcmc_x, data = exam…

Each target is responsible for a piece of the workflow.

The suffix _x comes from the base name of the model file, in this case x.jags. If you supply multiple model files to the jags_files argument, all the models share the same dataset, and the suffixes distinguish among the various targets.

The targets depend on one another: for example, example_mcmc_x takes example_data as input. targets can visualize the dependency relationships in a dependency graph, which is helpful for understanding the pipeline and troubleshooting issues.

tar_visnetwork(targets_only = TRUE)

Run the computation with tar_make().

tar_make()
#> • start target example_data
#> • built target example_data [0.002 seconds]
#> • start target example_file_x
#> • built target example_file_x [0.001 seconds]
#> • start target example_lines_x
#> • built target example_lines_x [0.001 seconds]
#> • start target example_mcmc_x
#> • built target example_mcmc_x [0.294 seconds]
#> • start target example_summary_x
#> • built target example_summary_x [0.102 seconds]
#> • start target example_dic_x
#> • built target example_dic_x [0.004 seconds]
#> • start target example_draws_x
#> • built target example_draws_x [0.002 seconds]
#> • end pipeline [0.664 seconds]

The output lives in a special folder called _targets/ and you can retrieve it with functions tar_load() and tar_read() (from targets).

tar_read(example_summary_x)
#> # A tibble: 2 × 11
#>   variable   mean median    sd   mad     q5    q95  rhat ess_b…¹ ess_t…² .join…³
#>   <chr>     <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>   <dbl>   <dbl>   <dbl>
#> 1 beta     -0.160 -0.161 0.438 0.434 -0.884  0.558  1.00   2451.   2689.  0.0275
#> 2 deviance 29.9   29.4   1.12  0.487 29.1   32.1    1.00   3081.   2930. NA     
#> # … with abbreviated variable names ¹​ess_bulk, ²​ess_tail, ³​.join_data

At this point, all our results are up to date because their dependencies did not change.

tar_make()
#> ✔ skip target example_data
#> ✔ skip target example_file_x
#> ✔ skip target example_lines_x
#> ✔ skip target example_mcmc_x
#> ✔ skip target example_summary_x
#> ✔ skip target example_dic_x
#> ✔ skip target example_draws_x
#> ✔ skip pipeline [0.073 seconds]

But if we change the underlying code or data, some of the targets will no longer be valid, and they will rerun during the next tar_make(). Below, we change the jags model file, so the MCMC reruns while the data is skipped. This behavior saves time and enhances reproducibility.

write(" ", file = "x.jags", append = TRUE)
tar_outdated()
#> [1] "example_summary_x" "example_dic_x"     "example_file_x"   
#> [4] "example_draws_x"   "example_mcmc_x"    "example_lines_x"
tar_visnetwork(targets_only = TRUE)
tar_make()
#> ✔ skip target example_data
#> • start target example_file_x
#> • built target example_file_x [0.002 seconds]
#> • start target example_lines_x
#> • built target example_lines_x [0.001 seconds]
#> • start target example_mcmc_x
#> • built target example_mcmc_x [0.08 seconds]
#> • start target example_summary_x
#> • built target example_summary_x [0.047 seconds]
#> • start target example_dic_x
#> • built target example_dic_x [0.004 seconds]
#> • start target example_draws_x
#> • built target example_draws_x [0.001 seconds]
#> • end pipeline [0.346 seconds]

At this point, we can add more targets and custom functions for additional post-processing. See below for a custom summary target (which is equivalent to customizing the summaries argument of tar_jags().)

# _targets.R
library(targets)
library(jagstargets)

generate_data <- function(n = 10) {
  true_beta <- stats::rnorm(n = 1, mean = 0, sd = 1)
  x <- seq(from = -1, to = 1, length.out = n)
  y <- stats::rnorm(n, x * true_beta, 1)
  out <- list(n = n, x = x, y = y)
}

list(
  tar_jags(
    example,
    jags_files = "x.jags",
    parameters.to.save = "beta",
    data = generate_data()
  ),
  tar_target(
    custom_summary,
    posterior::summarize_draws(
      dplyr::select(example_draws_x, -starts_with(".")),
      ~posterior::quantile2(.x, probs = c(0.25, 0.75))
    )
  )
)

In the graph, our new custom_summary target should be connected to the upstream example target, and only custom_summary should be out of date.

tar_visnetwork(targets_only = TRUE)

In the next tar_make(), we skip the expensive MCMC and just run the custom summary.

tar_make()
#> ✔ skip target example_data
#> ✔ skip target example_file_x
#> ✔ skip target example_lines_x
#> ✔ skip target example_mcmc_x
#> ✔ skip target example_summary_x
#> ✔ skip target example_draws_x
#> ✔ skip target example_dic_x
#> • start target custom_summary
#> • built target custom_summary [0.017 seconds]
#> • end pipeline [0.209 seconds]
tar_read(custom_summary)
#> # A tibble: 2 × 3
#>   variable    q25    q75
#>   <chr>     <dbl>  <dbl>
#> 1 beta     -0.447  0.136
#> 2 deviance 29.2   30.1

Multiple models

tar_jags() and related functions allow you to supply multiple models to jags_files. If you do, each model will run on the same dataset. Consider a new model, y.jags.

lines <- "model {
  for (i in 1:n) {
    y[i] ~ dnorm(x[i] * x[i] * beta, 1) # Regress on x^2 instead of x.
  }
  beta ~ dnorm(0, 1)
}"
writeLines(lines, "y.jags")

Below, we add y.jags to the jags_files argument of tar_jags().

# _targets.R
library(targets)
library(jagstargets)

generate_data <- function(n = 10) {
  true_beta <- stats::rnorm(n = 1, mean = 0, sd = 1)
  x <- seq(from = -1, to = 1, length.out = n)
  y <- stats::rnorm(n, x * true_beta, 1)
  out <- list(n = n, x = x, y = y)
}

list(
  tar_jags(
    example,
    jags_files = c("x.jags", "y.jags"),
    parameters.to.save = "beta",
    data = generate_data()
  ),
  tar_target(
    custom_summary,
    posterior::summarize_draws(
      dplyr::select(example_draws_x, -starts_with(".")),
      ~posterior::quantile2(.x, probs = c(0.25, 0.75))
    )
  )
)

In the graph below, notice how the *_x targets and *_y targets are both connected to example_data upstream.

tar_visnetwork(targets_only = TRUE)

More information

For more on targets, please visit the reference website https://docs.ropensci.org/targets/ or the user manual https://books.ropensci.org/targets/. The manual walks though advanced features of targets such as high-performance computing and cloud storage support.