Binomial.R
In StudyPrior: Distributions for Historical Information

### R code from vignette source 'Binomial.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: data
###################################################
library(StudyPrior)
x <- 21
n <- 37

#Full Bayes with delta~Be(1,1)
pp.fullbayes.11 <- binom.PP.FB.MC.BE(x = x, n = n, d.prior.a = 1, d.prior.b = 1) 
#Full Bayes with delta~Be(.5,.5)
pp.fullbayes.55 <- binom.PP.FB.MC.BE(x = x, n = n, d.prior.a = .5, d.prior.b = .5) 
#Power prior with fixed delta=0.8
pp.fix.08 <- binom.PP.FIX(x=x, n=n, d=0.8) 


###################################################
### code chunk number 2: EmpiricalBayes
###################################################
n.new <- 50
pp.eb <- binom.PP.EB(x=x, n=n, X=0:n.new, N=n.new)


###################################################
### code chunk number 3: plot
###################################################
plot(pp.fullbayes.11, col="red",
     xlim=c(0,1), ylim=c(0,6), ylab="Density", 
     xlab=expression(paste("Success Probability",~theta)))
curve(pp.fullbayes.55, add=TRUE, col="blue")
curve(pp.fix.08, add=TRUE, col="green")
curve(pp.eb(x,30), add=TRUE, col="orange")

