######################################################### # 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 over the priors for beta1, 12/16/08 # This is also for Task 3 ("to pool or not to pool") # # This is the BRugs code that calls three other files: # WinBUGS code: NN_BUGSdelta.txt # data: NN_datadelta.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; cee <- 1; Koverage <- .80 # SPECIFY THE TRUE CONTROL AND TRT VALUES IN THE HISTORICAL AND CURRENT GROUPS: truetheta0 <- 0; truetheta1 <- 0 #truebeta0 <- 0 truebeta0 <- 0 # or 2 #truebeta1 <- seq(0,2,by=.2); truebeta1 <- c(seq(0,5,by=1),10,17,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; tau = .316 abeta[3] <-1000; bbeta[3] <- 10 # lots of borrowing; tau = .1 powerdelta <- 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)){ rejectdelta <- 0; deltastat <- 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,cee=cee,atheta=atheta,btheta=btheta, abeta=abeta[j],bbeta=bbeta[j]),"NN_datadelta.txt") modelCheck("NN_BUGSdelta.txt") modelData(paste("NN_datadelta.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("withinc")) modelUpdate(1000) #samplesHistory("*", mfrow = c(3, 2)) #samplesAutoC("*", chain=1, mfrow = c(3, 2)) #samplesDensity("*", mfrow = c(3, 2)) #samplesStats(c("beta0,"beta1","delta","withinc")) # DECISION RULE FOR TRT EFFECT IN CURRENT TRIAL: if (samplesStats("withinc")$mean < Koverage) rejectdelta <- rejectdelta + 1 # Bind the summary statistics of the current iteration to beta1stat: deltastat <- rbind(deltastat, samplesStats("withinc")) } # end of Nrep loop # Write power summaries to the screen: cat("\nstats for withinc:\n"); print(deltastat) 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") powerdelta[j,k] <- rejectdelta/nrep cat("\nPr(reject H0: delta = beta1 - beta0 = 0) is: ",powerdelta[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("powerdelta_",truetheta0,"_",truetheta1,"_",truebeta0,".ps",sep="")) #pdf(paste("powerdelta_",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 <- .1 plot(lowess(truebeta1,powerdelta[1,],f=smooth),type="l",ylab="power",ylim=c(0,1), xlab="true beta_1",cex=1.5) #plot(truebeta1,powerdelta[1,],type="l",ylab="power",ylim=c(0,1), # xlab="true beta_1 value") lines(lowess(truebeta1,powerdelta[2,],f=smooth),lty=2) #lines(truebeta1,powerdelta[2,],lty=2) lines(lowess(truebeta1,powerdelta[3,],f=smooth),lty=3) #lines(truebeta1,powerdelta[3,],lty=3) #abline(h=0) #; abline(h=1) legend("center",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 to reject pooling based on delta = beta1 - beta0",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_delta_R PROGRAM # ###########early results for 0,0,0: #> truebeta1 # [1] 0 1 2 3 4 5 6 7 8 10 20 50 #> powerdelta # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] #[1,] 0.16 0.625 0.995 1 1 1 1 1 1 1 1 1 #[2,] 0.08 0.595 0.980 1 1 1 1 1 1 1 1 1 #[3,] 0.00 0.000 0.000 0 0 0 0 0 0 0 1 1 ###########early results for 0,0,2: #> truebeta1 # [1] 0 1 2 3 4 5 6 7 8 10 20 50 #> powerdelta # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] #[1,] 0.995 0.615 0.165 0.74 0.985 1 1 1 1 1 1.000 1 #[2,] 0.990 0.670 0.145 0.49 0.960 1 1 1 1 1 1.000 1 #[3,] 0.000 0.000 0.000 0.00 0.000 0 0 0 0 0 0.325 1