FAQ/powBayest

(taken from )

A power calculation is an important step in planning any experiment. In particular, virtually all clinical trials include a power analysis to determine the number of patients to enroll in the study. The power of the test is the (frequentist) probability of rejecting the null hypothesis of no treatment effect when there is in fact a treatment effect. Computing the power thus requires specifying a treatment effect (or at least a prior distribution on the treatment effect).

This is an example of using Monte Carlo simulation to find the power of a clinical trial. The statistical model is

Y 1 ,…,Yn ∼N(μ,σ 2 ) Y1,…,Yn∼N(μ,σ2)

and the hypotheses are H O :μ=0 versus H A :μ≠0. HO:μ=0 versus HA:μ≠0. In this example we assume the variance σ 2 σ2 is known, and use an uninformative prior for μ μ.

Simulation Settings

           <- 100  # Number of MC reps
n           <- 100  # Sample size
sigma       <- 1    # Error standard devation

true_mu     <- 0.2  # True value of mu

pri.mn.anal <- 0    # Prior for mu
pri.sd.anal <- Inf

Generate data sets and credible regions for posterior mean

L <- rep(0,N)
 U <- rep(0,N)

 set.seed(0820)

 for(rep in 1:N){

  #Generate data:
    Y  <- rnorm(n,true_mu,sigma)

   #Compute posterior 95% interval:
    post.var <- 1/(n/sigma^2+1/pri.sd.anal^2)
    post.mn  <- post.var*(pri.mn.anal/pri.sd.anal^2+sum(Y)/sigma^2)
    L[rep]   <- post.mn-1.96*sqrt(post.var)
    U[rep]   <- post.mn+1.96*sqrt(post.var)
 }

Plot the credible set for each dataset

plot(NA,xlim=c(-1,1),ylim=c(1,N),xlab="95% interval",ylab="Replication")
abline(v=0)

for(rep in 1:N){
    reject <- L[rep]>0 | U[rep]<0
    lines(c(L[rep],U[rep]),c(rep,rep),col=ifelse(reject>0,2,1))
}

Compute the power which is proportion of times 0 (the true mean) is outside the posterior credible region.

mean(L>0 | U<0)

Output

## [1] 0.53

Bayesian power calculation for a one sample t-test