# Priors for log(sigma) in distributional Gaussian model

I can fit a distibutional model with brms like this (example from Estimating Distributional Models with brms):

``````group <- rep(c("treat", "placebo"), each = 30)
symptom_post <- c(rnorm(30, mean = 1, sd = 2), rnorm(30, mean = 0, sd = 1))
dat1 <- data.frame(group, symptom_post)

fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),
data = dat1, family = gaussian())
``````

`prior_summary(fit1)` says that the default prior for `sigma_intercept` is `student_t(3, 0, 2.5)` and that for the group effect, it is `(flat)`.

How do I interpret these priors? brms fits log(sigma) with the linear predictor `~ 1 + group`. Do the priors correspond to log(sigma) or to sigma? and what about the prior for the `group` effect on sigma?

1 Like

just like any other family with a link function… the priors are for the linear model specified for log sigma.

What does this imply for the prior of sigma if the prior of log(sigma) is `student_t(3, 0, 2.5)`?

I suppose this puts maximum probability of sigma = 1, and equal probability of sigma = 0.5 and sigma = 2. Is this a reasonable prior for sigma?
I cant get my head around how to actually plot the density.

how about you call brm with the argument sample_prior=“only”. Then you get a plain sample from that thing.

2 Likes

Good idea. I had to set a prior for `b`

``````fit1 <- brm(bf(symptom_post ~ group, sigma ~ 1), prior = c(set_prior("normal(0,2)")),
data = dat1, family = gaussian(), sample_prior = "only")

plot(fit1)
`````` `sigma_intercept` looks `student_t` like…

`epred_draws` allow us to investigate what values are passed to the distribution.

``````expand(dat1, group) |>
+   add_epred_draws(fit1, dpar = "sigma") |>
+   median_qi()
# A tibble: 2 × 11
# Groups:   group 
group    .row .epred .epred.lower .epred.upper sigma sigma.lower sigma.upper .width .point
<chr>   <int>  <dbl>        <dbl>        <dbl> <dbl>       <dbl>       <dbl>  <dbl> <chr>
1 placebo     1  0.298        -8.29         8.85  1.04    0.000158       2083.   0.95 median
2 treat       2  0.187        -8.17         8.87  1.04    0.000158       2083.   0.95 median
# … with 1 more variable: .interval <chr>
``````

The 97.5th percentlile for sigma is 2083.

How can I set a reasonable prior for log(sigma)?

The basic strategy is to experiment with different choices of priors on log(sigma) until you arrive at one that yields a useful pushforward on `sigma`. To speed up this process, you can examine the pushforward by sampling a bunch of random draws from a distribution directly in R (e.g. `rnorm`) and then exponentiating those and summarizing the result in tabular or graphical form. To get started, note that normal priors on `log(sigma)` are by definition log-normal priors on `sigma`, so if you think a log-normal might be a reasonable prior then you should be able to achieve what you want with normal priors on `sigma`.

Thanks.

`plot(density(exp(rnorm(1000, 0, 1))), xlim = c(0, 10))` looks reasonable.

I have not even been able to draw the density of `exp(rt(1000, 3)*2.5)`.
e.g. `plot(density(exp(rt(1000, 3)*2.5)), xlim = c(0, 10))`

Is that default prior someting that was overlooked (or postponed) when `brms` was extended to fit distributional models, or is it sometimes meaningful? (@paul.buerkner )

1 Like

Your call to hist includes xlims that are applied after the breaks are computed, so you’re seeing only one bin really. Remove the limits to see the full distribution.

1 Like

If you want to quickly visualize (transformations of) distributions, my biased suggestion would be to try out {ggdist} with {distributional}, which can do transformations of analytical distributions for you (it uses numerical differentiation to apply the necessary Jacobian correction).

{ggdist} will also help you pick limits to see the bulk of the distribution. In your case its defaults won’t work well because the log-t distribution has such fat tails. The default for an infinite upper limit is the 0.999 quantile, which is nearly 30,000 for t(3, 0, 1), making the chart illegible:

``````library(ggplot2)
library(ggdist)
library(distributional)

data.frame(sigma = c(exp(dist_normal(0,1)), exp(dist_student_t(3, 0, 1)))) |>
ggplot(aes(xdist = sigma, y = as.character(sigma))) +
stat_halfeye()
`````` But you can use `p_limits` to change the upper limit to be the 0.975 quantile (i.e. the upper limit of the 95% central interval), which gives you a better view:

``````data.frame(sigma = c(exp(dist_normal(0,1)), exp(dist_student_t(3, 0, 1)))) |>
ggplot(aes(xdist = sigma, y = as.character(sigma))) +
stat_halfeye(p_limits = c(0, 0.975))
`````` 5 Likes

That is perfect. I had the idea that it should be trivial to draw the density of the exponent of a distribution (without simulation), but could not figure it out. I’m somehow glad that it does not seem to be entirely trivial :)

I am now able to draw the default prior (t(3, 0, 2.5):

``````data.frame(sigma = c(exp(dist_normal(0,1)), exp(dist_student_t(3, 0, 2.5)))) |>
ggplot(aes(xdist = sigma, y = as.character(sigma))) +
stat_slab(p_limits = c(0, 0.8))
`````` Notice that 20% of the density is > ~12.

Is that a reasonable default prior, or is it something I should consider opening and issue about?

If I remove the `sigma ~ group` term, the default prior for `sigma` (now not `log(sigma)`) is simply `student_t(3, 0, 2.5)` with `lb = 0`. This seems like very different default priors for the distributional and non-distributional models, but I do not know whether the difference is important.

(I would have added a plot, but I could not figure out how to draw a “half student t”)

1 Like