# Ggdist::marginalize\_lkjcorr() not working as intended

**URL:** <https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658>\
**Category:** Other\
**Tags:** visualisation\
**Created:** [October 12, 2026, 10:45pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658 "2026-10-12T22:45:59Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![blokeman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blokeman/32/16236_2.png) [@blokeman](https://discourse.mc-stan.org/u/blokeman)\
**Post date:** [October 12, 2026, 10:45pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658/1 "2026-10-12T22:45:59Z")

</div>

I can’t visualize the lkj priors of my correlation parameters because I can’t find how to convert the **lkj\_corr\_cholesky(eta)** priors into .dist objects that **stat\_slabinterval()** would understand. The **marginalize\_lkjcorr()** function should technically do this, but it seems not to be working:

```no-highlight
#Load libraries:
require(brms)
require(cmdstanr)
require(ggplot2)
require(ggdist)

# Create mock data:
N = 100
datums <- data.frame(y = rbinom(N, 1, 0.5), x1 = rnorm(N, 0, 1), x2 = rnorm(N, -5, 2), x3 = rnorm(N, 5, 3), g = sample(LETTERS, N, replace = T))

# Fit model:
mod <- brm(y ~ (x1|g) + x1 + x2 + x3, family = bernoulli, prior = prior(normal(0, 2.5), class = b) + prior(normal(0, 5), class = Intercept) + prior(normal(0, 1), lb = 0, class = sd) + prior(lkj(2), class = cor) 
  , data = datums, warmup = 2e3, cores = Cores, chains = Chains, iter = ceiling(25e3/Chains)+2e3, threads = threading(Cores/Chains, static = T), backend = "cmdstanr")

# Extract analytical prior distributions from model:
mypriors <- marginalize_lkjcorr(parse_dist(prior_summary(mod)), K = 2)

# Try to plot the priors:
ggplot(mypriors) + stat_slab(aes(xdist = .dist_obj, y = class))

```

```no-highlight
Error in `layer_slabinterval()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error in `get()`:
! object 'rlkj_corr_cholesky' of mode 'function' was not found
Run `rlang::last_trace()` to see where the error occurred.

```

While these simple analytical priors could of course be plotted using a more basic function such as **curve()**, that’s not so helpful when you want to combine both the analytical priors and the posterior densities in the same plot while maintaining the full functionality of the **ggdist** library.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [October 12, 2026, 10:59pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658/2 "2026-10-12T22:59:36Z")

</div>

I think this is probably needs to be fixed in ggdist, so I recommend opening an issue at [Issues · mjskay/ggdist · GitHub](https://github.com/mjskay/ggdist/issues). Until it’s fixed, does this avoid the problem for you?

```no-highlight
library(dplyr)

mypriors <- prior_summary(mod) |>
  as.data.frame() |>
  mutate(prior = sub("^lkj_corr_cholesky\\(", "lkj(", prior)) |>
  parse_dist(prior) |>
  marginalize_lkjcorr(K = 2) |>
  filter(!is.na(.dist))

ggplot(mypriors) + stat_slab(aes(xdist = .dist_obj, y = class))

```

I think the issue might be that brms is reporting the name as `lkj_corr_cholesky` not `lkj` in `prior_summary()`, but `ggdist::parse_dist()` is expecting `lkj`. I haven’t actually run the code above, but I think it’s approximately right and you can modify as needed.

---

<div class="post-metadata">

**Author:** ![blokeman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blokeman/32/16236_2.png) [@blokeman](https://discourse.mc-stan.org/u/blokeman)\
**Post date:** [October 20, 2026, 4:47pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658/3 "2026-10-20T16:47:35Z")

</div>

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](https://github.com/jkruohon/StatsMisc/blob/main/mydata.RData) is mock dataset _mydata,_ constructed for easy illustration of the problem.

Let’s fit a brms model to these data:

```no-highlight
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:

```no-highlight
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))

```

 ![Rplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/1/1/1129847bd67f26a292b28f1b6bbdf077284dfd5f.png)

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?

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [October 20, 2026, 5:49pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658/4 "2026-10-20T17:49:15Z")

</div>

I don’t actually think anything is necessarily going wrong here. I think your code and ggdist (after the parsing fix) are both doing the right thing, but the two priors are close enough that your plot can’t show the difference. I ran your script on your data (with fewer iterations, since the priors don’t depend on the posterior) and checked the objects:

```no-highlight
              Parameter dist dens_at_0 sd
1 cor_G1 __Intercept__ x1 lkjcorr_marginal(2, 6) 1.354 0.277
2 cor_G2 __Intercept__ x1 lkjcorr_marginal(4, 6) 1.466 0.258
... (all six G2 rows identical to row 2)

```

These match the marginal distributions if you do it analytically, which I think for each correlation is Beta(eta − 1 + K/2, eta − 1 + K/2) but rescaled to (−1, 1) instead of (0, 1). So that’s Beta(6, 6) for G1 and Beta(7, 7) for G2. When you go from K = 2 to K = 4 it only adds 1 to the shape parameters, so with eta = 6 the change is very small. I had to zoom in to see the difference:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/9/f/9f1b9929b08ece14ffea081fce3af761784840dd.png)

In the plot you shared I think two things are making this harder to see. The x-axis goes from -8 to 8, so a 0.02 difference in SD is less than a pixel. And `scale_thickness_shared()` scales the priors against the posteriors, so each prior curve is only a few pixels tall.

For the helper function, it may be helpful to give the correlations their own panel with `xlim = c(-1, 1)`. If you want to see the dependence on K more clearly, I would try with eta = 1. Then K = 2 gives a uniform prior and K = 4 is Beta(2, 2).

---

<div class="post-metadata">

**Author:** ![blokeman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/blokeman/32/16236_2.png) [@blokeman](https://discourse.mc-stan.org/u/blokeman)\
**Post date:** [October 20, 2026, 9:30pm UTC](https://discourse.mc-stan.org/t/ggdist-marginalize-lkjcorr-not-working-as-intended/41658/5 "2026-10-20T21:30:20Z")

</div>

Hot damn sir, you’ve resolved everything! There’s not even a hurry with that ggdist fix anymore, given that your/our script works around the issue just fine.