"LogBF = Inf" problem

[Masharipov]

Thanks for the useful package!

I was used bayestestR to compare sensetivity (TPR) of different methods (M1,M2,M3,M4,M5) averaged over the levels of SNR.
This is the code:

contrasts(Data$Method) <- contr.equalprior_pairs 
contrasts(Data$SNR) <- contr.equalprior_pairs 

fit_bayes <- stan_glm( TPR_prc ~ Method + SNR, data = Data, iter = 5000, warmup = 1000) 
pairs_bayes <- pairs(emmeans(fit_bayes, ~Method))
postprob <- describe_posterior(pairs_bayes,rope_range = c(-1,1))
bf <- bayesfactor_parameters(pairs_bayes,fit_bayes ,null=c(-1,1))
postprob$LogBF = bf$log_BF

My data:
Data.csv

The problem is that I get LogBF = Inf for pairwise comparisons:

I understand that the posterior distributions are very narrow and BF-->Inf.
The maximum floating-point number is 1.797693e+308.
Therefore, the maximum possible logBF is 709.7827.
However, JASP can report logBF > 709.7827 (e.g., logBF ≈ 4000-7000 for my data).

I as far as I understand, bayestestR first calculate BF and then takes log(BF). In contrast, JASP first calculate logBF and then calculate BF = exp(logBF).
Maybe this capability could be added to the bayestest R package?

[mattansb]

This issue is due to {logspline} not supporting the return of log-densities. I have sent the maintainer a feature request. If they will not respond, I will build this support into bayestestR.

[mattansb]

From the maintainer of logspline:

I can do this.
It will probably be another week or two before I have enough time to do this carefully, update helpfiles etc.

---

Once this is implemented, I will make the required changes in bayestestR.

[Masharipov]

Thanks!

[mattansb]

Okay, so this actually wasn't an issue in {logspline} (that's what happens when I read github issues on my phone...).

The "issue" is that the posterior probability for the defined null interval is infinitesimal - out of 16,000 posterior samples, not 1 falls within the null interval (except for M2 - M3 and M4 - M5):

library(bayestestR)
library(emmeans)

Data <- read.csv("C:\\Users\\user\\Downloads\\Data.csv", sep = ";")
Data$Method <- C(factor(Data$Method), contr.equalprior_pairs)
Data$SNR <- C(factor(Data$SNR), contr.equalprior_pairs)

fit_bayes <- stan_glm(TPR_prc ~ Method + SNR,
                      data = Data,
                      iter = 5000, warmup = 1000, cores = 4)

pairs_bayes <- pairs(emmeans(fit_bayes, ~ Method))

# how many samples in the [-1, +1] interval?
insight::get_parameters(pairs_bayes) |> 
  sapply(\(x) sum(-1 < x & x < 1))
#> M1 - M2 M1 - M3 M1 - M4 M1 - M5 M2 - M3 M2 - M4 M2 - M5 M3 - M4 M3 - M5 M4 - M5 
#>       0       0       0       0   16000       0       0       0       0   16000 

Not only that, but the smallest (of 16,000 samples!) is quite far from the null interval:

insight::get_parameters(pairs_bayes) |> 
  sapply(max)
#>     M1 - M2     M1 - M3     M1 - M4     M1 - M5     M2 - M3     M2 - M4     M2 - M5     M3 - M4     M3 - M5     M4 - M5 
#>  -2.4825258  -2.2830532 -29.6783006 -29.3708621   0.6660124 -26.8417212 -26.4680872 -27.0671175 -26.7161429   0.7777120 

And overall the whole mass is in a narrow band:

describe_posterior(pairs_bayes, test = NULL)
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |           95% CI
#> -------------------------------------
#> M1 - M2   |  -2.90 | [ -3.11,  -2.69]
#> M1 - M3   |  -2.66 | [ -2.86,  -2.45]
#> M1 - M4   | -30.13 | [-30.34, -29.92]
#> M1 - M5   | -29.76 | [-29.96, -29.55]
#> M2 - M3   |   0.24 | [  0.04,   0.45]
#> M2 - M4   | -27.23 | [-27.43, -27.02]
#> M2 - M5   | -26.86 | [-27.06, -26.65]
#> M3 - M4   | -27.47 | [-27.68, -27.26]
#> M3 - M5   | -27.10 | [-27.31, -26.89]
#> M4 - M5   |   0.37 | [  0.16,   0.57]

Thus, the posterior probability of that range is effectively 0:

# right tail probability
insight::get_parameters(pairs_bayes) |> 
  sapply(\(x) {
    lsp <- logspline::logspline(x) 
    1 - logspline::plogspline(-1, lsp)
  })
#> M1 - M2 M1 - M3 M1 - M4 M1 - M5 M2 - M3 M2 - M4 M2 - M5 M3 - M4 M3 - M5 M4 - M5 
#>       0       0       0       0       1       0       0       0       0       1

Since this is not the case for the prior differences (i.e., $p_{prior} > 0$), which represent the tested hypotheses (which look like this BTW:

library(ggplot2)
library(ggdist)

pairs_bayes_prior |> 
  insight::get_parameters() |> 
  datawizard::data_to_long() |>
  ggplot(aes(value, name)) + 
  # facet_wrap(~name) + 
  stat_slab(aes(fill = after_stat(abs(x) < 1)), n = 2000) + 
  scale_fill_discrete("Hypothesis", labels = c("Alternative", "Null")) + 
  labs(y = "Difference")

Does this look the the hypotheses you wish to compare?)

... and since 0 divided by any positive number is $\infty$, you get the BF = Inf (I'm skipping a step in the calculation, but the point is still true).

So, assuming this is the correct null interval for your hypothesis (is it?) and that the extremally wide default priors adequately represent the hypotheses you're trying to test (I very much doubt that), there is no issue at all - you have extreme evidence against the null!

---

I am not aware of JASP providing such BFs for contrasts.

[Masharipov]

Thanks for the detailed answer!

So, assuming this is the correct null interval for your hypothesis (is it?) and that the extremally wide default priors adequately represent the hypotheses you're trying to test (I very much doubt that), there is no issue at all - you have extreme evidence against the null!

I tried to change the priors. I supposed that prior difference in sensetivity (TPR) between methods should be around 0%, ranging from -100% to 100%. Intercept should be around 50%, ranging from 0% to 100%.

fit_bayes <- stan_glm(TPR_prc ~ Method + SNR,
                      data = Data,
                      prior = normal(0, rep(50,8), autoscale = FALSE),
                      prior_intercept = normal(50,50, autoscale = FALSE),
                      iter = 5000, warmup = 1000, cores = 4)

I am not aware of JASP providing such BFs for contrasts.

For pairwise comparisons in JASP, I used the default t-tests with a Cauchy (0, r = 1/sqrt(2)) prior. For instance, LogBF for M1-M4 difference = 1248.

As far as I understand, JASP uses this equation to calculate BF from t-values (Rouder et al., 2009):

Since JASP does not calculate probability densities (e.g. using logspline) but calculates BF analytically, 'inf' did not arise.

I believe "inf" will be unavoidable for my data if I use batestestR (calculate BF based on probability densities).

[mattansb]

I tried to change the priors. [...]

But a null range of ±1 points of change is reasonable?

Anyway, the formula used in JASP is also a probability density - just one derived analytically - BUT
as far as I know, JASP's contrasts are just model-less t-tests, and do not allow for null intervals (only points). So I see no reason those would be analogous to what you've attempted here.

Enjoy your $BF = \infty$!