Priors for log(sigma) in distributional Gaussian model

General
,
1

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
2

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

3

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.

4

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

2 Likes
5

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)

image

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 [2]
  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)?

6

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.

7

Thanks.

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

image
image881Ă—584 12 KB

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))

image
image620Ă—584 7 KB

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
8

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
9

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
10

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))

image

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”)

2 Likes