jomulder / BFpack

BFpack can be used for testing statistical hypotheses using the Bayes factor, a Bayesian criterion originally developed by sir Harold Jeffreys.
https://bfpack.info/
14 stars 4 forks source link

Problem when using the BF() command on the two sample t-test #25

Closed tpourmohamad closed 1 month ago

tpourmohamad commented 3 years ago

There are two issues I seem to have found when using the t-test. One is for the paired case, and the other is for the unpaired case.

For the paired case, it seems to be that the model parameter is called "mu", but when I try to specify mu in the hypothesis, the BF() command remarks that mu cannot be found. Is it supposed to be called something else? For example, running

library(BFpack)

set.seed(4)

n <- 100
x <- rnorm(n, 50)
y <- rnorm(n, 53)

model <-t_test(x, y, paired = TRUE)
get_estimates(model)

hypothesis <- "mu = 0; mu < 0; mu > 0"

BF(model, hypothesis = hypothesis)

leads to the following error message:

> BF(model, hypothesis = hypothesis)
Error in parse_hypothesis(names_coef, hypothesis) : 
  Some of the parameters referred to in the 'hypothesis' do not correspond to parameter names of object 'x'.
  The following parameter names in the 'hypothesis' did not match any parameters in 'x': mu
  The parameters in object 'x' are named: difference

Changing it to difference works, but then I would expect that get_estimates(model) should also return difference instead of mu.

In the case of an unpaired t-test, I always seem to get non-conformable arrays if I specify the hypothesis myself. For example,

library(BFpack)

set.seed(4)

n <- 100
x <- rnorm(n, 50)
y <- rnorm(n, 53)

model <-t_test(x, y)
get_estimates(model)

hypothesis <- "difference > 0"
BF(model, hypothesis = hypothesis)

yields the following error:

> BF(model, hypothesis = hypothesis)
Error in as.matrix(RrO[[h]][, 1]) %*% t(drawsN) - as.matrix(RrO[[h]][,  : 
  non-conformable arrays
cjvanlissa commented 3 years ago

@jomulder this is also a bug in BFpack, because you manually rename the estimates to mu without considering the fact that they have different names. @tpourmohamad 's code works in bain:

library(bain)

set.seed(4)

n <- 100
x <- rnorm(n, 50)
y <- rnorm(n, 53)

model <-t_test(x, y, paired = TRUE)
get_estimates(model)

hypothesis <- "difference = 0; difference < 0; difference > 0"

bain(model, hypothesis = hypothesis)

model <-t_test(x, y)
get_estimates(model)

hypothesis <- "x > y"
bain(model, hypothesis = hypothesis)

As you can see, the names of the estimates of t_test() change depending on the model being evaluated. The current BF.t_test() does not take this into account, assuming everything is named mu.

> get_estimates(t_test(1:10, y = c(7:20)))
     x      y 
 5.500 13.500 
> get_estimates(t_test(extra ~ group, data = sleep))
group1 group2 
 0.750  2.330 
> get_estimates(t_test(sleep$extra))
    x 
1.540 
> sleep2 <- reshape(sleep, direction = "wide", 
+                   idvar = "ID", timevar = "group")
> get_estimates(t_test(sleep2$extra.1, sleep2$extra.2, paired = TRUE))
difference 
    -1.580 
jomulder commented 3 years ago

Hello,

Thanks for spotting and reporting these bugs! I corrected both bugs. The first was in the naming of the variable for a 2 sample t test when using the get_estimates function. The second was in the number of draws for the computation of BFs in the case of a t test with unequal variances. When installing the developmental version via

remotes::install_github("jomulder/BFpack")

these bugs should be avoided.

Thanks again.

Best, Joris