demo/sim2026_plotResults.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

resultsPlot<-function(R=Results,doPct=TRUE,legloc=1,parm="D"){
  
  h<-ifelse(doPct, 100,1)
  # for the lines, we number the same as for the 
  # uncensored plots
  # 1=BPCP  (b_)
  # 2=mid-p BPCP (bmp_)
  # 3=adj hybrid Z-O (dah_)
  # 4=standard Z-O (d_)
  #LWD<-c(13,12,11,6,2)
  LWD<-c(13,10,7,4)
  LWD<-c(1,1,1,1)
  LWD<-4:1
  LTY<-c(1,1,1,2)
  COL<- gray(c(.5,.7,.3,.9,.1)-.1)
  COL<-palette.colors(4)[c(1,3,2,4)]

  X0<-c(0.8,1,1.2,1.4,1.8)
  sep<- 0.4
  x<- c(X0,1.8+sep+X0,2*(1.8+sep)+X0,3*(1.8+sep)+X0)
  
  
  # calculate the number of t per Scen.
  nps<- length(X0)
  if (4*nps !=nrow(R)) stop("rewrite program, R show have nrow(R)=4*nps")
  
  i1<- 1:length(X0)
  i2<- i1+length(X0)
  i3<- i1+2*length(X0)
  i4<- i1+3*length(X0)
  
  dolines<-function(ii=i1,var="_power_gr",X=x,yfactor=1){
    lines(x=X[ii],y=yfactor*R[ii,paste0("b",var)],lwd=LWD[1],lty=LTY[1],col=COL[1])
    lines(x=X[ii],y=yfactor*R[ii,paste0("bmp",var)],lwd=LWD[2],lty=LTY[2],col=COL[2])
    lines(x=X[ii],y=yfactor*R[ii,paste0("dah",var)],lwd=LWD[3],lty=LTY[3],col=COL[3])
    lines(x=X[ii],y=yfactor*R[ii,paste0("d",var)],lwd=LWD[4],lty=LTY[4],col=COL[4])
  }
  
  ###################################
  # Power Plot
  ###################################
  
  plot(c(0,max(x)+0.8),h*c(0,1),type="n",xlab="",ylab="Power (%)",
       axes=FALSE, cex.lab=1,
       main=expression(paste("Power to Show ",S[2],"(t) > ",S[1],"(t)")))
  AT<-c(1.3,1.3+1.8+sep,1.3+2*(1.8+sep),1.3+3*(1.8+sep))
  axis(1,at=AT,labels=c("Scen. 1","Scen. 2","Scen. 3","Scen. 4"),
       cex.axis=0.8,tck=0)
  axis(2,las=2)
  box()
  
  dolines(i1,var="_power_gr",yfactor=h)
  dolines(i2,var="_power_gr",yfactor=h)
  dolines(i3,var="_power_gr",yfactor=h)
  dolines(i4,var="_power_gr",yfactor=h)
  
  o<-c(4:1)
  if (legloc==1){
  legend("topleft",legend=c("meld BPCP","meld BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)")[o],
         lwd=LWD[o],lty=LTY[o],col=COL[o],cex=0.7)
  }
  ###################################
  # CI Width Plot
  ###################################
  
  yRange<-range(R[,c("b_width","bmp_width","dah_width","d_width")])
  
  if (parm=="D"){
    TITLE<-"95% Confidence Interval Width"
    YLAB<-"95% CI Width"
  } else {
    TITLE<-"Transformed 95% CI Width"
    YLAB<-"Trans. 95% CI Width"
  }
  
  plot(c(0,max(x)+0.8),yRange,type="n",xlab="",ylab=YLAB,
       axes=FALSE,cex.lab=1,
       main=TITLE)
  AT<-c(1.3,1.3+1.8+sep,1.3+2*(1.8+sep),1.3+3*(1.8+sep))
  axis(1,at=AT,labels=c("Scen. 1","Scen. 2","Scen. 3","Scen. 4"),
       cex.axis=0.8,tck=0)
  axis(2,las=2)
  box()
  
  dolines(i1,var="_width")
  dolines(i2,var="_width")
  dolines(i3,var="_width")
  dolines(i4,var="_width")
  
  o<-c(4:1)
  if (legloc==2){
    if (parm=="D"){
      legend("bottomright",legend=c("meld BPCP","meld BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)")[o],
             lwd=LWD[o],lty=LTY[o],col=COL[o],cex=0.7)
    } else if (parm=="E"){
      legend("topright",legend=c("meld BPCP","meld BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)")[o],
             lwd=LWD[o],lty=LTY[o],col=COL[o],cex=0.7)
    }

  }
  
  ###################################
  # CI Error Lo
  ###################################
  
  yRange<- h*range(R[,c("b_err_lo","bmp_err_lo","dah_err_lo","d_err_lo")])
  #yRange<-c(0,4)
  
  plot(c(0,max(x)+0.8),yRange,type="n",xlab="",ylab="Percent Error",
       axes=FALSE,cex.lab=1,
       #main=expression(paste("Percent of Lower Confidence Limit>",beta))
       main="Pct. of Lower Conf. Limit> beta"
       )
  AT<-c(1.3,1.3+1.8+sep,1.3+2*(1.8+sep),1.3+3*(1.8+sep))
  axis(1,at=AT,labels=c("Scen. 1","Scen. 2","Scen. 3","Scen. 4"),
       cex.axis=0.8,tck=0)
  axis(2,las=2)
  box()
  lines(c(-10,10),c(2.5,2.5),lty=2,col="red")
  
  dolines(i1,var="_err_lo",yfactor=h)
  dolines(i2,var="_err_lo",yfactor=h)
  dolines(i3,var="_err_lo",yfactor=h)
  dolines(i4,var="_err_lo",yfactor=h)
  
  
  o<-c(4:1)
  if (legloc==3){
    legend("topleft",legend=c("meld BPCP","meld BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)")[o],
           lwd=LWD[o],lty=LTY[o],col=COL[o],cex=0.7)
  }
  ###################################
  # CI Error Hi
  ###################################
  
  yRange<- h*range(R[,c("b_err_hi","bmp_err_hi","dah_err_hi","d_err_hi")])
  #yRange<-c(0,4)
  plot(c(0,max(x)+0.8),yRange,type="n",xlab="",ylab="Percent Error",
       axes=FALSE,cex.lab=1,
       #main=expression(paste("Percent of Upper Confidence Limit<",beta))
       main="Pct. of Upper Conf. Limit<beta"
       )
  AT<-c(1.3,1.3+1.8+sep,1.3+2*(1.8+sep),1.3+3*(1.8+sep))
  axis(1,at=AT,labels=c("Scen. 1","Scen. 2","Scen. 3","Scen. 4"),
       cex.axis=0.8,tck=0)
  axis(2,las=2)
  box()
  lines(c(-10,10),c(2.5,2.5),lty=2,col="red")
  
  dolines(i1,var="_err_hi",yfactor=h)
  dolines(i2,var="_err_hi",yfactor=h)
  dolines(i3,var="_err_hi",yfactor=h)
  dolines(i4,var="_err_hi",yfactor=h)
  
  
  o<-c(4:1)
  if (legloc==4){
    legend("topright",legend=c("meld BPCP","meld BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)")[o],
           lwd=LWD[o],lty=LTY[o],col=COL[o],cex=0.7)
  }
}

