The stantargets
package makes it easy to run a single Stan model and keep track of the results. cmdstanr
fits the models, and targets
manages the workflow and helps avoid unnecessary computation.
First, write a Stan model file.
lines <- "data {
int <lower = 1> n;
vector[n] x;
vector[n] y;
real true_beta;
}
parameters {
real beta;
}
model {
y ~ normal(x * beta, 1);
beta ~ normal(0, 1);
}"
writeLines(lines, "x.stan")
A typical workflow proceeds as follows:
- Prepare a list of input data to Stan, including vector elements
x
andy
. - Fit the Stan model using the list of input data.
- Use the fitted model object to compute posterior summaries and convergence diagnostics.
- Use the fitted model object to extract posterior draws of parameters and store them in a tidy data frame.
- Use the fitted model to compute Hamiltonian Monte Carlo (HMC) diagnostics.
stantargets
expresses this workflow using the tar_stan_mcmc()
function. To use it in a targets
pipeline, invoke it from the _targets.R
script of the project.
# _targets.R
library(targets)
library(stantargets)
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)
list(n = n, x = x, y = y, true_beta = true_beta)
}
# The _targets.R file ends with a list of target objects
# produced by stantargets::tar_stan_mcmc(), targets::tar_target(), or similar.
list(
tar_stan_mcmc(
example,
"x.stan",
generate_data(),
stdout = R.utils::nullfile(),
stderr = R.utils::nullfile()
)
)
Above, tar_stan_mcmc(example, ...)
only defines the pipeline. It does not actually run Stan, it declares the targets that will eventually run Stan. Run tar_manifest()
to show specific details about the targets.
tar_manifest()
#> # A tibble: 6 × 2
#> name command
#> <chr> <chr>
#> 1 example_data "generate_data()"
#> 2 example_file_x "\"x.stan\""
#> 3 example_mcmc_x "stantargets::tar_stan_mcmc_run(stan_file = example_fil…
#> 4 example_diagnostics_x "tibble::as_tibble(posterior::as_draws_df(example_mcmc_…
#> 5 example_summary_x "stantargets::tar_stan_summary_join_data(summaries = ti…
#> 6 example_draws_x "tibble::as_tibble(posterior::as_draws_df(example_mcmc_…
Each target listed above is responsible for a piece of the workflow.
-
example_file_x
: Reproducibly track changes to the Stan model file. -
example_data
: Run the code you supplied to thedata
argument oftar_stan_mcmc()
and return a dataset compatible with Stan. -
example_mcmc_x
: Run the MCMC and return an object of classCmdStanMCMC
. -
example_draws_X
: Return a friendlytibble
of the posterior draws fromexample
. Uses the$draws()
method. Suppress withdraws = FALSE
intar_stan_mcmc()
. -
example_summaries_x
: Return a friendlytibble
of the posterior summaries fromexample
. Uses the$summary()
method. Suppress withsummary = FALSE
intar_stan_mcmc()
. -
example_diagnostics_x
: Return a friendlytibble
of the sampler diagnostics fromexample
. Uses the$sampler_diagnostics()
method. Suppress withdiagnostics = FALSE
intar_stan_mcmc()
.
The suffix _x
comes from the base name of the model file, in this case x.stan
. If you supply multiple model files to the stan_files
argument, all the models share the same dataset, and the suffixes distinguish among the various targets.
targets
produces a graph to show the dependency relationships among the targets. Below, the MCMC depends on the model file and the data, and the draws, summary, and diagnostics depend on the MCMC.
tar_visnetwork(targets_only = TRUE)
Run the computation with tar_make()
.
tar_make()
#> • start target example_data
#> • built target example_data [0.001 seconds]
#> • start target example_file_x
#> • built target example_file_x [0 seconds]
#> • start target example_mcmc_x
#> • built target example_mcmc_x [8.572 seconds]
#> • start target example_diagnostics_x
#> • built target example_diagnostics_x [0.002 seconds]
#> • start target example_summary_x
#> • built target example_summary_x [0.053 seconds]
#> • start target example_draws_x
#> • built target example_draws_x [0.003 seconds]
#> • end pipeline [9.245 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 lp__ -5.86 -5.60 0.700 0.310 -7.31 -5.37 1.00 1799. 2381. NA
#> 2 beta -0.154 -0.153 0.442 0.444 -0.890 0.564 1.00 1356. 1857. 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_mcmc_x
#> ✔ skip target example_diagnostics_x
#> ✔ skip target example_summary_x
#> ✔ skip target example_draws_x
#> ✔ skip pipeline [0.067 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 Stan model file, so the MCMC reruns while the data is skipped. This behavior saves time and enhances reproducibility.
write(" ", file = "x.stan", append = TRUE)
tar_outdated()
#> [1] "example_summary_x" "example_diagnostics_x" "example_file_x"
#> [4] "example_draws_x" "example_mcmc_x"
tar_visnetwork(targets_only = TRUE)
tar_make()
#> ✔ skip target example_data
#> • start target example_file_x
#> • built target example_file_x [0.001 seconds]
#> • start target example_mcmc_x
#> • built target example_mcmc_x [8.564 seconds]
#> • start target example_diagnostics_x
#> • built target example_diagnostics_x [0.002 seconds]
#> • start target example_summary_x
#> • built target example_summary_x [0.051 seconds]
#> • start target example_draws_x
#> • built target example_draws_x [0.002 seconds]
#> • end pipeline [9.193 seconds]
#>
At this point, we can add more targets and custom functions for additional post-processing.
# _targets.R
library(targets)
library(stantargets)
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)
list(n = n, x = x, y = y, true_beta = true_beta)
}
list(
tar_stan_mcmc(
example,
"x.stan",
generate_data(),
stdout = R.utils::nullfile(),
stderr = R.utils::nullfile()
),
tar_stan_summary(
custom_summary,
fit = example_mcmc_x,
summaries = list(~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_mcmc_x
#> ✔ skip target example_diagnostics_x
#> ✔ skip target example_summary_x
#> • start target custom_summary
#> • built target custom_summary [0.028 seconds]
#> ✔ skip target example_draws_x
#> • end pipeline [0.217 seconds]
#>
tar_read(custom_summary)
#> # A tibble: 2 × 4
#> variable q25 q75 .join_data
#> <chr> <dbl> <dbl> <dbl>
#> 1 lp__ -6.02 -5.42 NA
#> 2 beta -0.447 0.153 NA
Multiple models
tar_stan_mcmc()
and related functions allow you to supply multiple models to stan_files
. If you do, each model will run on the same dataset. Consider a new model y.stan
.
lines <- "data {
int <lower = 1> n;
vector[n] x;
vector[n] y;
real true_beta;
}
parameters {
real beta;
}
model {
y ~ normal(x * x * beta, 1); // Regress on x^2 instead of x.
beta ~ normal(0, 1);
}"
writeLines(lines, "y.stan")
To include this y.stan
, we add it to the stan_files
argument of tar_stan_mcmc()
.
#> [1] FALSE
# _targets.R
library(targets)
library(stantargets)
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)
list(n = n, x = x, y = y, true_beta = true_beta)
}
list(
tar_stan_mcmc(
example,
c("x.stan", "y.stan"), # another model
generate_data(),
stdout = R.utils::nullfile(),
stderr = R.utils::nullfile()
),
tar_stan_summary(
custom_summary,
fit = example_mcmc_x,
summaries = list(~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)
Generated quantities
It is possible to use the CmdStanMCMC
object from one run to simulate generated quantities downstream. For example, the tar_stan_gq_rep_summaries()
function takes a single CmdStanMCMC
object, produces multiple replications of generated quantities from multiple models, and aggregates the summaries from each. The following pipeline uses this technique to repeatedly draw from the posterior predictive distribution.
lines <- "data {
int <lower = 1> n;
vector[n] x;
vector[n] y;
}
parameters {
real beta;
}
model {
y ~ normal(x * beta, 1);
beta ~ normal(0, 1);
}
generated quantities {
real y_rep[n] = normal_rng(x * beta, 1); // posterior predictive draws
}"
writeLines(lines, "gen.stan")
# _targets.R
library(targets)
library(stantargets)
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)
list(n = n, x = x, y = y, true_beta = true_beta)
}
list(
tar_stan_mcmc(
example,
"x.stan",
generate_data(),
stdout = R.utils::nullfile(),
stderr = R.utils::nullfile()
),
tar_stan_gq_rep_summary(
postpred,
stan_files = "gen.stan",
fitted_params = example_mcmc_x, # one CmdStanFit object
data = generate_data(), # Function runs once per rep.
batches = 2, # 2 dynamic branches
reps = 5, # 5 replications per branch
quiet = TRUE,
stdout = R.utils::nullfile(),
stderr = R.utils::nullfile()
)
)
Since we have defined many objects in the pipeline, it is extra important to check the graph to be sure everything is connected.
tar_visnetwork(targets_only = TRUE)
Then, we run the computation. The original MCMC is already up to date, so we only run the targets relevant to the generated quantities.
tar_make()
#> • start target postpred_batch
#> • built target postpred_batch [0.001 seconds]
#> • start target postpred_file_gen
#> • built target postpred_file_gen [7.942 seconds]
#> ✔ skip target example_data
#> ✔ skip target example_file_x
#> • start branch postpred_data_f1c15717
#> • built branch postpred_data_f1c15717 [0.011 seconds]
#> • start branch postpred_data_e945e52d
#> • built branch postpred_data_e945e52d [0.008 seconds]
#> • built pattern postpred_data
#> ✔ skip target example_mcmc_x
#> • start branch postpred_gen_0a30dca5
#> • built branch postpred_gen_0a30dca5 [3.548 seconds]
#> • start branch postpred_gen_6a425504
#> • built branch postpred_gen_6a425504 [3.337 seconds]
#> • built pattern postpred_gen
#> ✔ skip target example_diagnostics_x
#> ✔ skip target example_summary_x
#> ✔ skip target example_draws_x
#> • start target postpred
#> • built target postpred [0 seconds]
#> • end pipeline [15.939 seconds]
#>
Finally, we read the summaries of posterior predictive samples.
tar_read(postpred)
#> # A tibble: 100 × 16
#> varia…¹ mean median sd mad q5 q95 rhat ess_b…² ess_t…³ .join…⁴
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 y_rep[… 0.138 0.169 1.10 1.10 -1.70 1.92 1.00 3074. 3049. NA
#> 2 y_rep[… 0.111 0.127 1.04 1.01 -1.62 1.83 1.00 3317. 3649. NA
#> 3 y_rep[… 0.111 0.125 1.05 1.04 -1.65 1.82 1.00 3882. 3974. NA
#> 4 y_rep[… 0.101 0.0917 1.01 0.981 -1.54 1.84 1.00 4151. 3595. NA
#> 5 y_rep[… 0.0668 0.0201 1.02 1.01 -1.64 1.71 1.00 4111. 3816. NA
#> 6 y_rep[… -0.0459 -0.0367 1.00 0.994 -1.73 1.57 1.00 4339. 3531. NA
#> 7 y_rep[… -0.0421 -0.0397 0.997 0.984 -1.65 1.59 1.00 3743. 3719. NA
#> 8 y_rep[… -0.0813 -0.0746 1.02 1.00 -1.77 1.58 1.00 3101. 3531. NA
#> 9 y_rep[… -0.145 -0.145 1.06 1.07 -1.86 1.62 1.00 3430. 3336. NA
#> 10 y_rep[… -0.127 -0.121 1.13 1.11 -2.02 1.70 1.00 2856. 3227. NA
#> # … with 90 more rows, 5 more variables: .rep <chr>, .dataset_id <chr>,
#> # .seed <int>, .file <chr>, .name <chr>, and abbreviated variable names
#> # ¹variable, ²ess_bulk, ³ess_tail, ⁴.join_data
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.