## 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))