Multi-normal priors for non-linear formula syntax

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.

1 Like

Sorry, it looks like your question fell through.

I fear you cannot control it, although, you might look at make_standata (and maybe it’s code) to see how brms builds the b vector and use the same method.

What I would do, would be to force an empty prior and then add the multi_normal in the model block, i.e. (not tested, just a sketch):

stanvar(scode = "b ~ multi_normal(M, V);\n", block = "model")

(and use make_stancode to check that the code is built the way you expect it to be)

Best of luck!

Thanks for attempting to rescue this question from the dustbin of historical postings. I am not terribly optimistic there is a solution here (ie when you say “force an empty prior”, is there a way to prevent the brms error for unspecified priors on nonlinear parameters?) but I appreciate the go anyway. Forcing brms to accept these kind of contortions is probably like trying to take a Ferrari off-road when you should just take the raw Stan jeep.

1 Like

Actually, looking at docs for set_prior, you might be able to do prior_string("", check = FALSE), or - if empty string is problematic, just add something like target += 0…

But I agree that at some point using brms becomes too much of a hassle.

Hi @franzsf,

I am dealing with the same issue. In trying to reduce divergent transitions for a model discussed here, it was suggested that modeling correlations among population-level effects and group-level effects might help. However, I have not been able to implement the multi_normal() prior using brms.

Have you found a way to specify the multi-normal prior in the NL formula syntax? Would appreciate your feedback.

Thanks,
Meghna

I’m afraid I did not :( Ultimately it wasn’t critical for my purposes so I gave up. If someone has an answer, though, I’d still be interested in learning if this is possible.

OK - thanks for letting me know. I will update you if I manage to make it work in brms. Perhaps @paul.buerkner may have some insights? I came across this discussion, but am not able to implement the suggestion with the non-linear syntax.