#########################################################
#  R code to implement Basic Peter/Telba N/N Model
#    BC, 10/20/08
#  This one adds the loop over true beta1/mubeta, 12/3/08
#   and the loop over priors for beta1, 12/16/08
#  This is for Task 2 (test for trt effect in current trial)  
#
# This is the BRugs code that calls three other files:
#   WinBUGS code:  NN_BUGS.txt
#   data:          NN_data.txt (created by this program)
#   inits:         NN_inits.txt (created by this program)
#########################################################

#setwd("K:/brianho/Telba")    # change this to *your* working directory
library(BRugs) 
set.seed(400)

n <- 20;  nrep <- 200

#  SPECIFY THE TRUE CONTROL AND TRT VALUES IN THE HISTORICAL AND CURRENT GROUPS:
truetheta0 <- 0;  truetheta1 <- 30
truebeta0 <- 0
truebeta1 <- seq(0,2,by=.25)

#  SPECIFY THE HYPERPRIORS THAT CONTROL THE DEGREE OF BORROWING FOR THETA:
#atheta <- .1;  btheta <- .1    # not much borrowing
atheta <- 3;  btheta <- 1      # a little borrowing
#atheta <- 1000;  btheta <- 10  # lots of borrowing

#  SPECIFY THE HYPERPRIORS THAT CONTROL THE DEGREE OF BORROWING FOR BETA:
abeta <-c(0,0,0);  bbeta<-c(0,0,0)
abeta[1] <-.1;  bbeta[1] <- .1    # not much borrowing
abeta[2] <-40;  bbeta[2] <- 4     # a little borrowing
abeta[3] <-1000;  bbeta[3] <- 10  # lots of borrowing

powerbeta1 <- matrix(0,length(abeta),length(truebeta1))  
dput(pairlist(mutheta=0, mubeta=0, etatheta=1, etabeta=1),"NN_inits.txt")

x0 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
x1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
#Y0 <- c( 0.38, -0.89, -0.19,  0.16,  1.24,  1.95, -0.07, -0.78,  0.24,  0.03, -0.57,
#                -0.58,  0.01, -0.57, -0.54, -0.73,  0.27, -0.05,  0.60, -0.70)
#Y1 <- c(10.13, 11.19, 10.15,  9.46,  9.47,  9.42, 10.39,  9.39,  9.86,  8.67, 10.89,
#               10.45, 10.97,  9.06, 10.38, 10.10, 10.00,  9.64,  9.98, 10.46)

# HERE IS OUTERMOST LOOP OVER the J=3 PRIOR CHOICES:

