Simulate p-value fallacy

pvalue_fallacy.r

# simulate pvalue fallacy
# See https://stat.duke.edu/courses/Spring10/sta122/Labs/Lab6.pdf
# and https://stat.duke.edu/~berger/applet2/applet.pdf
sim.pval = function(pi = 0.5, theta.input = 1, sigma.input = 1, n = 100, max = 1000){
  alpha = 0
  beta = 0
  while(alpha + beta < max){
    rand = runif(1)
    if(rand > pi){
      theta = 0
      sigma = 1
      null.true = TRUE
    } else {
      theta = theta.input
      sigma = sigma.input
      null.true = FALSE
    }
    samp = rnorm(n,theta,sigma)
    t = sqrt(n) * abs(mean(samp))/sigma
    p.value = 2 * (1 - pnorm(t))
    if(p.value < 0.05 & p.value > 0.049){
      if(null.true){
        alpha = alpha + 1
      } else {
        beta = beta + 1
      }
    }
  }  
  return(list(alpha=alpha,beta=beta))  
}


##### Do simulation
y = sim.pval(theta=2,sigma=3,n=20,max=1000); 

# Graph stuff
y.mean = y$alpha/(y$alpha+y$beta)
y.sd = sqrt((y$alpha*y$beta) / ((y$alpha + y$beta)^2 * (y$alpha + y$beta + 1)))
curve(dbeta(x,y$alpha,y$beta),y.mean-y.sd*5,y.mean+y.sd*5,main=paste("Estimated type 1 error rate =",y.mean))