R/inst/app/functions/users/nzi_source.r

In BEACH: Biometric Exploratory Analysis Creation House

if(TRUE){#header
  #*****************************************************************************
  #CODE NAME (required)                : nzi_source.r
  #PROJECT NAME (required)             : Near-Zero-Informed Bayesian (NZIB)
  #PLATFORM                            : x86_64-w64-mingw32 
  #SOFTWARE/VERSION# (required)        : R Version 4.2.1 (2022-06-23)
  #INFRASTRUCTURE                      : R studio server
  #-----------------------------------------------------------------------------
  #     #Author                 Code History Description             Date Change
  #---- ------------       ---------------------------------------  ---------
  #1.0  Danni Yu                Last Revision                        06/30/2022
  #*****************************************************************************
}
# rm(list=ls())

#program name, track_info
if(TRUE){
  progfile<-'nzi_source.r'
}

#source functions for NZI Bayesian
if(TRUE){
  library(ggplot2)
  library(gridExtra)
  library(ggrepel)
  library(extraDistr)
  library(cowplot)
  library(dplyr)
  library(scales)
  library(ggthemes)
  
  
  #----------------------------------------------------------------------------#
  #### given percentile return number of events under beta binomail #####
  qbbinom<-function(n=10, a1=2, b1=3, pcnt=c(0.025, 0.975)){
    qq<-NULL
    Xs<-0:n
    for(pc in pcnt){
      all.p<-extraDistr::pbbinom(Xs, n, a1, b1)
      wh0<-which(all.p<=pc)
      if(length(wh0)>0){qq<-c(qq, Xs[max(wh0)])}else{qq<-c(qq, 0)}
    }
    names(qq)<-pcnt
    return(qq)
  }
  
  
  #----------------------------------------------------------------------------#
  #### Recommend the maximum P_a is at least ... #####
  recPa <- function(n, dp=6){
    if(is.null(n)){return("0.2")}
    if(is.character(n)){n<-round(as.numeric(n))}
    if(is.na(n)){return("0.2")}
    rp<-1/n
    for(i in 2:dp){
      rPa<-round(rp, i)
      if(rPa>0) break
    }
    if(i==dp & rPa==0){
      rPa<-as.numeric(paste0('0.',paste(rep(0,dp-1), collapse=''),'1'))
    }
    rPa<-min(rPa,1)
    return(rPa)
  }
  
  
  #----------------------------------------------------------------------------#
  #### compute CI (Wilson method) for p in binomial distribution ####
  binomCI <-function(y, n, alpha=.05){ #event, sample size, significance level
    phat<-y/n
    z<-qnorm(p=1-alpha/2, lower.tail=TRUE)
    tm0 <- (phat+z^2/(2*n))/(1+z^2/n)
    tm2 <- z*sqrt(phat*(1-phat)/n + z^2/(4*n^2)) / ((1+z^2/n))
    pL <- max(tm0 - tm2, 0)
    pH <- min(tm0 + tm2, 1)
    return(c(pL, pH))
  } #  binomCI(1, 20)
  
  
  #----------------------------------------------------------------------------#
  #~~~ check and plot NZI Beta distribution
  fbet<-function(bet, B=5000, p0=0.01, alp=NULL, xlim1=c(0,0.1),
                 ylim1=NULL, is.plot=F, ... ){
    if(is.null(alp)){
      alp<-2; tm.p<-min(0.05, p0); bet<-2*(0.1-tm.p)/tm.p+1;
    }
    ss<-rbeta(B, shape1=alp, shape2=bet) 
    if(is.plot){
      hist(ss, freq=F, ylim=ylim1, 
           xlim=xlim1, 
           main=paste0('p=', p0, ', beta=', bet), ...)
    }
    print(c(mean(ss), alp/(alp+bet)))
    oo<-density(ss, from=0, to=1)
    oo$y[oo$y<=0]<-min(oo$y[oo$y>0])
    return(oo)
  }
  
  #~~~ function to generate density values from the prior Beta distribution
  s.dbeta<-function(p,a,b){
    d<-dbeta(p,a,b)
    return(d)
  }
  
  #~~~ function to generate density values from the beta binomial distribution
  s.dbb<-function(x, n, a,b){
    d<- extraDistr::dbbinom(x, n, a, b)
    return(d)
  }
  
  #~~~ function to generate beta-binomial upper tail probability P[X>=x]
  #=== prior predictive
  s.pbb<-function(x, n, a,b){
    pp<- extraDistr::pbbinom(x, n, a, b, FALSE) +
      extraDistr::dbbinom(x, n, a, b)
    return(pp)
  }

  #~~~ function to generate posterior (Beta) upper tail probability
  s.pb<-function(x,n,a,b, q.p, lowt=T){
    pb<-pbeta(q.p, a+x, b+n-x, lower.tail=lowt)
    return(pb)
  }
  
  #~~~ convert string to vector
  c2v<-function(xx){
    if(is.character(xx)){xx<-eval(parse(text=paste0('c(',xx,')')))}
    return(xx)
  }
  
  #----------------------------------------------------------------------------#
  #### Re-set the structure for table output #####
  laytab<-function(ttb, ttb.th=0.99){
    colnames(ttb) <- paste0('n=', 1:ncol(ttb)-1)#NULL
    rownames(ttb) <- c("Number of Events",  #x
                       "Event Rate",        #orr
                       "Posterior Prob(Po.LT)",       #csr or po.lt
                       "FPR (PrP.UT)",      #prob.fp or prp.ut 
                       "FNR (PoP.LT)",     #prob.fn or pop.lt
                       "phat", 
                       "phat low", 
                       "phat high" )
    if(all(!is.na(ttb[3,]))){
      return( ttb[,(ttb[3,]<=ttb.th)] )
    }else{return(ttb) }
  }
  
  
  #----------------------------------------------------------------------------#
  ### count number of decimals ###
  ndecim <- function(x) {
    ndec0<-0
    if (!is.na(x) && (x %% 1) != 0) {
      tnm<-as.character(x)
      if(grepl('.', tnm, fixed=TRUE)){
        ndec0<-nchar(strsplit(tnm, ".", fixed=TRUE)[[1]][[2]])
      }else if(grepl('-', tnm, fixed=TRUE)) {
        ndec0<-as.numeric(strsplit(tnm, "-", fixed=TRUE)[[1]][[2]])
      }
    } 
    return(ndec0)
  }
  
  
  #----------------------------------------------------------------------------#
  # Bayesian approach for estimation and probabilities 
  nzi_Bayes_sim<-function(
    N=45,              #sample size
    pa=0.02,           #a maximally acceptable rate
    bet=NULL, a1=NULL, #the parameters in the prior Beta distribution
    ndec=3             #number of decimals
  ){
    #set up prior values
    if(is.null(pa) & (is.null(a1)|is.null(bet))){return(NULL)}
    if(!is.null(pa) & (is.null(a1)|is.null(bet))){ #if choose NZIB
      a1<-pa^2/(1-pa^2); 
      #a1<-pa/(1-pa); 
      bet<-1
    }else if(!is.null(a1) & is.null(pa) & !is.null(bet)){
      pa<-sqrt(a1/(a1+bet))
    }
    
    #list for all numbers of events
    if(!is.null(pa)){
      tab <- sapply(0:N, function(x, n=N){
        #--- occurrence rate
        orr <- round(x/n, ndec) 
        #--- critical success factor: posterior lower tail probabilities (beta) 
        po.lt <- round(pbeta(pa, a1+x, bet+n-x), ndec) #csf
        #--- Bayesian prior predictive upper tail probability (beta binomial)
        prp.ut<-round(s.pbb(x,n,a1,bet), ndec) #prob.fp
        #--- Bayesian posterior predictive lower tail probability 
        pop.lt<-1-round(s.pbb(x,n,a1+x,bet+n-x), ndec) # prob.fn
        #--- estimate and 95%CI for the occurrence rate
        ph <- (a1+x)/(a1+x+bet+n-x)
        ph.L <- qbeta(0.025, a1+x, bet+n-x)
        ph.H <- qbeta(0.975, a1+x, bet+n-x)
        return(c(x, orr, po.lt, prp.ut, pop.lt, ph, ph.L, ph.H))
      })
      mtt<-laytab(tab, 1.1) 
    }else{mtt<-NULL}
    
    return(laytab=mtt)
  } # nzi_Bayes_sim()
  
  
  #----------------------------------------------------------------------------#
  # Under Binomial model, get estimation and probabilities
  binom_sim<-function( 
    N=45,     #sample size
    p1=NULL,  #a lower referral rate for the null hypothesis
    pa=0.02,   #maximally acceptable rate
    getPlot=FALSE, getPnt=FALSE, getPnt.col='black', getPnt.lty=1, ndec=3,
    out='pp'  #or 'laytab'
  ){
    if(is.null(pa))return(NULL)
    
    #list for all numbers of events
    o<-0:N #number of events
    if(!is.null(p1)){
      tab <- sapply(o, function(x, n=N){
        orr <- round(x/n,ndec)
        csf <- NA
        prob.fp <- round(1-pbinom(x-1, n, pa), ndec)
        if(p1>pa){
          prob.fn <- round(pbinom(x-1, n, p1), ndec) 
        }else{prob.fn<-NA}
        ph <- x/n
        return(c(x, orr, csf, prob.fp, prob.fn, ph))
      })
      pLH <- matrix(unlist(lapply(o, binomCI, n=N, alpha=0.05)), nrow=2)
      mtt<-laytab(rbind(tab, pLH)) 
    }else{mtt<-NULL}
    
    #probability of observing at least #events: P(X>=x)
    if(!is.null(p1) && is.numeric(p1) && p1>pa){
      prob0<-pbinom(o-1, N, pa, lower.tail=F) #P(Y>=y)
      proba<-pbinom(o-1, N, p1, lower.tail=F)
      pp<-rbind(prob0, proba)
      colnames(pp)<-o
      rownames(pp)<-paste0('p=', c(pa, p1))
      if(getPlot){
        plot(prob0~o, type='l', xlim=c(0,9), ylim=c(0,0.6),
             ylab='Chance observing >= #events', 
             xlab='#events')
        points(proba~o, type='l', col=getPnt.col, lty=getPnt.lty)
        legend('topright', bty='n', legend=paste0(c('p1=','pa='), c(p1, pa)),
               lty=c(1,2), col=c('black', 'blue'), text.col=c('black', 'blue'))
      }
    }else{
      proba<-pbinom(o-1, N, pa, lower.tail=F) #P(Y>=y)
      pp<-proba
      names(pp)<-o
      if(getPlot){
        plot(proba~o, type='l',  col='blue', lty=2)
      }
      if(getPnt){
        points(proba~o, type='l',  col='blue', lty=2)
      }
    }
    #to get a normalized value relative to the center of data for model risk
    return(laytab=mtt)
  } # binom_sim(p1=0.01)
  
  
  #----------------------------------------------------------------------------#
  #visualize the prior Beta distribution and the Beta-Binomial distribution
  bbin<-function(
    prior.a="0.5,1", #parameters in p~Beta(a=, b=) prior distribution
    prior.b="0.5,1",  #parameters in p~Beta(a=, b=) prior distribution
    N11=30,           #planned sample size
    pa=0.2,           #expected average for event rate p
    p.xL='',          #xlim for beta distribution
    x.xL='',          #xlim for beta-binomial distribution
    showBinom=F,      #whether show binomial distribution in the plot
    drawPlot=TRUE,    #if F, only return the ggplot object for predictive distb
    valLeg=FALSE      #whether show the original parameter formular or value
  ){
    if(is.character(prior.a)){
      a<-c2v(prior.a); 
      pra<-strsplit(prior.a, ',')[[1]]
    }else{a<-pra<-prior.a}
    if(is.character(prior.b)){
      b<-c2v(prior.b)
      prb<-strsplit(prior.b, ',')[[1]]
    }else{b<-prb<-prior.b}
    if(valLeg){lab.a<-a; lab.b<-b}else{lab.a<-pra; lab.b<-prb}
    print(lab.a)
    if(is.character(p.xL)){xL.p<-c2v(p.xL)}
    if(is.character(x.xL)){xL.x<-c2v(x.xL)}
    if(is.character(pa)){pa<-as.numeric(pa)}
    if(is.character(N11)){N11<-round(as.numeric(N11))}
    if(length(a)!=length(b)){
      num<-min(length(a), length(b))
      a<-a[1:num]; b<-b[1:num];
    }else{num<-length(a)}
    col2<-rep(c('red', 'purple','magenta','pink', 'brown','valet'), each=4)
    col3<-c('black','blue')
    lin2<-rep(3:6,len=length(col2))
    lin3<-1:2
    
    #function to set x-axis limit
    xyL0<-function(pp, xl){
      if(is.numeric(xl) & length(xl)==2) pp<-pp+ xlim(xl)
      return(pp)
    }
    
    #generate prior distribution
    x1<-0:N11  #use exact probability instead of simulations
    u1<-x1/N11
    ppL<-bp2<-list()
    uu2<-seq(0, 1, len=50) #for beta prior density
    
    legLab1<-legLab2<-NULL #used for math expression in legend
    if(num>0){
      for(i in 1:num){
        #get density values
        bp2[[i]]<- data.frame(p=uu2, d=s.dbeta(uu2, a[i],  b[i]),
                             model=paste0('Beta(a=',lab.a[i],', b=',lab.b[i],')'))
        ppL[[i]]<- data.frame(p=u1, x=x1,
                              d=s.dbeta(u1, a[i],  b[i]),
                              dx=s.dbb(x1, N11,a[i], b[i]),
                              px=s.pbb(x1, N11,a[i], b[i]),
                              model=paste0('Beta(a=',lab.a[i],', b=',lab.b[i],')'))
        #add math expression in legend
        legLab1<-c(legLab1, 
                   bquote(paste("Beta(", alpha, "=", .(lab.a[i]),
                                ", ", beta, "=", .(lab.b[i]),')')))
        legLab2<-c(legLab2, 
                   bquote(paste("BetaB(N=", .(N11), ", ", alpha, "=", .(lab.a[i]),
                                ", ", beta, "=", .(lab.b[i]),')')))
      }
      pp <- do.call(rbind, ppL)
      bb <- do.call(rbind, bp2)
    }else{pp<-NULL}
    col1<-rep(col2, len=num) 
    lin1<-rep(lin2, len=num)
    
    #set NZIB label
    model.nzi1<- bquote(paste("Beta(", #alpha, "=", 
                              {.(pa)}^2, #{.(pa)},
                              "/(1-", {.(pa)}^2, #{.(pa)},
                              "), ", #beta, "=",
                              "1",')'))
    model.nzi2<- bquote(paste("BetaB(", .(N11), ", ", #alpha, "=", 
                              {.(pa)}^2, #{.(pa)},
                              "/(1-", {.(pa)}^2, #{.(pa)},
                              "), ", #beta, "=",
                              "1",')'))
    
    #combine 3 model data 
    tit.s<-NULL
    if(is.numeric(pa) && !is.na(pa)){
      nzi.a<-pa^2/(1-pa^2);  #nzi.a<-pa/(1-pa);  
      nzi.b<-1
      nbp2 <- data.frame(p=uu2, d=s.dbeta(uu2, nzi.a,  nzi.b),
                        model=paste0("NZIB Beta pa=", pa))
      tm0<-data.frame(p=u1, x=x1,
                      d=s.dbeta(u1, nzi.a, nzi.b),
                      dx=s.dbb(x1, N11, nzi.a, nzi.b),
                      px=s.pbb(x1, N11, nzi.a, nzi.b),
                      model=paste0("NZIB Beta pa=", pa) )
      pp <- rbind(tm0, pp)
      bb <- rbind(nbp2, bb)
      legLab1<-c(model.nzi1, legLab1)
      legLab2<-c(model.nzi2, 
                 #bquote(paste("NZIB(N=", .(N11),", ", p[a],"=", .(pa), ")")),
                 legLab2)
      col1<-c(col3[1], col1)
      lin1<-c(lin3[1], lin1)
      if(showBinom){
        tm<-data.frame(p=rep(pa,length(x1)), x=x1, 
                       d=rep(1, length(x1)),
                       dx=dbinom(x1, N11, pa),
                       px=1-pbinom(x1-1, N11, pa),
                       model=paste0('Binomial(N=',N11,', p=',pa, ')'))
        pp2<-rbind(pp, tm)
        legLab2<-c(legLab2, 
                   bquote(paste('Binomial(N=',.(N11),', p=', .(pa), ')')))
        col1<-c(col1, col3[2])
        lin1<-c(lin1, lin3[2])
        tit.s<-'or P(Y=y) for Binomial'
      }else{pp2<-pp}
    }else{pp2<-pp}
    
    #generate beta-binomial distribution
    pp2$model[grepl('NZIB', pp2$model, fixed=TRUE)]<-
      paste0("NZIB(N=", N11, ', pa=', pa, ')')
    pp2$model<-gsub('Beta(', paste0('BetaB(N=',N11, ', '), pp2$model, fixed=T)
    pp2$model<-factor(pp2$model, levels=unique(pp2$model))
    #for beta prior distribution
    pp$model<-factor(pp$model, levels=unique(pp$model))
    bb$model<-factor(bb$model, levels=unique(bb$model))
    
    #clean data for predictive distribution
    pp2<-unique(pp2[,c('x','dx', 'px','model')])
    
    #generate distribution plot for prior Beta
    p1<-ggplot(bb, aes(x=p, y=d, group=model,color=model, linetype=model)) + 
      geom_line() + ylab('Density') + 
      theme_bw() + ggtitle('Prior distribution') +
      xlab('p') + 
      scale_color_manual(values=col1, labels=legLab1)+
      scale_linetype_manual(values=lin1, labels=legLab1)
    p1<-xyL0(p1, xL.p)
    
    if(drawPlot){
      #generate distribution plot for Beta-Binomial
      p2<-ggplot(pp2, 
                 aes(x=x, y=dx, group=model, color=model, linetype=model)) + 
        geom_line() + ylab('Density') + theme_bw() + 
        ggtitle('Prior Predictive Distribution', subtitle=tit.s) +
        xlab('y: number of events') + 
        scale_color_manual(values=col1, labels=legLab2)+
        scale_linetype_manual(values=lin1, labels=legLab2)
      p2<-xyL0(p2, xL.x)
      
      p3<-ggplot(pp2, 
                 aes(x=x, y=px, group=model, color=model, linetype=model)) + 
        geom_line() + ylab('Probability') + theme_bw() + 
        ggtitle('Prior Predictive Upper Tail') +
        xlab('n: number of events') + 
        scale_color_manual(values=col1, labels=legLab2)+
        scale_linetype_manual(values=lin1, labels=legLab2)
      p3<-xyL0(p3, xL.x)
      
      grid.arrange(grobs=list(p1, p2), nrow=1, widths=c(1, 1.1))
    }else{
      #generate distribution plot for Beta-Binomial
      p2<-ggplot(pp2, aes(x=x, y=dx, group=model,color=model, linetype=model)) + 
        geom_line() + 
        theme_bw() + ggtitle('Prior Predictive distribution', subtitle=tit.s) +
        scale_color_manual(values=col1, labels=legLab2)+
        scale_linetype_manual(values=lin1, labels=legLab2)
      p2<-xyL0(p2, xL.x) 
      return(list(data=pp2, plot=p2))
    }
  } # bbin()

  
  #----------------------------------------------------------------------------#
  #source functions for operating characteristics and calculate probstop
  if(TRUE){
    ######Generating Toxicity Monitory Rule#####
    #compare to online shiny App with minor revision for small p
    #https://www.trialdesign.org/one-page-shell.html#BTOX
    #rule based on Prob(theta>th0)<stPr
    BayTox <- function(
      betaPrior="1, 1", #Prior distribution "alpha, beta"
      N=10,        #maximum number of subjects
      n0=3,        #minimum number of subjects if early stop by enough ADRs
      th0=0.1,     #maximum ADR rate
      stPr=0.05,   #stopping criterion probability for excessive toxicity
      coh.size=1,  #number of patients in each cohort run
      ruleTp=1,    #stopping rule type: 1=posterior>0.95; 2=priorPred<0.5; 
      #3=binomial p-value
      adj.pv='BH'  #the p-value adjustment method
    ){
      pa<-c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr")
      #if adj.pv is NULL, then no adjustment
      
      #check input
      if(!is.null(n0) && !is.null(N) && n0>N-1){n0<-N}
      if(!is.null(n0) && !is.null(N) && n0<1){n0<-1}
      if(is.null(th0)|is.na(th0) | is.null(stPr)|is.na(stPr)){return(NULL)}
      if(n0<2){stop('minimum number of subject n0 must be >1!')}
      
      #If Beta Prior is NULL or '', set it for NZIB
      if(is.null(betaPrior)||gsub(' ', '', betaPrior)==''){
        ab<-c(th0^2/(1-th0^2), 1) ;# ab<-c(th0/(1-th0), 1)
        pltp<-1
      }else{ #get hyper prior parameters
        ab<-as.numeric(unlist(strsplit(betaPrior, split=',', fixed=T)))
        if(length(ab)<2){return(NULL)}
        pltp<-3
      }
      bpr<-paste0("c(alpha=", ab[1], ", beta=", ab[2], ")")
      
      #set label for 1st row in the stopping rule tables
      lb<-max(n0-1, 1, coh.size-1)
      row01<-ifelse(lb>1, paste0('1-', lb), '1')
      row02<-'Not Applicable'
      
      #set a contain for all probabilities
      p0_ij<-matrix(NA, nrow=N, ncol=N)
      
      #set number of patients, number of ADRs
      max.j<- min(n0, coh.size) - 1
      i <- 1      #---i: number of ADR events; j: number of patients---#
      #generate all the probabilities
      while(max.j<=N & i<=N){
        j <- ceiling(i/coh.size)*coh.size
        while(j>=i & j<=N){#fix i, increase j
          if(ruleTp==3){#binomial p-value
            p0_ij[i,j]<- 1-pbinom(i-1, j, th0)
          } else if(ruleTp==2){#prior predictive p-value
            p0_ij[i,j]<- s.pbb(i, j, ab[1], ab[2])
          }else{#posterior probability
            p0_ij[i,j]<- pbeta(th0, ab[1]+i, ab[2]+j-i) 
          }
          j <- j+coh.size
        }
        i <- i + 1
      }
      
      #p-value adjust method
      if(!is.null(adj.pv) && adj.pv%in%pa){
        p_ij<-matrix(p.adjust(p0_ij, method=adj.pv), nrow=N)
      }else{p_ij<-p0_ij}
      
      #set number of patients, number of ADRs
      ns <- ets <- NULL 
      i <- 1      #---i: number of ADR events; j: number of patients---#
      prior.j<-n0
      while(max.j<=N & i<=N){
        j <- ceiling(i/coh.size)*coh.size
        tm.j<- NULL
        #search largest j(=#patients) for smallest p_ij (>stPr)
        while(j>=i  && j<=N && p_ij[i,j]<stPr){
          if(j>max.j){tm.j <- c(tm.j, j)}
          j <- j+coh.size
        }
        tm.j <- tm.j[!is.na(tm.j)&!is.infinite(tm.j) & is.numeric(tm.j)]
        if(length(tm.j)>0){
          max.j <- max(tm.j)
          ns0 <- unique(c(min(min(tm.j),N), min(max.j, N)))
          ns0 <- ns0[ns0>=n0]
          if(length(ns0)>0){
            ets <- c(ets, i)
            if(length(ns0)==1){ns0<-c(prior.j, ns0); prior.j<-NULL}
            ns <- c(ns, paste(unique(ns0), collapse='-'))
          }
        }
        i <- i + 1
      }#get: ns and ets
      
      if(n0>0){
        oT <- data.frame(x=c(row01,ns), y=c(row02, ets))
      }else{
        oT <- data.frame(x=ns, y=ets)
      }
      xx <- unique(oT$x)
      oT.od <- data.frame(od=1:length(xx), x=xx)
      oT <- aggregate(y~x, data=oT, FUN=min)
      oT <- merge(oT.od, oT, by='x', sort=F)[,c('x','y')]
      colnames(oT)<-c('Total Number of Treated Subjects', 
                      'Suspending Accrual if >= Number of Subjects with the ADR')
      print(oT)
      return(oT)
    }
    
    ######Run Operating Characteristics based on BayTox result#####
    BayToxOp<-function(
      #true ADR rates:
      trueProb='0.0000001, 0.001, 0.01,0.02,0.05, 0.1', 
      #get the threshold: 
      th=BayTox(), 
      #percentiles for number of patients:
      spct='0.10, 0.25, 0.5, 0.75, 0.90',
      only.early=FALSE, #if TRUE, return ProbStop as stop prob before reaching N
      est.n=FALSE #if TURE, return estimated n values at spct percentiles
    ){
      th2x<-th2y<-NULL
      for(i in 1:nrow(th)){
        n.rg <- unlist(strsplit(th[i,1], split='-'))
        suppressWarnings( n.rg<-unique(as.numeric(n.rg)) )
        suppressWarnings( evt <-as.numeric(th[i,2]) )
        if(is.na(evt)|any(is.na(n.rg))){next}
        if(length(n.rg)>1){
          nn <- min(n.rg): max(n.rg)
        }else if(length(n.rg)==1){
          nn <- n.rg[1]
        }else{nn<-NULL}
        if(!is.null(nn)){
          th2x<-c(th2x, nn) #sample sizes
          th2y<-c(th2y, rep(th[i,2], length(nn))) #number of events
        }   
      }
      #a dataframe with sample sizes and number of events
      th2 <- unique(data.frame(x=as.numeric(th2x), y=as.numeric(th2y)))
      th2 <- th2[!is.na(th2$x)&!is.na(th2$y),]
      th2 <- th2[order(th2$x),]
      
      #sample size percentiles included in results for estimation on n
      n.pct <- unique(unlist(strsplit(spct, split=',', fixed=T)))
      n.header <- gsub('0.', 'n.pct.', as.character(n.pct), fixed=T)
      n.pct <-as.numeric(n.pct)
      n.sel <- !is.na(n.pct)
      n.pct <- n.pct[n.sel]
      n.header <- n.header[n.sel]
      
      #convert true rates into a vector
      tps <- as.numeric(unlist(strsplit(trueProb, split=',', fixed=T)))
      tps <- unique(tps[!is.na(tps)])
      
      #setup operating table
      oT <- data.frame(x=1:length(tps), y=tps)
      if(nrow(th2)==0){
        oT <- cbind(oT, z=0)
        oT[,2]<-format(oT[,2], scientific=F)
        colnames(oT)<-c('Scenario', 'TrueRate', 'ProbStop')
        return(oT)
      }
      z <- ths <- NULL
      nSS <- matrix(NA, nrow=length(tps), ncol=length(n.header)+2)
      nSS.id<-1
      for(tp in tps){
        #prob of stopping at 1st check
        p0 <- pbinom(th2$y[1]-1, th2$x[1], tp, F) #Prob(stop.at.th2$x[1])
        evG <- 0:(th2$y[1]-1) #all numbers of events if not stop at 1st check
        #prob for each number of events at 1st check if not stop
        pG <- dbinom(evG, th2$x[1], tp)
        #prob of not stop at 1st check
        p2 <- sum(pG) #Prob(n>th2$x[1])
        if(nrow(th2)>1){ #start ith check after 1st check
          for(i in 2:nrow(th2)){ 
            evS2 <- th2$y[i]-1 - evG #all numbers of new events at ith check
            #prob for each number of new events at ith check
            pS <- pbinom(evS2, th2$x[i]-th2$x[i-1], tp, F) 
            #update prob of stopping at ith check
            p0 <- c(p0, sum(pG*pS)) #Prob(n=th2$x[i]) stop at th2$x[i]
            if(i<nrow(th2)){
              evG2 <- 0:(th2$y[i]-1) #all numbers if not stop at ith check
              pG2 <- NULL
              #prob for each number of events if not stop at ith~1st checks
              for(k in evG2){
                pG2<-c(pG2, sum(dbinom(k-evG, th2$x[i]-th2$x[i-1], tp)*pG))
              }
              evG <- evG2 #will be used for next iteration
              pG <- pG2
            }
            #update prob of not stop at ith check
            p2 <- c(p2, sum(pG)) #Prob(n>th2$x[i])
          }
        }
        names(p0)<-th2$x #probabilities stopping at the ith checks.
        if(only.early){
          z <- c(z, sum(p0[-length(p0)])) #prob.early.stop before reaching n 
        }else{
          z <- c(z, sum(p0)) #prob.stop before reaching or at n
        }
        #get sample size percentiles
        th2$p.stop<-p0
        th2$p.go<-p2
        th2$n.pct<-1-th2$p.go
        thS <- rbind(ths, th2)
        for(k in 1:length(n.pct)){
          dst <- th2$n.pct-n.pct[k]
          dst2 <- abs(dst)
          #wh2 <- which(dst2<=n.pct.dist)#nearest method
          wh0 <- which(th2$n.pct < n.pct[k]) #step up method
          if(length(wh0)>0){wh0 <-min(max(wh0) +1, nrow(th2))}else{wh0<-1}
          wh3 <- wh0#[wh0%in%wh2]
          if(length(wh3)<1){wh3<- max(which(dst2==min(dst2)))}
          nSS[nSS.id, k] <- th2$x[wh3]
        }
        if(nrow(th2)>1){
          th2$p.stop[nrow(th2)]<-th2$p.go[nrow(th2)-1]
        }
        n.mu<-sum(th2$x*th2$p.stop)
        n.var<- sum( ((th2$x-n.mu)^2)*th2$p.stop )
        nSS[nSS.id, k+1] <-round(n.mu,4) #get the average mean
        nSS[nSS.id, k+2] <-round(sqrt(n.var),4) #get the SD
        nSS.id <- nSS.id+1
      }
      
      oT$z <- round(z,4) #Prob.Early.Stop
      oT <- cbind(oT, nSS)
      oT[,2]<-format(oT[,2], scientific=F)
      colnames(oT)<-c('Scenario', 'TrueRate', 'ProbStop', n.header, 
                      'n.Mean', 'n.SD')
      if(!est.n){oT<-oT[,1:3]}
      return(oT)
    }
    
    ##################################################
    #Check whether a rule is successfully created
    IsRuleSet<-function(md=mods[i], hh=hh1, th.hh=th2){
      if(nrow(hh)==1 & hh[1,2]=="Not Applicable"){
        print(su1<-paste(md, 'fails to create stopping rule for', th.hh))
        return(list(is.fail=TRUE, subtitle=su1))
      }else{return(list(is.fail=FALSE))}
    }
    
    ##################################################
    #generate 1 set of operating characteristics
    # OC1 <- function(ab='0.5,0.5', NN1=NN, th1=th2, tp1=tp ) {
    #   modc<-rep(c(paste0(ab) ), each=2)
    #   mods<-paste0('Bayesian ', NN1, ",", modc)
    #   mod0<-rep(paste0('NZIB(',NN1, ', ', th1,')'), each=2)
    #   modLv<-paste(c(mods,mod0), c('by Po.LT', 'by PrP.UT'))
    #   nLv<-length(modLv)
    #   sz <- (1:nLv)*1.1/nLv
    #   subtt<-NULL #used for subtitle
    #   
    #   btL<-list(); id<-1;
    #   for(i in 1:length(modc)){
    #     if((i%%2)==1){rtp<-1}else{rtp<-2}
    #     hh1<- BayTox(betaPrior=modc[i], N=NN1, th0=th1, ruleTp=rtp)
    #     irs1<-IsRuleSet(md=modLv[i], hh=hh1, th.hh=th1)
    #     if(irs1$is.fail){subtt<-c(subtt, irs1$subtitle)}else{
    #       btL[[i]]<-data.frame(model=modLv[i], BayToxOp(th=hh1,trueProb=tp1) )
    #     }
    #   }
    #   for(j in 1:2){
    #     hh01<-BayTox(betaPrior="", N=NN1, th0=th1, ruleTp=j) 
    #     irs01<-IsRuleSet(md=modLv[i+j], hh=hh01, th.hh=th2)
    #     if(irs01$is.fail){subtt<-c(subtt, irs01$subtitle)}else{
    #       btL[[i+j]]<-data.frame(model=modLv[i+j], BayToxOp(th=hh01,trueProb=tp1) )
    #     }
    #   }
    #   
    #   #integrate datasets
    #   intd<-rbind( bind_rows(btL))
    #   intd$N<-NN1; intd$pa<-th1;
    #   intd$model<-factor(intd$model, levels=modLv)
    #   intd$TrueRate<-as.numeric(intd$TrueRate)
    #   
    #   return(intd)
    # } 
    
  }#end probstop functions
  
  
  #----------------------------------------------------------------------------#
  #---generate probability coverage plot---#
  prob_cov<-function(
    sample_size = c(45),  #number of samples in one trial
    show_freq = 1,        #observed event frequency in a trial
    p_a = c("0.0001"), #max acceptable occurrence rate p
    ab = c(0.5, 0.5),  #the alpha and beta in the other Bayesian model
    a.seq=seq(0.05, 0.95, by=0.01), #sequence of statistic significance levels
    lenLab=5, #number of breaks at Y-axis
    lab=c('NZIB', 'Bayesian', 'Binomial'),
    tt0='Posterior Probability Coverage' #title
  ){
    l.bay<-grepl('baye', tolower(lab))
    l.bin<-tolower(lab)=='binomial'
    if(any(l.bay)){lab[l.bay]<-paste('Bayesian', paste(ab, collapse=','))}
    
    if(is.character(p_a)){p_a<-as.numeric(unlist(strsplit(p_a, split=',')))}
    
    #set vectors
    bb<-nn<-aa<-yy <-pp <-ci<-pL<-pH<- NULL
    for(i in sample_size){
      for(k in show_freq){
        for(m in a.seq){
          for(b0 in lab){
            if(b0=='NZIB'){
              for(j in p_a){
                lb0<-paste(b0, j)
                bb<- c(bb, rep(lb0, 2))
                nn<- c(nn, rep(i, 2))
                yy<- c(yy, rep(k, 2))
                a2t<-c(m/2, 1-m/2)
                aa<- c(aa, a2t[1], a2t[1])
                j<-as.numeric(j)
                #get the posterior prob cov for NZIB
                ptt<-qbeta(a2t, j^2/(1-j^2)+k, 1+i-k)
                # ptt<-qbeta(a2t, j/(1-j)+k, 1+i-k)
                pp<-c(pp, ptt)
                ci<-c(ci, rep((1-m)*100, 2))
                pL<-c(pL, rep(ptt[1], 2))
                pH<-c(pH, rep(ptt[2], 2))
              }
            }else{
              bb<- c(bb, rep(b0, 2))
              nn<- c(nn, rep(i, 2))
              yy<- c(yy, rep(k, 2))
              a2t<-c(m/2, 1-m/2)
              aa<- c(aa, a2t[1], a2t[1])
              j<-as.numeric(j)
              if(b0=='Binomial'){ ptt<-binomCI(k, i, m)
              }else{ptt<-qbeta(a2t, ab[1]+k, ab[2]+i-k) }
              #get the prob cov for Bay or Binomial
              pp<-c(pp, ptt)
              ci<-c(ci, rep((1-m)*100, 2))
              pL<-c(pL, rep(ptt[1], 2))
              pH<-c(pH, rep(ptt[2], 2))
            }
          }
        }
      }
    }
    lev<-c(paste(lab[1], p_a), lab[-1])
    pcd<-unique(data.frame(lab=bb, n=nn, alpha=aa, y=yy, phat=pp, ci=ci))
    pcd$lab<-factor(pcd$lab, levels=lev)
    pcd <- pcd[order(pcd$y, pcd$phat),]
    pcd <- pcd[order(pcd$lab, pcd$n),]
    pcd2<-unique(data.frame(lab=bb, n=nn, alpha=aa, y=yy, ci=ci, pL=pL, pH=pH))
    if(TRUE){
      col1 <- c(rep("black", length(p_a)), 
                ifelse(any(l.bay), 'red', NA),
                ifelse(any(l.bin), 'blue', NA))
      col1 <- col1[!is.na(col1)]
      #get Y-labels
      a0<-sort(unique(pcd$alpha))
      len.a<-length(a0)
      a.rg<-a0[seq(1, len.a, len=min(lenLab, len.a))]
      u.a <- unique(pcd[,c('alpha', 'ci')])
      u.a <- u.a[u.a$alpha%in%a.rg, ]
      a.lb<-paste0(round(u.a$ci), '% CI')
      names(a.lb)<-as.character(u.a$alpha)
      #generate the data to show the plot
      pcd$Y<-paste0('y=', pcd$y)
      pcd$n<-paste0('N=', pcd$n)
      p0<-ggplot(data=pcd, 
                 aes(x=phat, y=alpha, color=lab, linetype=lab, fill=lab)) + 
        geom_ribbon(aes(x=phat, ymax=alpha, fill=lab), 
                    ymin=min(a0), alpha=0.2)+
        geom_line() + theme_bw() + 
        theme(plot.title=element_text(size = 11), 
              legend.position="bottom", legend.key.width=unit(0.008, 'npc') )+
        ggtitle(tt0) + ylab('') + 
        xlab('Estimated Rate') + scale_linetype('') +  
        scale_colour_manual('', values=col1) +
        scale_fill_manual('', values=col1) +
        scale_y_continuous(breaks=u.a$alpha, labels=a.lb)
      if(length(sample_size)>1){
        p0 <- p0 + facet_grid(n~Y, switch='y')
      }else{
        p0 <- p0 + facet_grid(.~Y, switch='y')
      }
    }
    return(list(pcd, pcd2, plot=p0))
  }  # prob_cov()
  
  
  #----------------------------------------------------------------------------#
  #All-in-one function for NIZB analysis
  nzib <- function( stopAt=NULL, #number of ADR for stopping rule
                    N1=45,         #max number of samples
                    maxIDR=0.0051, #max expected acceptable occurrence rate
                    normIDR=NULL,  #normal occurrence rate, default maxIDR/2 
                    plot.type=4,   #1. NZIB; 2. NZIB, Binomial; 
                    #3.NZIB, Bayesian; 4.all 3 lines
                    bay.prior="c(alpha=0.05, beta=0.05)", #hyper parameter in beta-binomial
                    col1=c("black",'red', 'blue'),     #line colors in graph
                    rsk.xlim=NULL,    #x-axis range for risk plot
                    rsk.ylim=NULL,    #y-axis range for risk plot
                    errI.xlim=c(0,5), #x-axis range for type I error plot
                    errI.ylim=NULL,   #y-axis range for type I error plot
                    errII.xlim=NULL,  #x-axis range for type II error plot
                    errII.ylim=NULL,  #y-axis range for type II error plot
                    sup.xlim=NULL,    #x-axis range for superiority plot, default=errI.xlim
                    sup.ylim=NULL,    #y-axis range for superiority plot
                    plot4=FALSE,      #show 4 plots
                    goTh="0.05, 0.2", #go/no-go thresholds for type I and type II errors
                    noplot=F,         #don't create any plot if T
                    showNZIB=TRUE,    #if plot.type!=1, show NZIB or not
                    ajpm = NULL       #method for p-adjust 
  ){
    if(is.null(maxIDR)){maxIDR<-""}
    pdig<-ndecim(as.numeric(maxIDR))
    
    #set prior according to plot type
    if(plot.type%in%c(1,2)){bay.prior<-NULL}
    if(plot.type%in%c(2,4)){showbm<-T}else{showbm<-F}
    
    #clean up parameters
    if(is.null(goTh)){
      goTh<-c(0.05, 0.2)
    }else{
      goTh<-as.numeric(strsplit(as.character(goTh), split=',', fixed=T)[[1]])
    }
    if(length(goTh[!is.na(goTh)])==0){goTh<-c(0.05, 0.2)}
    if(length(goTh[!is.na(goTh)])==1){goTh<-c(goTh, 0.2)}
    if(is.character(N1)){N1<-max(3, round(as.numeric(N1)))}
    if(!is.null(stopAt) & !is.numeric(stopAt)){
      stopAt<-gsub(' ', '', stopAt, fixed=T)
      if(stopAt!=''){
        stopAt<-unlist(strsplit(as.character(stopAt), split=',', fixed=T))
        stopAt<-as.numeric(stopAt)
        stopAt<-stopAt[!is.na(stopAt)]
        if(length(stopAt)>0){
          stopAt<-pmin(pmax(0, round(stopAt)),N1) #lowest bound
        }
      }
    }
    if(is.character(maxIDR)){maxIDR<-max(min(as.numeric(maxIDR),1),0)}
    if(is.character(normIDR)){normIDR<-min(max(as.numeric(normIDR),maxIDR),1)}
    
    bay.prior<-c2v(bay.prior)
    rsk.xlim<-c2v(rsk.xlim)
    rsk.ylim<-c2v(rsk.ylim)
    errI.xlim<-c2v(errI.xlim)
    errI.ylim<-c2v(errI.ylim)
    errII.xlim<-c2v(errII.xlim)
    errII.ylim<-c2v(errII.ylim)
    sup.xlim<-c2v(sup.xlim)
    sup.ylim<-c2v(sup.ylim)
    
    #get stop reference relying number of ADR events
    if(!is.null(stopAt)&is.numeric(stopAt)&length(stopAt)>0){
      sa <- round(stopAt) #require an integer
      saLab<-paste(sort(sa), collapse=',')
    }else{sa<-NULL; saLab<-NULL}
    
    #set default normIDR
    if(is.null(normIDR)){normIDR <- maxIDR}
    
    #get normal expected rate 
    if(plot.type%in%c(1,3)){ tP1<-NULL }else{ tP1<-maxIDR }
    
    #order Bayesian hyper parameters
    bay.prior.nm<-c('alpha','beta')
    if(all(bay.prior.nm%in%names(bay.prior))){
      bay.prior<-bay.prior[bay.prior.nm]
    }
    
    #label the methods
    mNm<-c('NZIB', 'Bayesian', 'Binomial')
    
    #colors
    if(plot.type==2 & length(col1)>2){
      if(!showNZIB){
        lab2<-col1[c(3)]; lin2<-c(3)
      }else{
        lab2<-col1[c(1,3)]; lin2<-c(1,3)
      }
    }else{ lab2<-col1; lin2<-1:3 }
    
    #----get get estimates and probabilities----#
    b.eg<-eb.eg<-nzi.eg<-NULL
    #: classic binomial
    if(plot.type%in%c(2,4)){
      b.eg<-binom_sim(N=N1, p1=normIDR, pa=maxIDR)
    }
    #: Bayesian
    if(plot.type%in%c(3,4)){
      eb.eg<-nzi_Bayes_sim(N=N1,a1=bay.prior[1], bet=bay.prior[2], pa=maxIDR)
    }
    #: NZI Bayesian
    if(plot.type==1 | (showNZIB & plot.type!=1)){
      nzi.eg<-nzi_Bayes_sim(N=N1, pa=maxIDR)
    }
    
    #adjust by plot type
    if(plot.type==4){#show all 3 lines
      eg<-list(nzi.eg, eb.eg, b.eg)
      names(eg)<-mNm
    }else if(plot.type==3){#show lines for NZIB, Bayesian
      eg<-list(nzi.eg, eb.eg)
      names(eg)<-mNm[1:2]
    }else if(plot.type==2){#show lines for NZIB, binomial
      eg<-list(nzi.eg, b.eg)
      names(eg)<-mNm[c(1,3)]
    }else{#show line for NZIB only
      eg<-list(nzi.eg)
      names(eg)<-mNm[1]
    }
    if(!showNZIB & plot.type!=1){ eg<-eg[names(eg)!=mNm[1]] }
    mNM0<-names(eg)
    
    #compare prior predictive and post post predictive
    set.errs<-function(egL=eg){
      numEvts<-fpr<-fnr<-sup<-grp<-phat<-phat.L<-phat.H <-NULL
      for(i in 1:length(egL)){
        numEvts<-c(numEvts, egL[[i]][1,]) 
        fpr<-c(fpr, egL[[i]][4,]) 
        fnr<-c(fnr, egL[[i]][5,])
        sup<-c(sup, egL[[i]][3,])
        grp<-c(grp, rep(names(egL)[i], ncol(egL[[i]])))
        phat<-c(phat, egL[[i]][6,])
        phat.L<-c(phat.L, egL[[i]][7,])
        phat.H<-c(phat.H, egL[[i]][8,])
      }
      ors<- data.frame( numEvts, fpr, fnr, sup, grp, phat, phat.L, phat.H)
      ors$sen<-1-ors$fpr;  
      ors$pwr<-1-ors$fnr
      ors$grp <- factor(ors$grp, levels=names(egL))
      return(ors)
    }
    errs<-set.errs(eg)
    errs.2<-errs[order(errs$grp, errs$numEvts),]
    errs.3<-errs[!is.na(errs$sup),]
    
    #internal functions for adding stopping thresholds into plot
    xyL<-function(pp, xl, yl, ssa, myLabel, leg='yes', 
                  addref=TRUE, anno.pos=0
    ){
      if(is.numeric(xl) & length(xl)==2) pp<-pp+ xlim(xl)
      if(is.numeric(yl) & length(yl)==2) pp<-pp+ ylim(yl)
      if(!is.null(ssa)){
        if(addref){
          pp<-pp+ geom_vline(xintercept=ssa,colour='darkgreen', linetype=2)
        }
        pp<-pp+geom_text_repel( data=myLabel[myLabel$numEvts%in%ssa,],
                                aes(label = mylab), 
                                position=position_dodge(anno.pos),
                                box.padding = unit(0.35, "lines"),
                                point.padding = unit(0.3, "lines"),
                                show.legend = FALSE)
      }
      if(leg=='no'){pp<-pp+theme(legend.position='none')}else{
        pp<-pp+guides(linetype="none")
      }
      return(pp)
    }
    
    #get posterior lower tail-area probability estimation
    #Bayesian
    eb.pst <- eb.eg[grepl("Po.LT", row.names(eb.eg)),]
    #NZIB
    nzi.pst <- nzi.eg[grepl("Po.LT", row.names(nzi.eg)),]
    #integrate posterior results: nzib, bayesian
    getid<-function(n1){if(n1<=0){return(NULL)}else{1:n1}}
    psto<-cbind(nzi.pst, eb.pst)
    psto<-data.frame(observe_number=(1:nrow(psto))-1, psto)
    colnames(psto)[-1]<-paste('posterior prob', mNM0[mNM0!='Binomial'])
    pst<-data.frame(
      numEvts=c(getid(length(eb.pst)), getid(length(nzi.pst))) - 1,
      post=c(eb.pst, nzi.pst),
      grp=c(rep(mNm[2], length(eb.pst)), rep(mNm[1], length(nzi.pst)) ) )
    pst$grp <- factor(pst$grp, levels=mNM0[mNM0!='Binomial'])
    pst$mylab<-as.character(round(pst$post, 3))
    
    #setup label? double check the use case locations
    if(is.na(normIDR)|plot.type%in%c(2,4)){
      p0lab.s<-bquote(paste("P(Y">="y"," | p=", .(normIDR), ") for Binomial" ))
      ppd.y<-bquote(paste("P(Y=y)"))
    }else{
      p0lab.s<- NULL
      #ppd.y<-bquote(paste("P(",tilde(Y),"=y)"))
      ppd.y<-bquote('P('~Y^rep~'=y)')
    }
    
    #1. posterior predictive upper tail probability
    errs.2$mylab<-as.character(round(errs.2$pwr, 3))
    p1<-ggplot(data=errs.2, aes(x=numEvts, y=pwr, color=grp, linetype=grp)) + 
      geom_line() + theme_bw() + theme(plot.title = element_text(size = 11))+
      ggtitle('Posterior Preditive Upper Tail (PoP.UT)') + 
      # ylab(bquote(paste("P(",tilde(Y),"">="y)"))) + 
      ylab(bquote(paste("P(",Y^rep,"">="y|observing y)"))) + 
      xlab('Data: y') + scale_linetype('') + scale_colour_manual('', values=lab2)
    pp.poput<-xyL(pp=p1, xl=errII.xlim, yl=c(0,1), ssa=sa, myLabel=errs.2, 
                  leg='no') 
    
    #2 posterior prob plot
    pp.po <- ggplot(data=pst, aes(x=numEvts, y=post, color=grp, linetype=grp))+
      geom_hline(yintercept=goTh[1], color='tan3') + 
      geom_line() + theme_bw() + theme(plot.title = element_text(size = 11)) +
      ggtitle( bquote(paste("Posterior Lower Tail (Po.LT)")))+ 
      xlab("Data: y") + ylab(bquote(paste("P(p"<="", .(maxIDR), ~"|"~"observing y)"))) +
      scale_linetype('') + scale_colour_manual('', values=lab2) 
    pp.po<-xyL(pp=pp.po, xl=errI.xlim, yl=c(0,1), ssa=sa, 
               myLabel=pst, leg='no')
    
    #3. prior predictive probability
    p0lab<-"Priopr Predictive Upper Tail (PrP.UT)"
    errs.2$mylab<-as.character(round(errs.2$fpr, 3))
    p2<-ggplot(data=errs.2, aes(x=numEvts, y=fpr, color=grp, linetype=grp)) + 
      geom_hline(yintercept=goTh[1], color='tan3') +
      geom_line() + theme_bw() + theme(plot.title = element_text(size = 11))+
      ggtitle(p0lab, subtitle=p0lab.s) + xlab("Data: y") +
      scale_colour_manual('', values=lab2) +
      ylab(bquote(paste("P(Y">="y)")))  
    pp.prput<-xyL(pp=p2, xl=errI.xlim, yl=c(0,1), ssa=sa, myLabel=errs.2,
                  leg='yes')

    #graph
    if(plot4){
      #-------------------set up operating characteristics---------------------# 
      truep<-paste(seq(0.000001, min(1, maxIDR*2.5), len=10), collapse=',')
      mod0<-rep(paste0('NZIB(',N1, ', ', maxIDR,')'), each=2)
      subtt<-NULL #used for subtitle
      btL<-stop_rule<-list(); id<-1;
      if('Bayesian'%in%mNM0){
        modc<-rep( paste(bay.prior, collapse=','), each=2)
        mods<-paste0('BetaB(', N1, ",", modc, ")",  c(' by Po.LT', ' by PrP.UT'))
        for(i.m in 1:length(modc)){ #for Bayesian
          if((i.m%%2)==1){rtp<-1}else{rtp<-2}
          print(goTh) ######
          hh1<- BayTox(betaPrior=modc[i.m], N=N1, th0=maxIDR, ruleTp=rtp, 
                       stPr=goTh[1], adj.pv=ajpm)
          stop_rule<-c(stop_rule, list(data.frame(model=mods[i.m], hh1)))
          irs1<-IsRuleSet(md=mods[i.m], hh=hh1, th.hh=maxIDR)
          if(irs1$is.fail){subtt<-c(subtt, irs1$subtitle)}else{
            btL[[i.m]]<-data.frame(model=mods[i.m], 
                                   BayToxOp(th=hh1,trueProb=truep) )
          }
        }
      }else{mods<-modc<-NULL; i.m <- 0}
      modLv<- c(mods, paste(mod0, c('by Po.LT', 'by PrP.UT')))
      nLv<-length(modLv)
      sz <- (1:nLv)*1.1/nLv
      for(j in 1:2){ #for NZIB
        hh01<-BayTox(betaPrior="", N=N1, th0=maxIDR, ruleTp=j, 
                     stPr=goTh[1], adj.pv=ajpm) 
        stop_rule<-c(stop_rule, list(data.frame(model=modLv[i.m+j], hh01)))
        irs01<-IsRuleSet(md=modLv[i.m+j], hh=hh01, th.hh=maxIDR)
        if(irs01$is.fail){subtt<-c(subtt, irs01$subtitle)}else{
          btL[[i.m+j]]<-data.frame(model=modLv[i.m+j], 
                                   BayToxOp(th=hh01, trueProb=truep) )
        }
      }
      #integrate data set
      atd<-rbind( bind_rows(btL))
      atd$N<-N1; atd$pa<-maxIDR;
      atd$TrueRate<-as.numeric(atd$TrueRate)
      nzlab<-c('NZIB by Po.LT', 'NZIB by PrP.UT')
      atd$model[grepl('NZIB', atd$model)& grepl('Po.LT', atd$model)]<-nzlab[1]
      atd$model[grepl('NZIB', atd$model)&grepl('PrP.UT', atd$model)]<-nzlab[2]
      atd$TrueRate<-as.numeric(atd$TrueRate)
      atd$N<-paste0('N=',atd$N)
      atd$pa<-paste0('pa=',atd$pa)
      u.mod<-unique(atd$model) 
      n.nzib<-sum(grepl('NZIB', u.mod))
      nLv<-length(unique(atd$model))
      a.0<-nrow(unique(atd[,c('N','pa')]))
      a.col<- rep(c(rep('red', nLv-n.nzib), rep('black', n.nzib)), a.0)
      legt<-''
      oc_legp <- c(min(atd$TrueRate) + diff(range(atd$TrueRate))*.1, 1 )
      #####~~~~~graph~~~~####
      pp.oc<-ggplot(atd, size=5, alpha=0.5,
                    aes(x=TrueRate, y=ProbStop, color=model, shape=model,
                        linetype=model, group=model)) + 
        ggtitle('Operating Characteristics') + xlab("True Rate: p") +
        ylab(paste0('P[stop before or at N=', N1, ']')) + 
        geom_point()+ geom_line() + theme_bw() + ylim(c(0,1))+
        theme(legend.justification=c(0,0.9),
              legend.background=element_rect(fill=alpha("yellow", 0)),
              axis.text.x=element_text(angle=90, vjust=0.5, hjust=1),
              panel.grid = element_blank()) +
        scale_color_manual(legt, values=a.col)+
        scale_linetype_manual(legt, values=rep(nLv:1, a.0))+
        scale_shape_manual(legt, values=rep(1:nLv, a.0))+
        annotate("text", x=maxIDR, y=0, label=list('paste(p[a])'), parse=TRUE) +
        geom_hline(yintercept=0.5, color='darkgray') +  
        geom_hline(yintercept=goTh[1], color='tan3') +  
        geom_vline(xintercept=maxIDR, color='gray') +
        theme(legend.position=oc_legp)

      # pp.oc
      
      #----------------------plot for probability coverage----------------------#
      pp.cov<- prob_cov(sample_size=N1, p_a=maxIDR, ab=bay.prior,
                        show_freq=ifelse(is.null(stopAt), 1, stopAt[1]),
                        lab=mNM0)$plot 
      # pp.cov
            
      #----------------------Estimate +/- 95%CI---------------------------#
      #7. plot the estimated values for p vs y #ADR with 95%CI. 
      if( is.numeric(errII.xlim)&length(errII.xlim)==2 ){
        errs.02 <- errs.2[errs.2$numEvts>=min(errII.xlim) & 
                            errs.2$numEvts<=max(errII.xlim),]
      }else{errs.02 <- errs.2}
      Lb0<-round(errs.02[,c('phat','phat.L', 'phat.H')], pdig+2)
      errs.02$mylab<-paste0(Lb0[,1],'(',Lb0[,2], ', ', Lb0[,3], ')')
      p7<-ggplot(data=errs.02, 
                 aes(x=numEvts, y=phat, color=grp, linetype=grp, fill=grp)) + 
        geom_errorbar(aes(ymin=phat.L, ymax=phat.H), 
                      width = 0.3, colour='darkgray',
                      position = position_dodge(0.25)) +
        geom_line( position = position_dodge(0.25)) + 
        theme_bw() + theme(plot.title = element_text(size = 11), 
                           legend.position='bottom')+
        ggtitle("Estimates") + 
        ylab('Posterior 95% CIs') + 
        xlab('Data: y') +
        scale_linetype('') +  
        scale_colour_manual('', values=lab2) +
        scale_fill_manual('', values=lab2) +
        geom_hline(yintercept=maxIDR, linetype='dashed', color='darkgray') +
        geom_text(x=max(errs.02$numEvts, na.rm=T), y=maxIDR, 
                  label=as.character(maxIDR), hjust=0.75, color='darkgray')
      pp.est<-xyL(pp=p7, xl=NULL, yl=NULL, ssa=sa, myLabel=errs.02, 
                  leg='no', addref=FALSE, anno.pos=0.25) 
      
      #----------------------show graphs in layout---------------------------#
      if(!noplot){
        x3loc<-0.3; x4loc <- 0.6
        show(ggdraw() + 
               draw_plot(pp.prput, x=0,     y=.5, width=x3loc, height=.5) + 
               draw_plot(pp.po, x=x3loc, y=.5, width=x4loc-x3loc, height=.5) +
               draw_plot(pp.poput, x=x4loc, y=.5, width=1-x4loc, height=.5) + 
               draw_plot(pp.oc, x=0,     y=0,  width=x3loc, height=.5) + 
               draw_plot(pp.cov, x=x3loc, y=0,  width=x4loc-x3loc, height=.5) + 
               draw_plot(pp.est, x=x4loc, y=0,  width=1-x4loc, height=.5) + 
               draw_plot_label(c('A','B', 'C', 'D', 'E', 'F'), 
                               x=c(0, x3loc+0.01, x4loc+0.01, 0, 
                                   x3loc+0.01, x4loc+0.01), 
                               y=c(1,1,1, 0.5,0.5,0.5), 
                               size=15))
      }
    }else{
      #layout the 2 graphs
      if(!noplot){
        x3loc<-0.3; x4loc <- 0.6
        if(plot.type==2){
          show( ggdraw() + 
                  draw_plot(pp.poput, x=0,     y=.5, width=x3loc, height=.5) + 
                  draw_plot(pp.po, x=x3loc, y=.5, width=x4loc-x3loc, height=.5) +
                  draw_plot(pp.prput, x=x4loc, y=.5, width=1-x4loc, height=.5) +
                  draw_plot_label(c('A','B', 'C'), 
                                  x=c(0, x3loc+0.01, x4loc+0.01), y=c(1,1,1), 
                                  size=15))
        }else{
          show( ggdraw() + 
                  draw_plot(pp.poput, x=0,     y=.5, width=x3loc, height=.5) + 
                  draw_plot(pp.po, x=x3loc, y=.5, width=x4loc-x3loc, height=.5) +
                  draw_plot(pp.prput, x=x4loc, y=.5, width=1-x4loc, height=.5) +
                  draw_plot_label(c('A','B', 'C'), 
                                  x=c(0, x3loc+0.01, x4loc+0.01), y=c(1,1,1), 
                                  size=15)) 
        }
      }
    }
    oth<-errs.2[,!colnames(errs.2)%in%'mylab']
    sel<-oth[oth$numEvts%in%sa, c('grp','numEvts', 'fpr', 'fnr')]
    if(nrow(sel)>0){
      sel$note<-"n/a"
      sel$note[sel$fnr<=goTh[2]]<-'go'
      sel$note[sel$fn>goTh[2]&sel$fpr<goTh[1]]<-'stop'
    }
    if(exists("stop_rule")){
      rule<-rbind( bind_rows(stop_rule))
      r2<-rule[rule[,3]!='Not Applicable',]
      rule2<-NULL
      for(rr in unique(r2$model)){
        r3<-r2[r2$model==rr,]
        rule2<-c(rule2, paste0('Using ', rr, ': '), 
                 paste("    stop if >=", r3[,3], 
                       'events in ', r3[,2], 'subjects.'), '  ')
      }
      rule2<-paste(rule2, collapse = '\n')
    }else{rule<-rule2<-NULL}
    return(list(posterior=psto, oth=oth, sel=sel, stop_rule=rule, rule=rule2))
  }# nzib(plot4=T)$rule
  
  
  #-----------------------function: Print result in table----------------------#
  Dis<-function(x='0.025 (0.01, 0.09)'){
    xx<-unlist(strsplit(gsub(')','',x), split='(', fixed=TRUE))
    xx<-as.numeric(unlist(strsplit(xx, split=',', fixed=TRUE)))
    sum(abs(xx[1]-xx[2:3]))
  }
  Tsh<-function(rr=r1, rg=1:3, bay.lab='', shdig=pdig){ 
    #show estimates
    est<-rr$oth[rr$oth$numEvts%in%rg,]
    if(shdig==0){shdig=3}
    x <- paste0(round(est$phat, shdig), ' (', round(est$phat.L, shdig), 
                ', ', round(est$phat.H, shdig),')')
    #distance
    y <- sapply(x, Dis)
    
    #table out    
    ot<-data.frame(numEvts=est$numEvts, mod=est$grp, est=x, dist=y, pv=est$fpr, 
                   po.lt=est$sup, pop.ut=est$fnr)
    ot<-ot[order(ot$mod, ot$numEvts),]
    ot$mod<-as.character(ot$mod)
    ot$mod[ot$mod=='Bayesian']<-paste('Bayesian',bay.lab)
    return(ot)
  }  
} #end for all source functions


