Different BF for interval Null from bayestestR and BayesFactor

[gbiele]

Edit: @mattansb has responded to this questions here: https://github.com/easystats/bayestestR/discussions/459. Maybe @richarddmorey also wants to weigh in?

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/easystats/bayestestR/issues.

[richarddmorey]

Are the hypotheses on the same space? I'm not sure about bayestestR, but in the case of BayesFactor the intervals are on the standardized effect size. I doubt this is the case for stan_glm.

[gbiele]

That's a good point (which I should have thought of myself).

But the BFs are still quite different if I set the scale to 1 before analysing the data (code below): BF_BayestesterR = 0.125 BF_BayesFactor = 2.108712.

Do you think this might have to do with difference in the models used by stan_glm and ttestBF (not as far as I understand it), or do you see other potential reasons for the disagreement?

One thing I can think of is that bayesfactor_parameters could run into problems because it samples from he prior, which I expect to be hard due to the fat tails of the cauchy.

Or am I totally missing something here?

code:

data = data.frame(y = scale(sleep$extra[1:10]-sleep$extra[11:20], center = F))

null.interval = c(-1,1)
scale = sqrt(2)/2

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

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

[richarddmorey]

Yes, the models are different. The stan code for the model the BayesFactor uses would be something like:

data {
   int N;
   vector[N] y;
   real<lower=0> prior_scale;
}
parameters {
   real<lower=0> sigma2;
   real delta;
}
transformed parameters {
  real mu;
  real sigma = sqrt(sigma2);
  mu = delta * sigma;

}
model {
  y ~ normal(mu, sigma); 
  delta ~ cauchy(0, prior_scale);
  target += -log(sigma2);
}
generated quantities {
  real prior_delta;
  prior_delta = cauchy_rng(0, prior_scale);
}

...but the stan code for the model stan_glm uses would be:

data {
   int N;
   vector[N] y;
   real<lower=0> prior_scale;
   real<lower=0> prior_scale_sigma;
}
parameters {
   real<lower=0> sigma;
   real mu;
}
transformed parameters {
  real delta;
  real sigma2 = sigma * sigma;
  delta = mu / sigma;
}
model {
  y ~ normal(mu, sigma); 
  mu ~ cauchy(0, prior_scale);
  sigma ~ exponential(prior_scale_sigma);

}
generated quantities {
  real prior_delta;
  real prior_mu;
  real prior_sigma;

  prior_mu = cauchy_rng(0, prior_scale);
  prior_sigma = exponential_rng(prior_scale_sigma);
  prior_delta = prior_mu / prior_sigma;
}

Also, the Bayes factor from {BayesFactor} is on the parameter 𝛿, while the Bayes factor from {bayestestR} is on the parameter μ. So the priors, models, and hypotheses are all different. We could compute a Bayes factor on 𝛿 using the same model as stan_glm uses, but due to the weakly informative, data dependent prior on σ this Bayes factor would be highly suspect.

If you'd like to explore this, you can use the R code below along with the stan code listings above.

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,
)

stan_glm_samples = do.call(
  rbind,
  rstan::As.mcmc.list(model$stanfit)
  )

stan_glm_delta = stan_glm_samples[,"alpha[1]"] / stan_glm_samples[,"aux"]
hist(stan_glm_delta)

stan_fit <- rstan::stan(
  file = 'stan_glm.stan',
  data = list(
    y = data$y,
    prior_scale = scale,
    prior_scale_sigma = 1/sd(data$y),
    N = length(data$y)
  )
)

stan_samples = rstan::extract(stan_fit, permute = TRUE)
hist(stan_samples$delta)

# Check to make sure they're the same

qqplot(stan_glm_delta, stan_samples$delta)
abline(0,1)

qqplot(stan_glm_samples$alpha[,1], stan_samples$mu)
abline(0,1)

qqplot(stan_glm_samples$aux, stan_samples$sigma)
abline(0,1)

# Estimate the Bayes factor on mu 
prior_prob = mean(stan_samples$prior_mu> null.interval[1] & stan_samples$prior_mu < null.interval[2])
prior_odds = prior_prob / (1-prior_prob)
posterior_prob = mean(stan_samples$mu > null.interval[1] & stan_samples$mu < null.interval[2])
posterior_odds = posterior_prob / (1-posterior_prob)

bf_est = posterior_odds / prior_odds
1/bf_est
# Confirm against bayestestR
bayesfactor_parameters(model, null = null.interval)

# Run the JZS model
jzs_fit <- rstan::stan(
  file = 'jzs_t.stan',
  data = list(
    y = data$y,
    prior_scale = scale,
    N = length(data$y)
    )
  )

jzs_samples = rstan::extract(jzs_fit, permute = TRUE)
hist(jzs_samples$delta)

# Estimate the JZS Bayes factor on delta
prior_prob_jzs = mean(jzs_samples$prior_delta > null.interval[1] & jzs_samples$prior_delta < null.interval[2])
prior_odds_jzs = prior_prob_jzs / (1-prior_prob_jzs)
posterior_prob_jzs = mean(jzs_samples$delta > null.interval[1] & jzs_samples$delta < null.interval[2])
posterior_odds_jzs = posterior_prob_jzs / (1-posterior_prob_jzs)

bf_est_jzs = posterior_odds_jzs / prior_odds_jzs
1/bf_est_jzs

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

[richarddmorey]

I see there's plenty of discussion over on https://github.com/easystats/bayestestR/discussions/459, so I'm going to consider this closed.

[gbiele]

Thanks a lot for the thorough response @richarddmorey!

I had actually tried to get around the problem of having the prior on mu in stan_glm vs. on delta in ttestBF by scaling the data, but I had messed up the scaling. When I do the scaling correctly like this

y = sleep$extra[1:10]-sleep$extra[11:20]
data = data.frame(y = y/sd(y))

The base factors get much more similar (BF_BayestesterR = 2.577, BF_BayesFactor = 2.108712.), though still not identical. I understand now better why this is the case, so thanks again.