Eight-schools: a real cmdstanr workflow • bayesgrove

Fit the centered eight-schools model with cmdstanr, review its sampler diagnostics, and rerun a non-centered version on a branch. The chunks use eval = FALSE with output captured from a real run, so the vignette builds without CmdStan on CI.

Setup

library(bayesgrove)

handle <- bg_init(
  path = file.path(tempdir(), "bg-eight-schools"),
  project_name = "Eight Schools"
)

# Activate the default review loop and the cmdstanr built-in executors.
bg_use_workflow_packs(handle, "bayesgrove.default_bayesian")
bg_use_cmdstanr(handle)

Build the graph

# Eight-schools data (Rubin, 1981; Gelman et al., 2013).
schools_data <- list(
  J = 8L,
  y = c(28, 8, -3, 7, -1, 1, 18, 12),
  sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
)

# Node params are data, never code. To preserve R types that JSON serialization
# would mangle (e.g. integers), attach the data object via the content-addressed
# store with bg_set_node_data().
n_data <- bg_add_node(
  handle,
  kind = "stan_data",
  label = "Schools data"
)
bg_set_node_data(handle, n_data, schools_data)

The centered parameterization of the hierarchical model is a standard example of funnel geometry: as tau approaches zero, the conditional distribution of theta contracts sharply, and HMC produces divergent transitions (Betancourt and Girolami, 2015). The package ships the example programs so the code below runs as pasted:

stan_file <- system.file(
  "stan/eight_schools_centered.stan",
  package = "bayesgrove",
  mustWork = TRUE
)

n_fit <- bg_add_node(
  handle,
  kind = "cmdstanr_fit",
  label = "Centered",
  inputs = n_data,
  params = list(
    stan_file = stan_file,
    chains = 4,
    iter_warmup = 1000,
    iter_sampling = 1000,
    seed = 42
  )
)

Run the model and inspect diagnostics

run <- bg_run(handle, targets = n_fit)

#> Starting run "run_179e87ae" with 1 node to execute.
#> Running node "node_9f8f0f45"...
#> Running node "node_55624d87"...
#> Running MCMC with 4 sequential chains...
#> All 4 chains finished successfully.
#> Warning: 81 of 4000 (2.0%) transitions ended with a divergence.
#> Warning: 2 of 4 chains had an E-BFMI less than 0.3.

The executor stores these diagnostics as a typed hmc_diagnostics summary:

summaries <- bg_read_summaries(handle)
summaries[[1]]$severity
#> [1] "error"
summaries[[1]]$metrics[c("divergences", "e_bfmi", "max_rhat")]
#> $divergences
#> [1] 81
#>
#> $e_bfmi
#> [1] 0.1574
#>
#> $max_rhat
#> [1] 1.0401

Because the summary carries severity = "error", the default pack raises a blocking review_computation_validity obligation and suggests two actions: record the review decision, or branch and modify the fit.

guide <- bg_next_actions(handle)
vapply(guide$obligations, `[[`, character(1), "kind")
#>                  obl_ef958c36
#> "review_computation_validity"
vapply(guide$actions, `[[`, character(1), "kind")
#>        act_2260dedb        act_663128f1
#>   "record_decision" "branch_and_modify"

Record a review

The suggested record_decision action carries the decision kind and the ids of the summaries under review in its payload. Recording a computation_review decision with those ids resolves the obligation; a decision of another kind, or one that does not reference the summaries, does not.

review <- Filter(
  function(a) identical(a$payload$decision_type, "computation_review"),
  guide$actions
)[[1]]

bg_record_decision(
  handle,
  scope = "project",
  prompt = "Is the centered fit acceptable for downstream use?",
  choice = "reject_and_reparametrize",
  rationale = paste(
    "81 of 4000 transitions diverged and E-BFMI is low: the sampler cannot",
    "explore the funnel between tau and theta. Switch to the non-centered",
    "parameterization."
  ),
  kind = "computation_review",
  metadata = review$payload[c("node_ids", "summary_ids")]
)

vapply(bg_next_actions(handle)$obligations, `[[`, character(1), "kind")
#> named character(0)

Branch to the non-centered parameterization

Create a branch from the centered fit, retaining the original fit and its diagnostics. Set an inferential goal, as required for each active branch, and select the non-centered program. bg_update_node() merges params, so the chains, iteration, and seed settings carry over.

branch <- bg_branch(handle, n_fit, label = "Non-centered")

bg_set_goal(
  project = handle,
  branch_id = branch$branch_id,
  kind = "latent_inference",
  label = "School effects, non-centered",
  rationale = "Same model; reparameterized so HMC can explore the posterior."
)

stan_file_nc <- system.file(
  "stan/eight_schools_noncentered.stan",
  package = "bayesgrove",
  mustWork = TRUE
)

bg_update_node(
  handle,
  branch$root_node_id,
  params = list(stan_file = stan_file_nc)
)

Run and review the non-centered fit

run2 <- bg_run(handle, targets = branch$root_node_id)

#> Starting run "run_a26583a8" with 1 node to execute.
#> Running node "node_4c430039"...
#> Running MCMC with 4 sequential chains...
#> All 4 chains finished successfully.

The non-centered fit has no sampler warnings. Both fits remain in the record, with their diagnostic summaries.

for (s in bg_read_summaries(handle, include_stale = FALSE)) {
  cat(s$node_id, "-", s$severity, "- divergences:", s$metrics$divergences, "\n")
}
#> node_55624d87 - error - divergences: 81
#> node_4c430039 - ok - divergences: 0

One obligation remains. Decisions are scope-specific (see vignette("concepts")): the project-scope review recorded above does not cover the new branch. Record a branch-scope review of the centered fit’s error summary in its lineage.

branch_guide <- bg_next_actions(
  handle,
  scope = "branch",
  branch_id = branch$branch_id
)
vapply(branch_guide$obligations, `[[`, character(1), "kind")
#>                  obl_e02b096b
#> "review_computation_validity"

review_nc <- Filter(
  function(a) identical(a$payload$decision_type, "computation_review"),
  branch_guide$actions
)[[1]]

bg_record_decision(
  handle,
  scope = branch$branch_id,
  prompt = "Is the centered fit acceptable as this branch's baseline?",
  choice = "superseded_by_non_centered",
  rationale = paste(
    "The centered fit's divergences are the reason this branch exists;",
    "the non-centered fit replaces it."
  ),
  kind = "computation_review",
  metadata = review_nc$payload[c("node_ids", "summary_ids")]
)

bg_next_actions(handle)$obligations
#> named list()

No model comparison is needed here: the two fits target the same posterior, and the reparameterization is a computational repair rather than a change of model. For distinct candidate models, the default pack requires a comparison decision and an accept-or-reject disposition for each branch. See vignette("case-study-roaches").

Export a report

print(bg_read_graph(handle))
#> ── Execution Graph ──
#>
#> └── Schools data <stan_data> [new]
#>     ├── Centered <cmdstanr_fit> [new]
#>     └── Non-centered <cmdstanr_fit> [new]

bg_export_report(handle, format = "md", path = "eight-schools-report.md")

The report contains the graph with both fits, the review decisions at each scope, and the artifact index. Diagnostic summaries remain in the project.