if(F){
  dev.off()

  # main analysis
  #input$text:  1. stopping threshold for number of Adverse Drug Reaction (#IADR)
  #input$text2: 2. max number of samples
  #input$text3: 3. max expected acceptable occurrence rate of ADR
  #input$text4: 4. expected occurrence rate of ADR considered as normal, 
  #               only useful for Binomial
  #input$text5: 5. Hyper parameters for traditional Bayesian
  #input$text6: 6. X-axis range in risk plot
  #input$text7: 7. Y-axis range in risk plot
  #input$text8: 8. X-axis range in type I error plot and Superiority plot
  #input$radio: 9. whether show all 4 available plots
  

  
  out.pdf<-'nzib_sim_n40_p03.pdf'
  out1<-gsub('.pdf', '_1.csv', out.pdf)
  out2<-gsub('.pdf', '_2.csv', out.pdf)
  
  #comparing the simulated results
  pdf(file=out.pdf, width=15, height=7)
  r1<- nzib(stopAt="1,2,3",   #1. stopping threshold
         N1=40,               #2. max number of samples
         maxIDR="0.03",       #3. max expected acceptable occurrence rate
         plot.type=1, 
         bay.prior=c(1, 1),   #5. hyper parameter in beta-binomial
         rsk.xlim="0,5",      #6. x-axis range for risk plot
         errII.xlim="0,5",    #8. x-axis range for type I error plot
         plot4=T,              #9. whether show all 4 available plots
         goTh="0.1",
         ajpm='BH'
    )
  #get the result for Bayesian model only  
  r2<-nzib(stopAt="1,2,3", N1=40, maxIDR="0.03",  plot.type=4, 
         bay.prior=c(alpha=0.03092784, beta=1), 
         rsk.xlim="0,3",  errII.xlim="0,3", plot4=T, ajpm='BH', goTh='0.05')
  r2<-nzib(stopAt="20", N1=40, maxIDR="0.05",  plot.type=4, 
           bay.prior=c(alpha=1, beta=0.5), 
           plot4=T, ajpm='BH', goTh='0.05')
  
  r2<-nzib(stopAt="3", N1=400, maxIDR="0.001",  plot.type=4, 
           bay.prior=c(alpha=0.5, beta=0.5), 
           rsk.xlim="0,3",  errII.xlim="0,3", plot4=T, ajpm='BH', goTh='0.05')
  dev.off()
  
  #print result
  t1<-Tsh(r1, 1:3, '1,1')
  t2<-Tsh(r2, 1:3, '0.5,0.5')
  write.csv(t1, file=out1, row.names=F)
  write.csv(t2, file=out2, row.names=F)
  
  
  #generate the density plot of prior distribution
  pdf(file='prior_density.pdf', width=10, height=6.1)
  bbin(
    prior.a="1,0.5, 0.0051", #1. parameters in p~Beta(a=, b=) prior distribution
    prior.b="1,0.5,1", #2. parameters in p~Beta(a=, b=) prior distribution
    N1=45,            #3. planned sample size
    pa='0.01',           #4. expected average for event rate p
    p.xL='',        #5. xlim for beta distribution
    x.xL=''         #6. xlim for beta-binomial distribution
  )

  
  bbin(
    prior.a="0.01/(1-0.01), 0.05", #1. parameters in p~Beta(a=, b=) prior distribution
    prior.b="1, 0.05", #2. parameters in p~Beta(a=, b=) prior distribution
    N1=45,            #3. planned sample size
    pa='0.01',           #4. expected average for event rate p
    p.xL='',        #5. xlim for beta distribution
    x.xL='',         #6. xlim for beta-binomial distribution
    valLeg=TRUE
  )  

  bbin(
    prior.a="0.01/(1-0.01), 0.05", #1. parameters in p~Beta(a=, b=) prior distribution
    prior.b="1, 0.05", #2. parameters in p~Beta(a=, b=) prior distribution
    N1=45,            #3. planned sample size
    pa='0.01',           #4. expected average for event rate p
    p.xL='',        #5. xlim for beta distribution
    x.xL=''         #6. xlim for beta-binomial distribution
  )
  
  #input$text:  1. parameters in p~Beta(a=, b=) prior distribution
  #input$text2: 2. parameters in p~Beta(a=, b=) prior distribution
  #input$text3: 3. planned sample size
  dev.off()
  
}