## Sample size for two binomial proportions. Target parameter is relative risk. ##### Predictive probability (Table 5.7) library(R2WinBUGS) SampleSizeProportions <- function(m,n1,n2,theta1T,theta2T,TrueRR,alpha){ power1 <- rep(NA,m) coverage1 <- rep(NA,m) length1 <- rep(NA,m) for(i in 1:m) { y1 <- rbinom(1,n1,theta1T) y2 <- rbinom(1,n2,theta2T) binomialdata <- list("n1","n2","y1","y2") parameters <- c("RR", "HT") inits <- list( list(theta1=0.5,theta2=0.5)) SampleSize <- bugs(binomialdata, inits, parameters,"SampleSizeProportions.txt", working.directory = "F:\\MyDocuments\\BAYES\\OneandTwoSample\\Code", n.chains=1, n.iter=2000, n.thin=1, n.burnin=0, save.history=FALSE) power1[i] <- SampleSize\$summary[2,1] length1[i] <- SampleSize\$summary[1,7] - SampleSize\$summary[1,3] coverage1[i]<-(TrueRR>=SampleSize\$summary[1,3])&(TrueRR<=SampleSize\$summary[1,7]) } AveragePower <- mean(power1) PredProb <- c( mean(power1>alpha[1]), mean(power1>alpha[2]), mean(power1>alpha[3]) ) AverageLength <- mean(length1) AverageCoverage <- mean(coverage1) return(c(PredProb, AveragePower,AverageLength,AverageCoverage)) } # RR=1.25 SampleSizeProportions(m=1000, n1=50,n2=50, theta1T=0.25,theta2T=0.2,TrueRR=1.25, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=100,n2=100, theta1T=0.25,theta2T=0.2,TrueRR=1.25, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=150,n2=150, theta1T=0.25,theta2T=0.2,TrueRR=1.25, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=200,n2=200, theta1T=0.25,theta2T=0.2,TrueRR=1.25, alpha=c(0.80,0.90,0.95)) ##### Average Power library(R2WinBUGS) SampleSizeProportions <- function(m,n1,n2,a1,b1,a2,b2,alpha){ power1 <- rep(NA,m) coverage1 <- rep(NA,m) length1 <- rep(NA,m) for(i in 1:m) { theta1 <- rbeta(1,a1,b1) theta2 <- rbeta(1,a2,b2) y1 <- rbinom(1,n1,theta1) y2 <- rbinom(1,n2,theta2) TrueRR <- theta1/theta2 binomialdata <- list("n1","n2","y1","y2") parameters <- c("RR", "HT") inits <- list( list(theta1=0.5,theta2=0.5)) SampleSize <- bugs(binomialdata, inits, parameters,"SampleSizeProportions.txt", working.directory = "F:\\MyDocuments\\BAYES\\OneandTwoSample\\Code", n.chains=1, n.iter=2000, n.thin=1, n.burnin=0, save.history=FALSE) power1[i] <- SampleSize\$summary[2,1] #power1[i] <- mean(SampleSize\$sims.list\$RR > nullvalue) length1[i] <- SampleSize\$summary[1,7] - SampleSize\$summary[1,3] coverage1[i]<-(TrueRR>=SampleSize\$summary[1,3])&(TrueRR<=SampleSize\$summary[1,7]) } AveragePower <- mean(power1) PredProb <- c( mean(power1>alpha[1]), mean(power1>alpha[2]), mean(power1>alpha[3]) ) AverageLength <- mean(length1) AverageCoverage <- mean(coverage1) return(c(AveragePower,PredProb, AverageLength,AverageCoverage)) } # RR=1.75 SampleSizeProportions(m=1000, n1=50,n2=50, a2=23.36,b2=90.45, a1=8.50,b1=14.92, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=100,n2=100, a2=23.36,b2=90.45, a1=8.50,b1=14.92, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=150,n2=150, a2=23.36,b2=90.45, a1=8.50,b1=14.92, alpha=c(0.80,0.90,0.95)) SampleSizeProportions(m=1000, n1=200,n2=200, a2=23.36,b2=90.45, a1=8.50,b1=14.92, alpha=c(0.80,0.90,0.95))