Value
stanvars object to be used as argument to the
stanvars argument of a brm model.
Details
To declare mixture priors in a brm
model requires two steps: First, the mixture densities need to
be converted by the adapter function mixstanvar into a
stanvars object which is passed to the stanvars
argument of the brm function. Doing so
extends the Stan code and data generated by
brm such that the mixture densities can be
used as priors within the brm model. The
second step is to assign parameters of the
brm model to the mixture density as prior
using the set_prior command of brms.
The adapter function translates the mixture distributions as
defined in R to the respective mixture distribution in
Stan. Within Stan the mixture distributions are named in
accordance to the R functions used to create the respective
mixture distributions. That is, a mixture density of normals
created by mixnorm is referred to as
mixnorm_lpdf in Stan such that one can refer to the density
as mixnorm within the set_prior
functions (the suffix _lpdf is automatically added by
brm). The arguments to these mixture
distributions depend on the specific distribution type as follows:
| Density | Arguments |
mixbeta(w, a, b) | w weights, a shapes, b shapes |
mixgamma(w, a, b) | w weights, a shapes, b inverse scales |
mixnorm(w, m, s) | w weights, m means, s standard deviations |
mixmvnorm(w, m, sigma_L) | w weights, m means, sigma_L cholesky factors of covariances |
These arguments to the mixture densities refer to the different
density parameters and are automatically extracted from the
mixtures when converted. Important here are the argument names as
these must be used to declare the mixture prior. For each mixture
to convert as part of the ... argument to mixstanvar
a label is determined using the name given in the list. In case no
name is given, then the name of the R object itself is
used. To declare a prior for a parameter the mixture distribution
must be used with it's arguments following the convention
label_argument. Please refer to the examples section for an
illustration.
Note: Models created by brm do use by
default a data-dependent centering of covariates leading to a shift
of the overall intercept. This is commonly not desirable in
applications of this functionality. It is therefore strongly
recommended to pass the option center=FALSE as argument to
the brms formula created with the bf
function as demonstrated with the example below.
Examples
if (FALSE) { # \dontrun{
# The mixstanvar adapter requires the optional packages brms and glue
stopifnot(require("brms"), require("glue"))
# Assume we prefer a logistic regression MCMC analysis rather than a
# beta-binomial analysis for the responder endpoint of the ankylosing
# spondylitis (AS) example. Reasons to prefer a regression analysis is
# to allow for baseline covariate adjustments, for example.
map_AS_beta <- mixbeta(c(0.62, 19.2, 57.8), c(0.38, 3.5, 9.4))
# First we need to convert the beta mixture to a respective mixture on
# the log odds scale and approximate it with a normal mixture density.
map_AS_samp <- rmix(map_AS_beta, 1E4)
map_AS <- mixfit(logit(map_AS_samp), type="norm", Nc=2)
# Trial results for placebo and secukinumab.
trial <- data.frame(n=c(6, 24),
r=c(1, 15),
arm=factor(c("placebo", "secukinumab")))
# Define brms model such that the overall intercept corresponds to the
# placebo response rate on the logit scale. NOTE: The use of
# center=FALSE is required here as detailed in the note above.
model <- bf(r | trials(n) ~ 1 + arm, family=binomial, center=FALSE)
# to obtain detailed information on the declared model parameters use
# get_prior(model, data=trial)
# declare model prior with reference to mixture normal map prior...
model_prior <- prior(mixnorm(map_w, map_m, map_s), coef=Intercept) +
prior(normal(0, 2), class=b)
# ... which must be made available to brms using the mixstanvar adapter.
# Note that the map_AS prior is labeled "map" as referred to in the
# previous prior declaration.
analysis <- brm(model, data=trial, prior=model_prior,
stanvars=mixstanvar(map=map_AS),
seed=365634, refresh=0)
# Let's compare the logistic regression estimate for the probability
# of a positive treatment effect (secukinumab response rate exceeding
# the response rate of placebo) to the direct beta-binomial analysis:
hypothesis(analysis, "armsecukinumab > 0")
post_secukinumab <- postmix(mixbeta(c(1, 0.5, 1)), r=15, n=24)
post_placebo <- postmix(map_AS_beta, r=1, n=6)
pmixdiff(post_secukinumab, post_placebo, 0, lower.tail=FALSE)
# The posterior probability for a positive treatment effect
# is very close to unity in both cases.
} # }