Prior/link selection for sigma in distributional Gaussian model

I want to fit a distributional model with brms along the following lines:

df <- data.frame() # some columns include: value (continuous, zero-mean), afactor (3 levels), bfactor (7 levels), ...
frml <- bf(
  value ~ 0 + Intercept,
  sigma ~ 0 + Intercept + afactor * bfactor
)
priors <- c(
  set_prior("normal(0,0.1)", class = "b", dpar = ""),
  set_prior("normal(0,0.1)", class = "b", coef = "Intercept", dpar = ""),
  set_prior("normal(0,20)", class = "b", dpar = "sigma", lb = 0),
  set_prior("normal(0,40)", class = "b", coef = "Intercept", dpar = "sigma")
)
bfit <- brm(
  frml,
  data = df,
  family = gaussian(),
  prior = priors
)

Importantly, I want to set a prior so that sigma is biased to be zero.

The gaussian family in brms uses a log link for sigma. My understanding is that putting this together with my last two priors above results in a sigma whose prior biases sigma to be 1. However, I want the prior on sigma to bias it towards 0.

That is, if the identity link could be used for sigma (which I realize can create issues with negative sigma), I would want to set the priors to be half normal distributions. Alternatively, I’m open to using a log link on sigma while specifying a prior which is the log of a half normal distribution. Any ideas about how to achieve this?

Maybe think of this by how a normal prior on the log space is the same as a lognormal prior on the identity space. Let’s compare how two priors, prior(normal(0, 1)) and prior(normal(-1, 1)) on the log space would look like on the identity space using their lognormal equivalents. We can use functions from {ggdist} and {ggplot2} to show them in a plot.

library(ggdist)
library(ggplot2)

c(prior(lognormal(0, 1)),
  prior(lognormal(-1, 1))) |> 
  parse_dist() |> 

  ggplot(aes(y = prior, dist = .dist_obj)) +
  stat_halfeye(.width = 0.95, linewidth = 1, n = 501 * 2) +
  scale_y_discrete(NULL, expand = expansion(mult = 0.1)) +
  coord_cartesian(xlim = c(0, 10))

Both are pretty close to what you want, but prior(normal(-1, 1)) is closer.

If you go all the way to prior(normal(-2, 1)) on the log space (i.e., prior(lognormal(-2, 1)) on the identity space), the entire inner 95% prior density is between 0 and 1.

qlnorm(c(0.025, 0.975), meanlog = -2, sdlog = 1)
[1] 0.0190638 0.9607548

Make a ggdist plot like the one above to see what it looks like. You might like it.

But anyway, you’ll never get a normal prior on the log space to go all the way to zero on the identity space, but you can get pretty dang close with priors as simple as normal(-1, 1) or normal(-2, 1). IMO, this is preferable to fitting a model with the idenitty link for sigma, or to experimenting with normal priors on the log space with extremely large values on their sd parameters (such as your normal(0, 40), above).

I think what I want is a normal prior on the identity space, but only the upper half since sigma must be positive.

I want to test whether the variable afactor has an effect on the variability of the variable value. I want the prior to bias the results towards no effect. I’m concerned that with what is essentially a lognormal prior, the prior will bias the results towards there being an effect of afactor, and at that point I won’t be able to tell whether any results I see are really coming from the data.