Posterior predictions for sigma in a distributional model

Interfaces brms
1

Hi,

I’m trying to get posterior predictions for a distributional model using posterior_predict. It’s easy to get a prediction for the means this way, however I couldn’t find a way to get predictions for the sigmas.

The code below reproduces this problem: pp ends up with values of symptom_post but not of sigma (code based on the vignette).

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)
head(dat1)

fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),  data = dat1, family = gaussian())
pp = posterior_predict(fit1, newdata = data.frame(group = c('treat','placebo')))
  • Operating System: Ubuntu
  • brms Version: 2.15.0

Thanks :)

1 Like
2

I’m having the same problem. I’d also be very interested for an answer to this!

3

@eladzlot @atsentinella

This query has come up on a few previous occasions on this forum. And it is indeed not obvious in the online documentation.

Anyway, …

str(fit1) will show the specific distributional parameters for your specific likelihood (family).

And in your case those are “mu” and “sigma”.

So each component or family-specific parameter could be extracted for instance by specifying each distributional parameter using dpar as follows:

sigma_component<-conditional_effects(fit1,dpar="sigma",method="posterior_epred",robust=FALSE)

mu_component<-conditional_effects(fit1,dpar="mu",method="posterior_epred",robust=FALSE)

Now both family-specific components combined for the expected mean in this instance:

both<-conditional_effects(fit1,method="posterior_epred",robust=FALSE),

which is saying that dpar=TRUE for both components.

You could of course use (…, method=“posterior_predict”) if you think that was appropriate,

So, your …

posterior_predict(fit1, newdata = data.frame(group = c(‘treat’,‘placebo’)))

would be the predicted response for new data with both components (mu, sigma) combined as expected.

Hope this helps.

4

Thank Milani!

So if I understand correctly, to use new data to predict only sigma I would use predict(fit1, newdata, dpar = "sigma")?

Except, when I try this (predict/posterior_predict) there is no difference in results between dpar = "sigma",dpar = "mu" and no dpar (while setting the same seed).

Also, with the example you gave, I’m failing to see the difference between the “mu_component” and “both”. Sigma is clearly different, but the others appear the same. Also, dpar=TRUE also doesn’t seem to be accepted in any function here.

Sorry if I’m missing something obvious.

5

Not sure why you would be using posterior_predict in this way but are you sure that the following doesn’t work?

posterior_predict(fit1, dpar="sigma",newdata = data.frame(group = c('treat','placebo')))

and …

posterior_predict(fit1, dpar="mu",newdata = data.frame(group = c('treat','placebo')))

6

Both lines give the same results (with the same seed). It seems like the “dpar” argument is ignored here.

The reason I am trying this is to make predictions of sigma with new data.

7

I think it should be fitted(…, dpar = “sigma”) or equivalently posterior_epred(…, dpar = “sigma”).

8

Thanks Paul! Yes, that works as expected.