Noting the stanvar
documentation has a brief illustration of how this is done in the usual brms syntax, and noting that this is probably more straightforward to execute in raw Stan, I was hoping for some pointers (probably from @paul.buerkner) on how a defined covariance matrix can be used to construct a multi-normal prior using brms non-linear formula syntax.
Example code below works well for the normal brms syntax to create a prior-only model that can recover the correlation between the distribution of intercepts and slopes, though I’m not sure how the ordering of the population-effect vector in “b” is (or can be) specified.
library(data.table)
library(brms)
library(ggplot2)
# Create fake data that is appx. linear
fake_data <- data.table("x" = c(0,1,2,3,4,5,6),
"y" = c(-1,3,4,1,2,6,7))
ggplot(fake_data) +
geom_point(aes(x=x,y=y))
# Fit model
fake_data_model <- brm(bf(y ~ 0 + Intercept + x),
data = fake_data,
family = "gaussian",
prior = c(prior(normal(0,10), class = "b", coef = "Intercept"),
prior(normal(0,1), class = "b"),
prior(normal(0,1), class = "sigma")),
chains = 4,
cores = 4,
iter = 1000,
warmup = 500)
# Parameters and covariance matrix
summary(fake_data_model)
cov_matrix <- vcov(fake_data_model)
# Prior-only model that draws intercepts and slopes from multi-normal
# with known covariance matrix
prior_only_model <- brm(bf(y ~ 0 + Intercept + x),
data = fake_data,
family = "gaussian",
prior = c(prior(multi_normal(M, V), class = "b"),
prior(normal(1.72,0.39), class = "sigma")),
stanvars = stanvar(c(0.41,0.91),
name = "M",
scode = " vector[K] M;") +
stanvar(cov_matrix,
name = "V",
scode = " matrix[K, K] V;"),
sample_prior = "only",
chains = 4,
cores = 4,
iter = 1000,
warmup = 500)
It’s unclear, however, how this can be done when the same model is written in brms’ non-linear formula syntax, i.e., how can stanvars and priors linked to allow the same multi-normal prior here:
nlf_syntax_model <- brm(bf(y ~ yIntercept + ySlope*x,
yIntercept ~ 1,
ySlope ~ 1,
nl = TRUE),
data = fake_data,
family = "gaussian",
prior = c(prior(normal(0.41,1.12), nlpar = "yIntercept"),
prior(normal(0.91,0.31), nlpar = "ySlope"),
prior(normal(1.72,0.39), class = "sigma")),
chains = 4,
cores = 4,
iter = 1000,
warmup = 500)
Thanks for any guidance you can provide.