## Playing with simulations

I hope it is now clear the reason why p-values and Bayes factors have in principle nothing to do with each other, and why p-values are not only responsible of unjustified claims of discoveries, but might also relegate genuine signals to the level of fluke, or reduce their significance', the word now used as normally understood and not with the technical meaning' of statisticians. But since I know that many might not be used with the reasoning just shown, I made a little R script[56], so that those who are still sceptical can run it and get a feeling of what is going on.
# initialization
mu.H0 <- 0; sigma.H0 <- 1
mu.H1 <- 0; sigma.H1 <- 1e-3
p.H1  <- 1/2
mu    <- c(mu.H0, mu.H1)
sigma <- c(sigma.H0, sigma.H1)

# simulation function
simulate <- function() {
M <- rbinom(1, 1, p.H1); x <- rnorm(1, mu[M+1], sigma[M+1])
x <- rnorm(1, mu[M+1], sigma[M+1])
p.val <- 2 * pnorm(mu[1] - abs(x-mu[1]), mu[1], sigma[1])
BF <- dnorm(x, mu[2], sigma[2]) / dnorm(x, mu[1], sigma[1])
lBF <- dnorm(x, mu[2], sigma[2], log=TRUE) - dnorm(x, mu[1], sigma[1], log=TRUE)
cat(sprintf("x = %.5f  =>  p.val = %.2e,  BF = %.2e  [ log(BF) = %.2e ]\n",
x, p.val, BF, lBF))
return(M)
}

By default is simply a standard Gaussian distribution ( and ), while is still a Gaussian centered in 0, with a very narrow width (). The prior odds are set at 1 to 1, i.e. . Each call to the function simulate() prints the values that we would get in a real experiment (x, p-value, Bayes factor and its log) and returns the true model (0 or 1), stored in a vector variable for later check. In this way you can try to infer what was the real cause of x before knowing the truth' (in simulations we can, in physics we cannot!). Here are the results of a small run, with x = 12 chosen in order to fill the page, thus postponing the solution to the next one.

> set.seed(150914); n=12; M <- rep(NA, n); for(i in 1:n) M[i] <- simulate()
x = -0.00079  =>  p.val = 9.99e-01,  BF = 7.29e+02  [ log(BF) = 6.59e+00 ]
x = -0.62293  =>  p.val = 5.33e-01,  BF = 0.00e+00  [ log(BF) = -1.94e+05 ]
x = -0.00029  =>  p.val = 1.00e+00,  BF = 9.57e+02  [ log(BF) = 6.86e+00 ]
x = -0.00162  =>  p.val = 9.99e-01,  BF = 2.68e+02  [ log(BF) = 5.59e+00 ]
x = -0.39258  =>  p.val = 6.95e-01,  BF = 0.00e+00  [ log(BF) = -7.71e+04 ]
x = -0.82578  =>  p.val = 4.09e-01,  BF = 0.00e+00  [ log(BF) = -3.41e+05 ]
x = 0.00073  =>  p.val = 9.99e-01,  BF = 7.69e+02  [ log(BF) = 6.64e+00 ]
x = -0.00012  =>  p.val = 1.00e+00,  BF = 9.93e+02  [ log(BF) = 6.90e+00 ]
x = 0.22295  =>  p.val = 8.24e-01,  BF = 0.00e+00  [ log(BF) = -2.48e+04 ]
x = -0.00022  =>  p.val = 1.00e+00,  BF = 9.76e+02  [ log(BF) = 6.88e+00 ]
x = 0.00117  =>  p.val = 9.99e-01,  BF = 5.07e+02  [ log(BF) = 6.23e+00 ]
x = -1.03815  =>  p.val = 2.99e-01,  BF = 0.00e+00  [ log(BF) = -5.39e+05 ]

And the winners are:

> M
[1] 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 0 0 0 1 0 1 1

It should not be any longer a surprise that the best figure to discriminate between the two models is the Bayes factor and not the p-value.31 You can now play with the simulations, varying the parameters. If you want to get a situation yielding Bayes factors of you can keep the standard parameters of , fixing instead at and at . Then you can choose p.H1 at wish and run the simulation. (You also need to change the numbers of digits of x, replacing %.5f'' by `%.11f'' inside sprintf().)

Giulio D'Agostini 2016-09-06