# Power_BRugs runs a Bayesian simulation for calculating the power
# of a clinical trial using weibull regression
#
#  This is the joe version of "Power.R", modified by BPC 12/11/05

# NOTE: Don't forget to change line 2 to your own working directory:

oldwd<-getwd() #save your old working directory
setwd("K:\\book\\CL3\\BRugs\\PowerCalcs") #change your working directory

					
library(BRugs)
set.seed(400)

N <- 30						# N: sample size per group
nrep <- 5        				# nrep: number of iteratio
X <- rep(0:1, c(N, N))
beta1stat <- NULL

medsurv <- 36						# median survival time
A <- 2							# shape parameter
imprv <- 1.5						# improvement for treatment group
censmean <- 80						# mean censoring time
censsd<- 10							# sd of censoring time

hypbeta0 <- log(medsurv^A/log(2))
hypbeta1 <- log(imprv^A)

acccontrol <- 0; acctrt <- 0;  rejcontrol <-0;  rejtrt <- 0
equiv <- 0;  nodec <- 0

#pairlist and dput outputs initial values or data in WinBUGS format

initialize <- pairlist(beta0 = hypbeta0, beta1 = hypbeta1)	
dput(initialize, "powerinits.txt", control=NULL)

#  HERE IS MAIN LOOP OVER FAKE DATASETS:

loopstart <- Sys.time()
for (i in 1:nrep) {

#Sampling true parameters from priors

truebeta0 <- rnorm(1, mean = hypbeta0, sd = .2) 	# beta0 prior
truebeta1 <- rnorm(1, mean = hypbeta1, sd = .17) 	# optimistic prior
#truebeta1 <- rnorm(1, mean = hypbeta1/2, sd = .05) 	# equivalence prior
#truebeta1 <- rnorm(1, sd = .17) 			# skeptical prior 

surv <- rep(0,2*N)					
cens <- rep(0,2*N)					
T <- rep(0,2*N)						
T.cens <- rep(0,2*N)					
yinits <- rep(0, 2*N)

# The following steps translate beta1 and beta2 into the scale parameter
# for the weibull distribution

b1 <- sqrt(1/exp(-truebeta0))				 
b2 <- sqrt(1/exp(-(truebeta0 + truebeta1)))		 

surv[1:N] <- rweibull(N, shape = A, scale = b1)		# Generating survival times
surv[N+1:(2*N)] <- rweibull(N, shape = A, scale = b2)	# for the two groups
cens <- rnorm(2*N, mean = censmean, sd = censsd)	# Generating censoring times

# The 'j' loop places values in T and T.cens
# T: survival times for individuals with events and NA for censored
# T.cens: censoring times for individuals who were censored and NA for events

length <- 2*N
for(j in 1:length) {

if (surv[j] > cens[j]){
yinits[j] <- surv[j]
T[j] <- NA
T.cens[j] <- max(0, cens[j])
}

if (surv[j] <= cens[j]) {
T[j] <- surv[j]
yinits[j] <- NA
}

}

#pairlist and dput output the current data set in WinBUGS format

mydata <- pairlist(t = T, t.cens = T.cens, x = X, n = 2*N)
dput(mydata, "powerdata.txt", control=NULL)

#The next several steps are BRugs commands

#modelCheck: checks the model file, written in usual WinBUGS syntax

#modelCheck("optimistmodel.txt")
#modelCheck("skepticmodel.txt")
modelCheck("refmodel.txt")

#modelData: reads in the data file, written in the usual WinBUGS syntax

modelData("powerdata.txt")

#modelCompile: compiles the model

modelCompile()

#modelInits: enters initial values, in the usual WinBUGS format
#modelGenInits: generates initial values

modelInits("powerinits.txt")
modelGenInits()

#update, set, and stats are more BRugs commands

modelUpdate(100)
samplesSet("beta0")
samplesSet("beta1")
dicSet()

modelUpdate(500)
samplesAutoC("*",chain=1)
samplesHistory("*")


dicStats()

modelUpdate(500)
samplesDensity("*")

LL <- samplesStats("beta1")$val2.5pc
UL <- samplesStats("beta1")$val97.5pc 

if (UL < 0) acccontrol <- acccontrol + 1
  else
if (LL > hypbeta1) acctrt <- acctrt + 1
  else
if ((UL < hypbeta1) & (LL > 0))  equiv <- equiv + 1
  else
if ((UL < hypbeta1) & (LL < 0))  rejtrt <- rejtrt + 1
  else
if ((UL > hypbeta1) & (LL > 0))  rejcontrol <- rejcontrol + 1
  else
nodec <- nodec + 1

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

}   #  end of Nrep loop

# Output final results:
dput(beta1stat,'100skepref.txt', control=NULL)

# Plot the true survival and censoring densities:
time <- seq(from=0,to=200,length=100)
survdens <- dweibull(time, shape = A, scale = b1)
censdens <- dnorm(time, mean = censmean, sd = censsd)
plot(time,survdens,type="l",ylab="density",ylim=c(0,max(survdens,censdens)))
lines(time,censdens,lty=2)
# abline(v=medsurv)
title("True survival (solid) and censoring (dotted) distributions")

# Write power summaries to the screen:

cat("\n\nHere are simulated outcome frequencies for N=",N,", nrep=",nrep,":\n")

cat("\naccept control:   ",acccontrol/nrep)
cat("\nreject treatment: ",rejtrt/nrep)
cat("\nequivalence:      ",equiv/nrep)
cat("\nreject control:   ",rejcontrol/nrep)
cat("\naccept treatment: ",acctrt/nrep)
cat("\nno decision:      ",nodec/nrep)

tottime <- Sys.time()-loopstart
cat("\n\nTotal time elapsed: ");  print(tottime)

setwd(oldwd) #set working directory back to the original one

cat("\n\nEnd of BRugs power simulation\n\n")