#resultsPlot(R=Results)
#dev.print(pdf,file="./simBPCP/plotResults1.pdf")

########################################
# Difference simulation results
#######################################
load(file="./simBPCP/simResults/Results_difference_Run_1.RData")
pdf("./simBPCP/simResults/RD01.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()


load(file="./simBPCP/simResults/Results_difference_Run_2.RData")
pdf("./simBPCP/simResults/RD02.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()


load(file="./simBPCP/simResults/Results_difference_Run_3.RData")
pdf("./simBPCP/simResults/RD03.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_4.RData")
pdf("./simBPCP/simResults/RD04.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_5.RData")
pdf("./simBPCP/simResults/RD05.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_6.RData")
pdf("./simBPCP/simResults/RD06.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_7.RData")
pdf("./simBPCP/simResults/RD07.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_8.RData")
pdf("./simBPCP/simResults/RD08.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_9.RData")
pdf("./simBPCP/simResults/RD09.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_10.RData")
pdf("./simBPCP/simResults/RD10.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_11.RData")
pdf("./simBPCP/simResults/RD11.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_12.RData")
pdf("./simBPCP/simResults/RD12.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_13.RData")
pdf("./simBPCP/simResults/RD13.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=4)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_14.RData")
pdf("./simBPCP/simResults/RD14.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=4)
dev.off()


load(file="./simBPCP/simResults/Results_difference_Run_15.RData")
pdf("./simBPCP/simResults/RD15.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results, legloc=4)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_16.RData")
pdf("./simBPCP/simResults/RD16.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=4)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_17.RData")
pdf("./simBPCP/simResults/RD17.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=4)
dev.off()

load(file="./simBPCP/simResults/Results_difference_Run_18.RData")
pdf("./simBPCP/simResults/RD18.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=4)
dev.off()


########################################
# efflogs simulation results
#######################################
load(file="./simBPCP/simResults/Results_efflogs_Run_1.RData")
pdf("./simBPCP/simResults/RE01.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()


load(file="./simBPCP/simResults/Results_efflogs_Run_2.RData")
pdf("./simBPCP/simResults/RE02.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()


load(file="./simBPCP/simResults/Results_efflogs_Run_3.RData")
pdf("./simBPCP/simResults/RE03.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_4.RData")
pdf("./simBPCP/simResults/RE04.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_5.RData")
pdf("./simBPCP/simResults/RE05.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_6.RData")
pdf("./simBPCP/simResults/RE06.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results)
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_7.RData")
pdf("./simBPCP/simResults/RE07.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2, parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_8.RData")
pdf("./simBPCP/simResults/RE08.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_9.RData")
pdf("./simBPCP/simResults/RE09.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2, parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_10.RData")
pdf("./simBPCP/simResults/RE10.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_11.RData")
pdf("./simBPCP/simResults/RE11.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_12.RData")
pdf("./simBPCP/simResults/RE12.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_13.RData")
pdf("./simBPCP/simResults/RE13.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_14.RData")
pdf("./simBPCP/simResults/RE14.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()


load(file="./simBPCP/simResults/Results_efflogs_Run_15.RData")
pdf("./simBPCP/simResults/RE15.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_16.RData")
pdf("./simBPCP/simResults/RE16.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()

load(file="./simBPCP/simResults/Results_efflogs_Run_17.RData")
pdf("./simBPCP/simResults/RE17.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()


load(file="./simBPCP/simResults/Results_efflogs_Run_18.RData")
pdf("./simBPCP/simResults/RE18.pdf", width = 7, height = 5)
par(mfrow=c(2,2),mar=c(2,4,3,1)+0.1)
resultsPlot(R=Results,legloc=2,parm="E")
dev.off()