# Is there an easy non-modelling way for visualizing different prior distributions?

**URL:** <https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757>\
**Category:** Modeling\
**Tags:** rstanarm, brms\
**Created:** [August 1, 2021, 1:22pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757 "2021-08-01T13:22:30Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![di4oi4](https://avatars.discourse-cdn.com/v4/letter/d/779978/32.png) [@di4oi4](https://discourse.mc-stan.org/u/di4oi4)\
**Post date:** [August 1, 2021, 1:22pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/1 "2021-08-01T13:22:30Z")

</div>

Normally distributed priors are easy to be simulated and visualized without doing any modelling (which takes some extra time). For example, with ggplot:

```
dist = rnorm(1000, 15, 8)
ggplot(NULL, aes(dist)) + geom_density()

```

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/28dd0dc07be98777a7355d237ed363672f076bcb.png)

**Is there an easy and similar option for checking cauchy, laplace, student t and other prior distributions?**

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [August 1, 2021, 5:22pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/2 "2021-08-01T17:22:05Z")

</div>

You can plot simulations from any available probability distribution, or alternatively, you can directly plot the density function for a given distribution. For example:

```r
library(tidyverse)
library(VGAM) # For laplace distribution
library(patchwork) # For plot layouts
theme_set(theme_bw())

# Sample from distribution
p1=ggplot(tibble(x=rt(1000, df=4)), aes(x)) +
  geom_density()

p2=ggplot(tibble(x=rlaplace(1000, location=0, scale=4)), aes(x)) +
  geom_density()

p3=ggplot(tibble(x=rlnorm(1000, meanlog=0, sdlog=0.5)), aes(x)) +
  geom_density()

# Plot multiple sets of random draws from a distribution
p4=ggplot() +
  replicate(10, geom_density(data=tibble(x=rlaplace(100, location=0, scale=4)), aes(x)))

list(p1,p2,p3,p4) %>% wrap_plots(nrow=2)

```

