bayprior handles prior specification, elicitation, pooling, and
diagnostics; it does not fit posterior models itself. This vignette
shows how to take a bayprior prior object’s fitted hyperparameters and
hand them to rstanarm or brms for the
model-fitting step, using the same historical-trial data as the MAP
Priors example in the package’s main vignette and companion paper.
Both rstanarm and brms are Suggested, not
Imported, dependencies: this vignette’s code only runs if both are
installed.
library(bayprior)
map <- map_prior(
y = c(-0.85, -0.62, -1.10),
se = c(0.25, 0.30, 0.28),
outcome_type = "single_arm_log_odds",
label = "Historical control (3 trials)"
)
mu <- map$fit_summary$mean
sigma <- map$fit_summary$sd
c(mu = mu, sigma = sigma)
#> mu sigma
#> -0.8618678 0.3349206map$dist is "normal", and
mu/sigma are its mean and SD on the
log-odds scale – exactly the scale a
logistic-regression intercept uses. This is what makes the handoff below
possible without any distributional-family conversion: both packages
already expect a Normal prior on that scale for an intercept term.
rstanarm::normal() evaluates its
location/scale arguments as ordinary R
expressions, so mu and sigma can be passed
directly. Setting autoscale = FALSE is important:
rstanarm otherwise may rescale the SD you supplied relative
to the data before fitting, which would silently change the prior
bayprior computed.
library(rstanarm)
# A toy current-trial dataset, sized only to keep this vignette's build
# time short -- substitute your own trial data in practice.
set.seed(1)
n_trial <- 40
p_trial <- 0.32
dat <- data.frame(y = rbinom(n_trial, 1, p_trial))
# Wrapped in tryCatch so an incomplete local C++/Stan toolchain degrades
# this vignette gracefully (a build failure here would otherwise fail
# R CMD check on such a machine) rather than as a claim this code is
# untested -- it has been run successfully end to end where the toolchain
# is complete.
fit_rstanarm <- tryCatch(
stan_glm(
y ~ 1,
data = dat,
family = binomial(link = "logit"),
prior_intercept = normal(location = mu, scale = sigma, autoscale = FALSE),
chains = 1,
iter = 500,
refresh = 0,
seed = 1
),
error = function(e) {
message("Model fit skipped (local Stan toolchain issue): ", conditionMessage(e))
NULL
}
)
if (!is.null(fit_rstanarm)) print(fit_rstanarm, digits = 3)
#> stan_glm
#> family: binomial [logit]
#> formula: y ~ 1
#> observations: 40
#> predictors: 1
#> ------
#> Median MAD_SD
#> (Intercept) -0.794 0.211
#>
#> ------
#> * For help interpreting the printed output see ?print.stanreg
#> * For info on the priors used see ?prior_summary.stanregbrms::prior() is not a drop-in
replacement here: it captures its argument by non-standard evaluation
and stores it as literal text, so
brms::prior(normal(mu, sigma), class = "Intercept") would
store the string "normal(mu, sigma)" – the variable
names, not their values – and fail at fit time.
brms::prior_string() builds the prior specification from an
already-evaluated R string instead, and is the reliable way to pass
numeric values programmatically:
library(brms)
prior_spec <- prior_string(
paste0("normal(", mu, ",", sigma, ")"),
class = "Intercept"
)
prior_spec
#> Intercept ~ normal(-0.861867808691721,0.334920636951119)
# See the rstanarm chunk above for why this is wrapped in tryCatch.
fit_brms <- tryCatch(
brm(
y ~ 1,
data = dat,
family = bernoulli(link = "logit"),
prior = prior_spec,
chains = 1,
iter = 500,
refresh = 0,
seed = 1,
silent = 2
),
error = function(e) {
message("Model fit skipped (local Stan toolchain issue): ", conditionMessage(e))
NULL
}
)
if (!is.null(fit_brms)) summary(fit_brms)
#> Family: bernoulli
#> Links: mu = logit
#> Formula: y ~ 1
#> Data: dat (Number of observations: 40)
#> Draws: 1 chains, each with iter = 500; warmup = 250; thin = 1;
#> total post-warmup draws = 250
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept -0.84 0.23 -1.32 -0.45 1.01 97 66
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).| Package | Correct call | Common mistake |
|---|---|---|
rstanarm |
normal(location = mu, scale = sigma, autoscale = FALSE) |
Omitting autoscale = FALSE lets rstanarm silently
rescale the SD |
brms |
prior_string(paste0("normal(", mu, ",", sigma, ")"), class = "Intercept") |
prior(normal(mu, sigma), ...) stores the variable names
as text, not their values |
Both handoffs apply specifically to an intercept
term representing the same quantity map_prior() (or an
elicited prior on the log-odds scale) was fit to – not to
slope/coefficient priors (class = "b" in brms,
prior = in rstanarm), which are on a different
parameter and should not receive these values unmodified.