write posterior sample vignette

[n8thangreen]

---
title: "Using CEdecisiontree with Bayesian Posterior Samples (Stan)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using CEdecisiontree with Bayesian Posterior Samples (Stan)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

While CEdecisiontree provides internal functions for sampling from standard distributions (e.g., Beta, Gamma) to conduct Probabilistic Sensitivity Analysis (PSA), health economic models often require more complex evidence synthesis. When your transition probabilities or unit costs are derived from full Bayesian probability models—such as hierarchical models or network meta-analyses fit in Stan, JAGS, or WinBUGS—you can bypass the internal samplers and pass your posterior draws directly into the decision tree.

This vignette demonstrates how to extract posterior samples from a Stan model and format them for use as PSA inputs in CEdecisiontree.

1. Fitting the Bayesian Model in Stan

Suppose we are evaluating a new intervention against a standard of care. We have a simple decision tree where patients either experience a successful treatment or a failure.

First, we fit our evidence synthesis model in Stan. For this example, we assume a simple model where treatment success probabilities ($p$) and treatment costs ($c$) are estimated from trial data.

// example_model.stan
data {
  int<lower=0> N_standard;
  int<lower=0> successes_standard;
  int<lower=0> N_new;
  int<lower=0> successes_new;
  // ... cost data omitted for brevity ...
}
parameters {
  real<lower=0, upper=1> p_success_standard;
  real<lower=0, upper=1> p_success_new;
  real<lower=0> cost_new;
  real<lower=0> cost_standard;
}
model {
  // Priors
  p_success_standard ~ beta(1, 1);
  p_success_new ~ beta(1, 1);
  cost_new ~ gamma(100, 0.1);
  cost_standard ~ gamma(50, 0.1);

  // Likelihood
  successes_standard ~ binomial(N_standard, p_success_standard);
  successes_new ~ binomial(N_new, p_success_new);
  // ... cost likelihood omitted ...
}

We fit this model in R using rstan (or cmdstanr).

library(rstan)

# Fit the model (assuming data list 'stan_data' is prepared)
fit <- stan(file = 'example_model.stan', data = stan_data, 
            iter = 2000, warmup = 1000, chains = 4)

2. Extracting and Formatting Posterior Samples

Once the model has converged, we extract the joint posterior distribution. These draws represent our PSA parameter sets.

# Extract posterior samples as a data frame
post_draws <- as.data.frame(extract(fit))

# Inspect the first few draws
head(post_draws)
#>   p_success_standard p_success_new cost_standard cost_new     lp__
#> 1          0.621           0.754          512.4    1015.3  -152.1
#> 2          0.605           0.781          498.1     988.7  -150.3
#> 3          0.643           0.712          530.2    1054.1  -153.8

CEdecisiontree expects PSA inputs to be structured appropriately, typically as a data frame or a list of matrices where rows correspond to simulation runs (PSA iterations) and columns correspond to the nodes or branches of the decision tree.

Let's format these samples to match the labels of our tree branches.

n_sim <- nrow(post_draws)

# Create probability inputs for the decision tree
psa_probs <- data.frame(
  p_success_std = post_draws$p_success_standard,
  p_fail_std    = 1 - post_draws$p_success_standard,
  p_success_new = post_draws$p_success_new,
  p_fail_new    = 1 - post_draws$p_success_new
)

# Create cost inputs for the decision tree
psa_costs <- data.frame(
  c_success_std = post_draws$cost_standard,
  c_fail_std    = post_draws$cost_standard + 500, # e.g., added cost of failure
  c_success_new = post_draws$cost_new,
  c_fail_new    = post_draws$cost_new + 500
)

3. Defining the Decision Tree

Next, we define the structure of our decision tree using a transition matrix or a parent-child list, as is standard in CEdecisiontree.

library(CEdecisiontree)

# Define the tree structure (parent-child list)
tree_structure <- list(
  "Root" = c("Standard", "New"),
  "Standard" = c("Success_Std", "Fail_Std"),
  "New" = c("Success_New", "Fail_New")
)

# Convert to transition matrix
trans_mat <- child_list_to_transmat(tree_structure)

4. Evaluating the Model over the Posterior

Instead of calling internal sampling functions, we pass our pre-calculated psa_probs and psa_costs directly into the evaluation functions. We can iterate over the posterior draws to calculate the expected values for each MCMC iteration.

# Initialize an empty matrix to store expected costs for each intervention
expected_costs <- matrix(NA, nrow = n_sim, ncol = 2)
colnames(expected_costs) <- c("Standard", "New")

# Fold back the tree for each posterior draw
for (i in 1:n_sim) {
  
  # Extract the ith posterior draw for probabilities and costs
  current_probs <- as.numeric(psa_probs[i, ])
  current_costs <- as.numeric(psa_costs[i, ])
  
  # Calculate expected values for this specific MCMC draw
  ev <- dectree_expected_recursive(
    transmat = trans_mat,
    probs = current_probs,
    vals = current_costs
  )
  
  # Store the expected values for the decision nodes
  expected_costs[i, "Standard"] <- ev["Standard"]
  expected_costs[i, "New"]      <- ev["New"]
}

# The result is a distribution of expected costs that preserves the 
# joint posterior correlation structure estimated in Stan
summary(expected_costs)

Conclusion

By mapping the output of rstan::extract() (or the posterior package) to the branch labels of your tree, CEdecisiontree seamlessly integrates with Bayesian workflows. This allows you to leverage the sophisticated evidence synthesis capabilities of Stan while utilizing CEdecisiontree for rapid, transparent, spreadsheet-style model evaluation.