loopstart <- Sys.time() 
for (j in 1:length(abeta)){

# HERE IS LOOP OVER TRUE BETA1 VALUES (x-axis of power curve):

for (k in 1:length(truebeta1)){
#  rejectmubeta <- 0;  mubetastat <- NULL;  
  rejectbeta1 <- 0;   beta1stat <- NULL

# HERE IS MAIN LOOP OVER FAKE DATASETS: 

for (i in 1:nrep) {
Y0 <- round(rnorm(n,truetheta0+truebeta0*x0, 1),3)
Y1 <- round(rnorm(n,truetheta1+truebeta1[k]*x1, 1),3)

dput(pairlist(x0=x0,x1=x1,Y0=Y0,Y1=Y1,n=n,atheta=atheta,btheta=btheta,
  abeta=abeta[j],bbeta=bbeta[j]),"NN_data.txt")
modelCheck("NN_BUGS.txt")
modelData(paste("NN_data.txt",sep=""))
modelCompile(numChains=1)
modelInits("NN_inits.txt")
modelGenInits()    # to generate starting values for the theta's and beta's

modelUpdate(100)
#samplesSet(c("theta0","theta1","beta0","beta1","mutheta","mubeta",
#  "etatheta","etabeta","tausqtheta","tausqbeta"))
samplesSet(c("beta1"))

modelUpdate(1000)
#samplesHistory("*", mfrow = c(3,2))
#samplesAutoC("*", chain=1, mfrow = c(3,2))
#samplesDensity("*", mfrow = c(3,2))
#samplesStats(c("theta0","theta1","beta0","beta1","mutheta","mubeta",
#  "etatheta","etabeta","tausqtheta","tausqbeta"))

#  DECISION RULE FOR OVERALL TRT EFFECT:
#LL <- samplesStats("mubeta")$val2.5pc 
#UL <- samplesStats("mubeta")$val97.5pc 
#if (LL > 0) rejectmubeta <- rejectmubeta + 1 
#  else 
#if (UL < 0) rejectmubeta <- rejectmubeta + 1

#  DECISION RULE FOR TRT EFFECT IN CURRENT TRIAL:
LL <- samplesStats("beta1")$val2.5pc 
UL <- samplesStats("beta1")$val97.5pc 
if (LL > 0) rejectbeta1 <- rejectbeta1 + 1 
  else 
if (UL < 0) rejectbeta1 <- rejectbeta1 + 1

# Bind the summary statistics of the current iteration to beta1stat: 
#mubetastat <- rbind(mubetastat, samplesStats("mubeta")) 
beta1stat <- rbind(beta1stat, samplesStats("beta1")) 

} # end of Nrep loop

# Output final results: 
#dput(mubetastat,'mubetastat.txt', control=NULL) 

# Write power summaries to the screen: 
#cat("\nstats for mubeta:\n"); print(mubetastat)
cat("\nstats for beta_1:\n"); print(beta1stat)

cat("\n\nHere are simulated results for nrep=",nrep,", prior no.",j,", beta1 no.",k,
  ",\n  true theta0, theta1, beta0, beta1 are:  ",truetheta0,truetheta1,truebeta0,truebeta1[k],
  ",\n  and true atheta, btheta, abeta, bbeta are:  ",atheta,btheta,abeta[j],bbeta[j],"\n")

#powermubeta[j,k] <- rejectmubeta/nrep
powerbeta1[j,k] <- rejectbeta1/nrep

#cat("\nPr(reject H0: mu_beta=0) is:",powermubeta[k]) 
cat("\nPr(reject H0: beta_1=0) is: ",powerbeta1[j,k]) 

} # end of k loop (truebeta1's)
} # end of j loop (priors)

tottime <- Sys.time()-loopstart 
cat("\n\nTotal time elapsed: "); print(tottime) 
cat("\n\nEnd of Peter/Telba N/N power simulation\n\n") 

##############
#  now plot the power curve:
##############
postscript(paste("powerbeta1_",truetheta0,"_",truetheta1,"_",truebeta0,".ps",sep=""))
#pdf(paste("powerbeta1_",truetheta0,"_",truetheta1,"_",truebeta0,".pdf",sep=""),paper="letter")

par(lwd=1.5,cex=1.5)
#mar=c(1,1,1,1))
par(mfrow=c(1,1));  smooth <- .465
plot(lowess(truebeta1,powerbeta1[1,],f=smooth),type="l",ylab="power",ylim=c(0,1),
  xlab="true beta_1",cex=1.5)
#lines(truebeta1,powerbeta1[2,],lty=2)
lines(lowess(truebeta1,powerbeta1[2,],f=smooth),lty=2)
#lines(truebeta1,powerbeta1[3,],lty=3)
lines(lowess(truebeta1,powerbeta1[3,],f=smooth),lty=3)
abline(h=0)      #;  abline(h=1)

legend("topleft",cex=1.5,legend=c(paste("low shrink, G(",abeta[1],",",bbeta[1],")"),
paste("mid shrink, G(",abeta[2],",",bbeta[2],")"),
paste("high shrink, G(",abeta[3],",",bbeta[3],")")),lty=1:3)
#mtext("Power curve for trt effect in current trial (beta_1)",line=2,cex=1.2)
#mtext(paste("2 studies x 2 arms;  nrep=",nrep,";  true theta0, theta1, and beta0 are:  ",
#  truetheta0,truetheta1,truebeta0),line=1)

dev.off()
# 
#  END OF NN_beta1_R PROGRAM
#