Can `powerscale()` be used to obtain prior draws?

General
1

In {bayestestR} we use prior draws together with posterior draws to estimate Savage-Dickey ratios, Order restricted Bayes factors, and support intervals.

Currently, the way to do this, for example with a brmsfit object is to re-run the model with prior="only" (this is what bayestestR::unupdate() does). This is rather time consuming, as the full HMC must be run again.

So I was wondering if somehow powerscale() can be used to get prior samples? My thinking was that setting alpha=0 for the likelihood would “cancel” it’s effect, and so we would get prior samples, but this doesn’t seem to do the job.

library(brms)
library(priorsense)

mod <- brm(mpg ~ hp, data = mtcars, 
           prior = set_prior("normal(1, 2)", class = "b"),
           backend = "cmdstanr", silent = TRUE)


ps1 <- powerscale(mod, component = "likelihood", alpha = 0, resample = TRUE)

# expecting a normal distribution centered around 1 with a sd of 2:
hist(ps1$b_hp)

image
image483×299 3.64 KB

mean(ps1$b_hp)
#> [1] -0.07144655

sd(ps1$b_hp)
#> [1] 0.01998533

Thanks!

2

powerscale() uses importance weighting and re-sampling. Here the posterior is much narrower than the prior, and by doing weighted resampling of the posterior draws you can’t get draws in the range where the prior has most of it’s mass. In general importance sampling based power scaling works when alpha is close to 1. You should have got some warnings about high Pareto-k’s.

Not that you asked, but bridge-sampling has lower variability than Savage-Dickey ratios and don’t require re-running MCMC, so I don’t see a reason to use Savage-Dickey ratios. But then, I also don’t see reasons to use Bayes factors :D

3

Thanks Aki!

(I personally don’t use Bayes factors to test point null hypotheses, so we’re in agreement about Savage-Dickey ratios ;) but I do like order restricted Bayes factors!)