Different Bayes factors from rstanarm vs BayesFactor

[gbiele]

I am not sure if I am doing something wrong, but I think the following analyses should result in the same BF calculated with the bayestestR and BayesFactor packages:

library(rstanarm)
library(bayestestR)
library(BayesFactor)

data = data.frame(y = (sleep$extra[1:10]-sleep$extra[11:20]))
null.interval = c(-1,1)
scale = sqrt(2)/2

model <- stan_glm(
  formula = y ~ 1,
  data = data,
  prior_intercept = cauchy(0, scale, autoscale = FALSE), iter = 20000,
)

BFa = bayesfactor_parameters(model, null = null.interval)
bfInterval = ttestBF(data$y, nullInterval = null.interval, rscale = scale)
BFb = bfInterval[2] / bfInterval[1]
BFa
BFb

However, the results are noticeably different, such that BF_BayestesterR = 6.282 and BF_BayesFactor = 2.108712.
[BayestesterR can also have difficulties with sampling from the prior].

The BFs are also not identical if I use a point Null:

BFa = bayesfactor_parameters(model, null = 0)
BFb = ttestBF(data$y, rscale = scale)
BFa
BFb

Here, BF_BayestesterR = 16.058 and BF_BayesFactor = 17.25888. This is obviously much closer.

Did I miss something that makes that the piors for the two analysis are not the same?

Any hint is very much appreciated.

I'll also post this Issue at https://github.com/richarddmorey/BayesFactor/issues.

[mattansb]

Hi @gbiele,

