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 $H_0$ is simply a standard Gaussian distribution ($\mu=0$ and $\sigma=1$), while $H_1$ is still a Gaussian centered in 0, with a very narrow width ($\sigma=1/1000$). The prior odds are set at 1 to 1, i.e. $P(H_1)=P(H_0)=1/2$. 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.³¹ You can now play with the simulations, varying the parameters. If you want to get a situation yielding Bayes factors of ${\cal O}(10^{10})$ you can keep the standard parameters of $H_0$, fixing instead $\mbox{{\tt mu.H1}}$ at $1.7$ and $\mbox{{\tt sigma.H1}}$ at $\approx 4\times 10^{-10}$. 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().)