legend("topleft", col=c("red","blue","green", "orange"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Fixed PP d=0.8", "Empirical Bayes PP"))


###################################################
### code chunk number 4: MSE
###################################################
MSE <- sapply(list(pp.fullbayes.11, pp.fullbayes.55, pp.fix.08, pp.eb),
              function(PRIOR) calc.MSE(prior=PRIOR, 
          prob.range = c(0.3,0.8), length = 11, n.binom=n.new))

p <- seq(0.3,0.8,len=11)
plot(p,MSE[,1],  col="red" , pch=16, type='l',
     ylab="MSE", xlab=expression(paste("Success Probability ",theta)),
     ylim=c(0.005,0.015))
lines(p,MSE[,2],  col="blue", pch=16)
lines(p,MSE[,3],  col="green", pch=16)
lines(p,MSE[,4],  col="orange", pch=16)

legend("top", col=c("red","blue","green", "orange"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Fixed PP d=0.8", "Empirical Bayes PP"))



###################################################
### code chunk number 5: sigmat
###################################################

SM <- lapply(list(pp.fullbayes.11, pp.fullbayes.55, pp.fix.08, pp.eb),
       function(PRIOR) sig.matrix(n.new, n.new, 0.95, prior=PRIOR, mc.cores = 1))


###################################################
### code chunk number 6: Power
###################################################
Difference <- 0.2

POWER <- sapply(SM, function(SIGMAT) calc.power(sig.mat=SIGMAT, 
      prob.rang=c(0.3,0.8), length=11, treatment.difference = Difference,
      n.binom.control = n.new, n.binom.treatment = n.new))

plot(p,POWER[,1],  col="red", type='l', ylab="Power",
     xlab=expression(paste("Success Probability ",theta)))
lines(p,POWER[,2],  col="blue")
lines(p,POWER[,3],  col="green")
lines(p,POWER[,4],  col="orange")

legend("topleft", col=c("red","blue","green", "orange"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Fixed PP d=0.8", "Empirical Bayes PP"))


###################################################
### code chunk number 7: TypeIError
###################################################
T1E <- sapply(SM, function(SIGMAT) calc.power(sig.mat=SIGMAT, 
         prob.rang=c(0.3,0.8), length=11, treatment.difference = 0,
         n.binom.control = n.new, n.binom.treatment = n.new))


plot(p,T1E[,1],  col="red", type='l', ylab="Type I Error",
     xlab=expression(paste("Success Probability ",theta)),
     ylim=c(0,0.15))
lines(p,T1E[,2],  col="blue")
lines(p,T1E[,3],  col="green")
lines(p,T1E[,4],  col="orange")


legend("topleft", col=c("red","blue","green", "orange"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Fixed PP d=0.8", "Empirical Bayes PP"))



###################################################
### code chunk number 8: table
###################################################
historical <- data.frame(Study=1:4, x=c(21,33,23,27), n=c(37,50,40,50), 'Success'= paste0(round(c(21,33,23,27)/c(37,50,40,50)*100),'%'))

print(xtable::xtable(historical, digits=c(0,0,0,0,2)), include.rownames = FALSE)


###################################################
### code chunk number 9: MAP
###################################################
x.multi <- c(21,33,23,27)
n.multi <- c(37,50,40,50)

MAP <- binom.MAP.FB(x.multi, n.multi, tau.prior = list(prior= "logtnormal", param=c(0,1)))


###################################################
### code chunk number 10: PPFB
###################################################
PP.FB.11 <- binom.PP.FB.MC.BE(x=x.multi, n=n.multi, d.prior.a = 1, d.prior.b = 1)
PP.FB.55 <- binom.PP.FB.MC.BE(x=x.multi, n=n.multi, d.prior.a = .5, d.prior.b = .5)


###################################################
### code chunk number 11: PPCOR
###################################################
PP.FB.COR <- binom.PP.FB.MC.COR(x=x.multi, n=n.multi, d.prior.cor = 0.5)


###################################################
### code chunk number 12: PPEB
###################################################
PP.EB <- binom.PP.EB(x=x.multi, n=n.multi, X=0:n.new, N=n.new)


###################################################
### code chunk number 13: approx
###################################################
MAP.approx <- conj.approx(MAP, type="beta", max.degree = 3)

curve(MAP, n=200, lty=1, lwd=2, ylab="Density", xlab=expression(theta))
plot(MAP.approx,  col="cyan", lwd=2, lty=2, lines.only=TRUE)

legend("topleft", lty=1:2, col=c("black","cyan"), lwd=2, legend=c("MAP","Approximation"))


###################################################
### code chunk number 14: bias
###################################################
BIAS.PP <- sapply(list(PP.FB.11, PP.FB.55, PP.FB.COR, PP.EB), 
                  function(PRIOR) calc.bias(PRIOR,
                                            prob.range = c(0.3,0.8),
                                            length=11,
                                            n.binom=50)
)


###################################################
### code chunk number 15: biasmap
###################################################
system.time(BIAS.MAP <- calc.bias(MAP, prob.range = c(0.3,0.8),length=11,n.binom = 50))
system.time(BIAS.APPROX <- calc.bias(MAP.approx, prob.range = c(0.3,0.8),length=11,n.binom = 50))


###################################################
### code chunk number 16: biasplot
###################################################
plot(p,BIAS.PP[,1],  col="red", type='l', ylab="Bias",
     xlab=expression(paste("Success Probability ",theta)),
     ylim=c(-0.15,0.15), lwd=2)
lines(p,BIAS.PP[,2],  col="blue", lwd=2)
lines(p,BIAS.PP[,3],  col="green", lwd=2)
lines(p,BIAS.PP[,4],  col="orange", lwd=2)

lines(p,BIAS.MAP,  col="black", lwd=2)
lines(p,BIAS.APPROX,  col="cyan", lty=2, lwd=2)

legend("topright", col=c("red","blue","green", "orange", "black","cyan"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Full Bayes PP with Corr=.5", "Empirical Bayes PP","MAP","MAP Approx"))



###################################################
### code chunk number 17: cover
###################################################
COVER <- sapply(list(PP.FB.11, PP.FB.55, PP.FB.COR, PP.EB, MAP), 
                  function(PRIOR) calc.coverage(PRIOR, level = 0.95, 
                                                smooth = 0.05, n.control = 50))
p2 <- seq(0.3,.8, len=501)

plot(p2,COVER[300:800,1],  col="red", type='l', ylab="Coverage",
     xlab=expression(paste("Success Probability ",theta)),
     ylim=c(0,1), lwd=2)

lines(p2,COVER[300:800,2],  col="blue", lwd=2)
lines(p2,COVER[300:800,3],  col="green", lwd=2)
lines(p2,COVER[300:800,4],  col="orange", lwd=2)
lines(p2,COVER[300:800,5],  col="black", lwd=2)

legend("bottom", col=c("red","blue","green", "orange", "black"), lty=1, lwd=2,
       legend=c("Full Bayes PP B(1,1)", "Full Bayes PP B(.5,.5)",
                "Full Bayes PP with Corr=.5", "Empirical Bayes PP","MAP"))