This isn't really an issue with bayestestR, but has to do with differences between rstanarm and BayesFactor.
I will try to give the gist nonetheless (:

This all has to do with Stan (via rstanarm, or brms) using a very different parameterization than BayesFactor.

Just to give an example:

• In rstanarm the priors are scaled, but they are scaled by SD_y = the unconditional SD of y.
• In BayesFactor the prios are also scaled, but they are scaled by σ_y = the conditional SD of y.

(There are many many other differences - how estimation of sigma is done, how random effects are treated... etc.)

In other words, this isn't a bug - just different flavors of priors.

Hope this helps!

[gbiele]

See also this nice explanation of the issue at hand by @richarddmorey: richarddmorey/BayesFactor#156 .
Richard shows that the main reason for the diverging results is that ttestBF puts the chauchy prior on the standardized effect, whereas the stan_glm model I specified did put the same chauchy prior on the non-standardized effect size.

[jhmaindonald]

It is worth pointing out that, not only are the Bayes factors different between the two functions, the Bayes factors from bayesfactor_parameters() vary (at least for the simple example that I have tried) quite dramatically from one run to the next.

btrfun <- function(pairedDiffs, seed=NULL, data=sleep1){
set.seed(seed)
model <- stan_glm(
  formula = pairedDiffs ~ 1, data=data, 
)
BF0 <- bayestestR::bayesfactor_parameters(model, null = 0)
## BF <- bayestestR::bayesfactor_parameters(model, null = c(-.75,.75))
fac.BF0 <- round(exp(as.data.frame(BF0)[,"log_BF"]),1)
fac.BF0
## fac.BFint <- round(exp(as.data.frame(BF)[,"log_BF"]),1)
## fac.BFint
}
sleep1 <- data.frame(pairedDiffs=with(sleep, extra[group==2]-extra[group==1]))
btrBF <- numeric(10)
for(i in 1:10)btrBF[i] <-btrfun(pairedDiffs=pairedDiffs, data=sleep1)
btrBF
## [1] 28.9 56.4 31.8 23.7 40.7 41.5 12.7 27.3 49.1 21.5

Such comparing of different priors as I tried did not lead to more consistent output.

BayesFactor::tttestBF() returns a consistent 17.26.

With the NULL is an interval, I obtained more consistent results from bayesfactor_parameters() , e..g.,
with an interval that extends 45 minutes either side of zero.

## [1] 9.8 8.4 9.7 5.9 8.2 6.5 6.6 7.0 6.8 8.3

[mattansb]

@jhmaindonald Yes - should also have mentioned a 3rd difference between BayesFactor and Stan + bayestestR:

1. Different parameterization of the model
2. Different type of BF (marginal vs univariate likelihood)
   And
3. How the likelihood is quantified:
   • BayesFactor actually integrates over the prior parameter space - that is why the results are fixed
   • Stan + bayestestR is based on fitting a density function to the prior and posterior samples.

Since the prior/posterior are sampled, and since empirical density functions are unstable (esp. at the tails), this leads to the instability of BFs using this method.
This is also why interval BFs are more stable as they are based on probability functions (CDF = cumulative distribution function), which are more stable.

To increase the stability of the point null BF, one can sample more than the default 1,000 samples * 4 chains - say 10 times more (e.g., 4,000 samples * 10 chains).

[richarddmorey]

@jhmaindonald As @mattansb mentions, this will be due to the approximation used. Such Bayes factors are notoriously difficult to estimate. BayesFactor does a lot of work (analytical integration; restrictions on priors) to keep the integrals easy and one-dimensional, at the expense of generality.

You can read about some of the issues in our paper "Using MCMC chain outputs to efficiently estimate Bayes factors" (https://www.sciencedirect.com/science/article/abs/pii/S0022249611000666) which gives a more efficient way of doing it than a density estimate under some circumstances (the conditional-marginal density estimation, CMDE). See Figure 2 for evidence of the extreme unreliability of the density-fitting approach (left panel) versus the CMDE approach (right panel).

The CMDE can also be used for intervals (section 6.1), but most of the time it is not worth it because those are fairly stable anyway.

[mattansb]

Thanks @richarddmorey
I will examine is we can implement the CMDE method for point-null BFs.

[gbiele]

hi @mattansb,

first thanks for the extremely fast response!

So rstanarm will scale the prior even though the prior is defined by cauchy(0, scale, autoscale = FALSE)?

I'l try to look more into how the error variance is estimated in BayesFactor and see if I can set prior_aux in rstanarm appropriately.

Or were you implying that it won't be possible to set up an rstanmodel that matches what BayesFactor does?

I am aware that priors have an undesirably large effect on BF, but I am surprised how large the differences became when I played a bit around with other d.

[mattansb]

So rstanarm will scale the prior even though the prior is defined by cauchy(0, scale, autoscale = FALSE)?

No, if you set autoscale = FALSE it won't scale them. But if autoscale = TRUE it will scale differently then BayesFactor (which does scale them).

I'l try to look more into how the error variance is estimated in BayesFactor and see if I can set prior_aux in rstanarm appropriately.
Or were you implying that it won't be possible to set up an rstanmodel that matches what BayesFactor does?

Exactly. The parameterization makes the priors and resulting Bayes factors uncomparable, even for the simplest of models.

[gbiele]

Exactly. The parameterization makes the priors and resulting Bayes factors uncomparable, even for the simplest of models.

OK, even if this seems rather undesirable to me.* If you have the time, could you expand a bit with what you mean with "The parameterization"? (I'll also look at the source code for ttestBF to hopefully learn what I am missing here.)

Edit: *To be clear, I think this is an undesirable property of the Bayes factor, not of its implementation in any package.

[mattansb]

If you have the time, could you expand a bit with what you mean with "The parameterization"?

Simply put: a correlation can be computed in two ways:

1. Cov(x,y)/(S_x * S_y)
2. b_x * S_y / S_x

In an OLS of ML frame work you will get the exact same correlation in both parameterizations.
But in a Bayesian framework, you would be putting different priors on different parameters, so the resulting correlations from both methods can be different, reflecting differences in their priors. (If you know some Stan, you can try this out yourself.)

---

I think this is an undesirable property of the Bayes factor, not of its implementation in any package.

I don't like the framing as Bayes factors as "very sensitive". Bayes factors do one very specific thing: they ask "which of these two hypotheses is more supported by the data". What BFs are sensitive too, is the hypothesis, which in the Bayesian world is denoted by the priors and the parameterization of the model - this is a desirable property.

(The same is true for p-values. I don't see anyone complaining that p-values are too sensitive to the likelihood function that was used...)

The issue is that researchers often don't know how to translate their hypothesis into a formal mathematical model.
Or they have non-hypotheses like "there will some difference between some or all of these 9 groups", which can really only be tested by showing that there isn't no difference (note the double negative) which is a NHST.

</rant>

[jhmaindonald]

While I see lots of problems with p-values, I'd judge that if different parameterizations give different values, one needs to be looking carefully to see what gives the best sense. There are very effective diagnostic plots that one can use (but checking for independence is another matter!)

[mattansb]

Yes, I strongly agree - prior, parameterization and hypotheses should be considered very carefully.

Often ROPE tests or loo-IC / WAIC comparisons might be better suited for the "test" at hand.

[gbiele]

Here richarddmorey/BayesFactor#156 Richard Morey shows that the key thing I did not recognize was that if one uses stan_glm the prior is on the un-scaled effect size, whereas on ttestBF the prior is on the standardized effect (d).

[mattansb]

By default, both rstanarm and BayesFactor use a prior on the standardized effect size - but they standardized differently, which was the point I made above about different parameterizations.

But yes, in you sepecific example you "turned off" the standardization in rstanarm by setting autoscale = FALSE.

[gbiele]

But yes, in you specific example you "turned off" the standardization in rstanarm by setting autoscale = FALSE.

and despite this one gets somewhat different BF even if one "pre-scales" the outcome variable. This is likely due to the fact that the sigma estimated with stan_glm is not equal to the sd of the data (even when using a practically non-informative prior for sigma).

[jhmaindonald]

Thanks for the various comments. Surely both the “Bayes Factor” vignette and the help page for bayesfactor_parameters() should, when they give an example, explicitly give the 1,000 samples * 4 chains parameters when the model is specified. As it is, for anyone not familiar with rstanarm and bayestestR, these are not easy to chase down. Also, it would be useful to provide a warning in cases where the Bayes factor estimate is highly unstable. The estimate might be checked against the upper limit given in the Sellke at al 2001 American Statistician paper (~22.4 for p=0.0028 as in the paired differences p-value for the sleep data); the BayesFactor::ttestBF() value of 17.26 is consistent with this. Values that are widely different at the lower end might also be flagged. Also, encourage users to run the calculations more than once. The publishers have now, after stalling for quite a time on going ahead with a 4th end of my book with John Braun “Data Analysis and Graphics Using R”, indicated that they wish to proceed. One last-minute inclusion comments briefly on the Bayes Factor connection to t-tests that go beyond highlighting the implications of the Sellke et al result. Our stance is pretty much an acceptance that p-values will rule for a while yet — but please understand that they do not have many of the implications that are commonly given to them. It is well to point out that different priors will give different Bayes Factors, but wildly different Bayes Factors that result from the way that the prior/posterior are sampled risk encouraging a “you can prove anything with statistics” mentality.

…

On 3/10/2021, at 23:00, Mattan S. Ben-Shachar ***@***.***> wrote:

[mattansb]

@jhmaindonald Thanks for these suggestions.

I will add the chains and samples explicitly to the vignettes.

The estimate might be checked against the upper limit given in the Sellke at al 2001 American Statistician paper

This would require computing the equivalent p-values - something that is easy for simple designs, but becomes exponentially harder with model complexity, meaning is would be complex and computationally time consuming - all for just a warning.

Instead, I think it would be a good idea to add a warning when testing point nulls if the number of samples is too small (we already do this for bayesfactor_models()).

[...] wildly different Bayes Factors that result from the way that the prior/posterior are sampled risk encouraging a “you can prove anything with statistics” mentality.

I'm sure you agree that this is not unique to Bayes factors - "When a measure becomes a target, ...".
Anyone who sets out to prove something with statistics shouldn't use (and doesn't really need) statistics.