Bayesian power calculation for a one sample t-test
(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