Many thanks to @jonah , and apologies for the delayed response (I had to take it all in and then write code). Based on those pointers, I have written a script that takes a brms model and a parameter-specifying regex as input, extracts the analytical priors and sampled posteriors, and plots both in the same plot using ggplot2 and ggdist. It’s a script I’d eventually like to turn into a handy helper function, given I like visualizing analytical priors together with empirical posteriors. But there’s a problem.
As far as I can understand, the marginal density of any pairwise correlation in a prior correlation matrix R should depend not only on the eta parameter of the LKJ distribution, but also on K, namely the dimension of R. Yet, even though my script (which renames the distribution and applies marginalize_lkjcorr() as shown by Jonah) appears to parse the distribution correctly, the plotted prior curves fail to reflect differences in K between correlation matrices with the same eta. Here is mock dataset mydata, constructed for easy illustration of the problem.
Let’s fit a brms model to these data:
require(brms); require(cmdstanr)
Cores <- 6 # your mileage may vary
Chains <- Cores
brmsmod <- brm(y ~ (1+x1|G1) + (1+x1+x2+x3|G2) + x1 + x2 + x3, family = bernoulli, data = mydata, seed = 2026, refresh = 2e3,
prior = prior(lkj(6), class = cor) + prior(normal(0, 2.5), class = b) + prior(normal(0, 5), class = Intercept) + prior(normal(0, 0.5), class = sd)
, backend = "cmdstanr", warmup = 2e3, cores = Cores, chains = Chains, iter = ceiling(25e3/Chains)+2e3, threads = threading(Cores/Chains, static = T))
And below is my plotting script, followed by the graphical output evidencing the problem:
require(ggplot2); require(ggdist); require(stringi)
pars_regex <- "^(b|sd|cor)"
# Extract posterior draws for the specified parameters:
as.data.frame(head(posttab.raw <- as_draws_df(brmsmod)[grep(pars_regex, names(as_draws_df(brmsmod)))], 3))
posttab <- data.frame(Parameter = rep(names(posttab.raw), each = nrow(posttab.raw)), postval = unlist(posttab.raw))
# Extract the table of priors:
(priors <- prior_summary(brmsmod))
for(i in 2:nrow(priors)){
for(j in which(names(priors) %in% c("prior","lb","ub"))){
if(priors[i, j] == ""){priors[i, j] <- priors[i-1, j]}else{next} } }
(priors$prior <- gsub("_corr_cholesky\\(", "(", priors$prior)) #Gabry changes the name like this.
# Convert LKJ priors into a format ggdist supposedly understands:
(Priors <- parse_dist(priors))
if(any(startsWith(names(posttab.raw), "cor")) ){
(groupingVars <- unique(priors$group[priors$class == "L" & priors$group != ""]))
(corrpairs <- character())
for(g in groupingVars){
(Rdim <- length(unique(priors$coef[priors$group == g]))-1)
(eta <- as.integer(stringi::stri_extract(priors$prior[priors$class == "L" & priors$group == g][1], regex = "\\d+(?=\\)$)"), ")"))
(Priors <- marginalize_lkjcorr(Priors, K = Rdim, predicate = (class == "L" & group == g)))
(varying <- c("Intercept",
unlist(strsplit(gsub("^1[\\W]+", "",stringi::stri_extract(as.character(formula(brmsmod))[1], regex = paste0("(?<=\\()[^|)]+(?= \\| ", g,"\\))") ), perl = T), split = " + ", fixed = T))))
newpairs <- character()
for(i in 1:(length(varying)-1)){for(j in (i+1):length(varying)){newpairs <- c(newpairs, paste0(varying[i], "__", varying[j]))}}
(corrpairs <- c(corrpairs, newpairs))
(Priors <- rbind(Priors, Priors[rep(which(Priors$class == "L" & Priors$group == g)[1], times = length(newpairs)),]))
(Priors[tail(1:nrow(Priors), length(newpairs)), c("class","coef")] <- data.frame(rep("cor", length(newpairs)), newpairs))
rownames(Priors)[tail(1:nrow(Priors), length(newpairs))] <- paste0(g, "_cor", 1:length(newpairs))
}
}
# Finish a dataframe with all the analytical priors for ggplotting:
(tab.pr <- Priors[Priors$class == "Intercept" | Priors$coef != "",])
(tab.pr$Parameter <- gsub("^Intercept$", "b_Intercept", gsub("^b_+","b_", gsub("_+$", "", paste0(tab.pr$class, "_", tab.pr$group, "__", tab.pr$coef), perl = T), perl = T), perl = T))
(tab.pr <- tab.pr[grep(pars_regex, tab.pr$Parameter, perl = T),])
# Plot all priors and posteriors:
ggplot(posttab, aes(postval, Parameter)) + geom_vline(xintercept = 0, lwd = 0.5, color = "gray80") + #
stat_slab(alpha = 1, fill = "#b3cde0", lwd = 0.6, , color = "black", point_size = 1) + geom_hline(yintercept = 1:length(unique(posttab$Parameter))) +
stat_slab(aes(xdist = .dist_obj, y = Parameter), data = tab.pr, fill = NA, lwd = 0.5, color = "black", lty = 5, inherit.aes = F) +
scale_y_discrete(limits = rev(names(posttab.raw)), rev(names(posttab.raw))) + scale_thickness_shared() + theme(axis.title = element_blank()) +
coord_cartesian(xlim = 8 * c(-1,1)) + scale_x_continuous(name = NA, breaks = seq(-8,8,2))
As you can see, grouping factor G1’s single correlation parameter is portrayed as having the same prior density as those of G2, even though G2’s prior correlation matrix has dimension 4 as opposed to only 2, and both share the same eta value of 6.
The problem appears not to be caused by a misparse ― as can easily be verified by running the code step by step, the .dist_obj vector (as well as the .args and .dist ones) does show the correct argument values, reflecting the difference in matrix dimension between G1 and G2. Yet somehow stat_slab plots the exact same prior curve for G1 and G2. Is it a software bug or am I missing something?
