← Back to blog

From a bayprior MAP prior to a fitted model: handing priors to rstanarm and brms

Bayesian bayprior rstanarm brms

bayprior handles prior specification, elicitation, pooling and diagnostics. It does not fit posterior models itself. For that step most teams already use rstanarm or brms, so the practical question is how to pass a prior from one to the other without it being altered on the way.

This post walks through the handoff for a common case: a historical control-arm response rate, combined across trials into a meta-analytic-predictive (MAP) prior, then used as the intercept prior in a logistic model for the current trial.

rstanarm and brms are suggested packages of bayprior, not dependencies, and historical_effect_sizes() needs metafor. Install what you need: install.packages(c("rstanarm", "brms", "metafor")). Both modelling packages need a working Stan toolchain.

Step 1: Derive the prior from historical trials

Start from raw event counts per historical trial. historical_effect_sizes() converts them into the estimate and standard error that map_prior() expects. For a control-arm response rate the outcome type is "single_arm_log_odds", and the columns must be named xi (events) and ni (patients):

library(bayprior)

trials <- data.frame(xi = c(12, 8, 15), ni = c(40, 35, 50))
es     <- historical_effect_sizes("single_arm_log_odds", trials)

map <- map_prior(
  y            = es$y,
  se           = es$se,
  outcome_type = "single_arm_log_odds",
  label        = "Historical control (3 trials)"
)

print(map)
plot_tau_posterior(map)

If you already have per-trial estimates and standard errors on the log-odds scale, you can pass them to map_prior() directly as y and se.

The outcome_type choice matters. It selects the default prior on the between-trial heterogeneity tau, and "log_or" and "single_arm_log_odds" get different default scales. plot_tau_posterior() shows how much the data say about heterogeneity, which is worth a look before you borrow strength from historical trials.

Prefer to explore without writing code? The Shiny app includes a MAP prior module. Try it in the live demo, or run it locally with bayprior::run_app().

Step 2: Read off the hyperparameters

The MAP prior is a Normal distribution on the log-odds scale. Its mean and SD are stored in fit_summary:

mu    <- map$fit_summary$mean
sigma <- map$fit_summary$sd
c(mu = mu, sigma = sigma)

A logistic-regression intercept is also on the log-odds scale, and both modelling packages accept a Normal prior for it. No change of distribution family is needed, only the right call.

Step 3: Handoff to rstanarm

rstanarm::normal() evaluates location and scale as ordinary R expressions, so mu and sigma can be passed directly. Set autoscale = FALSE. Otherwise rstanarm may rescale the SD you supply relative to the data, which silently changes the prior bayprior computed.

library(rstanarm)

# Toy current-trial data: 40 patients, binary response
set.seed(1)
dat <- data.frame(y = rbinom(40, 1, 0.32))

fit_rstanarm <- stan_glm(
  y ~ 1,
  data            = dat,
  family          = binomial(link = "logit"),
  prior_intercept = normal(location = mu, scale = sigma, autoscale = FALSE),
  chains          = 4,
  iter            = 2000,
  refresh         = 0,
  seed            = 1
)

print(fit_rstanarm, digits = 3)

Step 4: Handoff to brms

brms::prior() is not a drop-in replacement. It captures its argument by non-standard evaluation and stores it as text, so prior(normal(mu, sigma), class = "Intercept") stores the variable names, not their values, and fails at fit time. prior_string() builds the specification from an already-evaluated string, which is the reliable way to pass numbers programmatically:

library(brms)

prior_spec <- prior_string(
  paste0("normal(", mu, ",", sigma, ")"),
  class = "Intercept"
)
prior_spec

fit_brms <- brm(
  y ~ 1,
  data    = dat,
  family  = bernoulli(link = "logit"),
  prior   = prior_spec,
  chains  = 4,
  iter    = 2000,
  refresh = 0,
  seed    = 1
)

summary(fit_brms)

Two limits to keep in mind

Only a single Normal prior hands off this way. A mixture, such as one from robust_prior(), is not a single Normal, so it is not covered by this recipe.

Keep the justification with the model

The prior in your model and the prior in your documentation should be the same object. Pass the same map to prior_report() so the report that goes in the submission package records exactly what was handed to Stan:

prior_report(
  prior         = map,
  output_format = "html",
  trial_name    = "TRIAL-001",
  sponsor       = "Example Pharma Ltd",
  author        = "N.P., Biostatistician"
)

Summary: rstanarm needs normal(location = mu, scale = sigma, autoscale = FALSE); brms needs prior_string() with the evaluated values.

References


Questions or corrections — contact@reprostats.org.