![](https://i.imgur.com/xZsrBcF.png)

```r
# Plot density function
p5 = ggplot() +
  stat_function(fun=dcauchy, args=list(location=1, scale=2), 
                xlim=c(-10,10))

p6=ggplot() +
  stat_function(fun=dnorm, args=list(mean=2, sd=3), 
                xlim=c(-10,15))

p7=ggplot() +
  stat_function(fun=dt, args=list(df=4), 
                xlim=c(-20,20), n=1000)

p8=ggplot() +
  stat_function(fun=dlaplace, args=list(location=0, scale=4), 
                xlim=c(-20,20), n=1000)

list(p5,p6,p7,p8) %>% wrap_plots(nrow=2)

```

![](https://i.imgur.com/5VvYrz5.png)

Created on 2021-08-01 by the [reprex package](https://reprex.tidyverse.org) (v2.0.0)

---

<div class="post-metadata">

**Author:** ![di4oi4](https://avatars.discourse-cdn.com/v4/letter/d/779978/32.png) [@di4oi4](https://discourse.mc-stan.org/u/di4oi4)\
**Post date:** [August 14, 2021, 8:12am UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/3 "2021-08-14T08:12:27Z")

</div>

This is very useful! Does the same approach allow plotting different lkj\_corr\_cholesky priors?

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [August 14, 2021, 10:00pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/4 "2021-08-14T22:00:59Z")

</div>

For univariate density plots, you can plot the density for an LKJ prior for 2x2 correlation matrices, or marginal plots of LKJ priors for larger correlation matrices. Below are examples.

```r
library(tidyverse)
theme_set(theme_bw() +
            theme(plot.title=element_text(size=rel(1), hjust=0.5)))
library(ggdist)

# Plot LKJ Density
lkj_dfunc = function(K, etas=c(0.1, 0.5, 1, 1.5, 2, 5)) {
  
  etas = set_names(etas)
  cor = seq(-1,1,length=200)
  
  etas %>% 
  map_df(
    ~ dlkjcorr_marginal(cor, K=K, eta=.x)
  ) %>% 
  bind_cols(x=cor) %>% 
  pivot_longer(cols=-x, names_to='eta', values_to='density') %>% 
  ggplot(aes(x=x, y=density)) + 
    geom_area(colour=hcl(240,100,30), fill=hcl(240,100,65)) +
    facet_wrap(vars(eta)) +
    scale_y_continuous(expand=expansion(c(0,0.05))) +
    coord_cartesian(ylim=c(0,2)) +
    labs(title=paste0(ifelse(K==2, "", "Marginal "),
                      "LKJ correlation distribution for K=", K,
                      " and various values of eta"),
         x="Correlation") 
}

lkj_dfunc(2)

```

![](https://i.imgur.com/00FLMaZ.png)

```r
lkj_dfunc(5)

```

![](https://i.imgur.com/6jHoARD.png)

```r
# Plot density of random draws from LKJ distribution
lkj_rfunc = function(K, N=1000, etas=c(0.1, 0.5, 1, 1.5, 2, 5)) {
  
  etas = set_names(etas)
  
  etas %>% 
    map_df(
      ~ rlkjcorr_marginal(N, K=K, eta=.x)
    ) %>% 
    pivot_longer(cols=everything(), names_to='eta', values_to='corr') %>% 
    ggplot(aes(x=corr)) + 
      geom_density(colour=hcl(240,100,30), fill=hcl(240,100,65)) +
      facet_wrap(vars(eta)) +
      scale_y_continuous(expand=expansion(c(0,0.05))) +
      labs(title=paste0(ifelse(K==2, "D", "Marginal d"),
                        "ensity of random draws from LKJ distribution for K=", K,
                        " and various values of eta"),
           x="Correlation")
}

lkj_rfunc(2)

```

![](https://i.imgur.com/46Iq0wq.png)

```r
lkj_rfunc(5)

```

![](https://i.imgur.com/ebEXdnM.png)

Created on 2021-08-14 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)

---

<div class="post-metadata">

**Author:** ![andymilne](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andymilne/32/21046_2.png) [@andymilne](https://discourse.mc-stan.org/u/andymilne)\
**Post date:** [August 15, 2021, 4:06am UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/5 "2021-08-15T04:06:21Z")

</div>

A related question: Is there an easy way to plot what priors look like after an inverse link function has been applied – e.g., what does a normal(0, 1) prior on the logit scale look like on the probability scale?

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [August 16, 2021, 6:50pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/6 "2021-08-16T18:50:52Z")

</div>

Below is an example for normally distributed log(odds) priors centered at zero and with a range of standard deviations.

Note that as the standard deviation moves downward from one towards zero, most of the log(odds) density piles up near zero, (i.e., odds near 1:1), so most of the probability density is near 0.5.

On the other hand, when the standard deviation is greater than one, most of the probability quickly starts piling up near 0 and 1. For example, for SD=5, 69% of the density is in log(odds) less than -2 or greater than 2 (`1 - diff(pnorm(c(-2,2), 0, 5))`), meaning 69% of the probability density is less than 0.12 or greater than 0.88. For SD=10, 75% of the log(odds) density is less than -3 or greater than 3, meaning 75% of the probability density is less than 0.05 or greater than 0.95.

```r
library(LaplacesDemon) # For invlogit function
library(tidyverse)
theme_set(ggthemes::theme_tufte(base_size=13) +
            theme(axis.text.y=element_blank(),
                  axis.ticks.y=element_blank(),
                  panel.border=element_rect(fill=NA)))
library(patchwork) # For plot layout

set.seed(10)
d = c(0.1, 0.2, 0.5, 1, 2, 5, 10) %>% set_names() %>%  
  map_df(~ tibble(`log(odds)` = rnorm(1000, 0, .x), 
                  probability = invlogit(`log(odds)`)),
         .id="sd") %>% 
  mutate(sd = factor(sd, levels=unique(sd))) 

p1 = d %>% 
  ggplot() +
    geom_density(aes(x=`log(odds)`, y=..scaled..)) +
    facet_grid(rows=vars(sd), scales="free_y", switch="y") +
    scale_y_continuous(expand=expansion(c(0,0.05))) +
    labs(x="log(odds)", y="Standard Deviation") +
    theme(strip.text.y.left=element_text(angle=0))

p2 = d %>% 
  ggplot() +
    geom_density(aes(x=probability, y=..scaled..), adjust=0.5) +
    facet_grid(rows=vars(sd), scales="free_y") +
    scale_y_continuous(expand=expansion(c(0,0.05))) +
    scale_x_continuous(expand=expansion(c(0.001,0.001)), breaks=c(0,0.5,1)) +
    labs(x="probability", y=NULL) +
    theme(strip.text.y=element_text(angle=0))

p1 + p2 + 
  plot_layout(widths=c(5,1)) + 
  plot_annotation(title="Normally distributed log(odds) priors and transformation to probability scale", 
                  theme=theme(plot.title=element_text(hjust=0.5)))

```

![](https://i.imgur.com/5zLQ9hx.png)

Created on 2021-08-16 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)

---

<div class="post-metadata">

**Author:** ![andymilne](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andymilne/32/21046_2.png) [@andymilne](https://discourse.mc-stan.org/u/andymilne)\
**Post date:** [August 16, 2021, 8:52pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/7 "2021-08-16T20:52:47Z")

</div>

That’s beautiful! Thank you.

---

<div class="post-metadata">

**Author:** ![andymilne](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andymilne/32/21046_2.png) [@andymilne](https://discourse.mc-stan.org/u/andymilne)\
**Post date:** [August 16, 2021, 9:30pm UTC](https://discourse.mc-stan.org/t/is-there-an-easy-non-modelling-way-for-visualizing-different-prior-distributions/23757/8 "2021-08-16T21:30:59Z")

</div>

Part of the reason for asking was because I was interested in what prior would result in a flat probability distribution in a probit regression. Substituting this

```
d = c(0.1, 0.2, 0.5, 1, 2, 5, 10) %>% set_names() %>%  
    map_df(~ tibble(`log(odds)` = rnorm(1000000, 0, .x), 
                    probability = pnorm(`log(odds)`,0,1)),
           .id="sd") %>% 

```

into your code gives this lovely plot:

 ![Probit_priors](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/072151c189900650c0939290760f860743dcc5c6.png)

So a standard normal prior it is! I expect this is well-known analytically, but I had been unable to find a reference for it. Note that for a logistic regression (logit link), the unit logistic distribution prior is flat on the probability scale, and this distribution can be approximated by `student_t(7.763, 0, 1.566)` (Northrup and Gerber, 2018) – also easy to demo with @joels code. I mention all this simply because it may be of general interest.