demo/Simulation2026_discrete_functions.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

#DIR<- "/data/matejaam/"

#source(paste0(DIR, "bpcp2sample.R"))
#source(paste0(DIR, "kmciFunctions.R"))
#source(paste0(DIR, "delta2samp.R"))


library(bpcp)
packageVersion("bpcp")


createData<-function(n,p,pi.c,prop.treated=0.5){
  pA<-p[1]
  pE<-p[2]
  pN<- 1-pA-pE
  if (pA>1 | pA<0 | pE>1 | pE<0 | pN>1 | pN<0) stop("incorrect p")
  cnt<- rmultinom(1,n,prob=c(pA,pE,pN))
  x<-rep(c(1,2,7),times=cnt)
  x<- sample(x,replace=FALSE)
  n.treated<- round(n*prop.treated)
  n.cntl<- n-n.treated
  group<- c(rep(1,n.cntl),rep(2,n.treated))
  # treatment changes t=2 to t=4 in effected population
  x[group==2 & x==2]<- 4
  cens<- rbinom(n,1,pi.c)
  cens[cens==1]<- 3
  cens[cens==0]<- 6
  y<- pmin(x,cens)
  status<-rep(0,n)
  status[y==x]<- 1
  #data.frame(group=group,x=x,cens=cens,y=y,status=status)
  # a list is faster to use in simulations
  list(group=group,x=x,cens=cens,y=y,status=status)
}


sim<-function(nsim,N,P,PI.C,PROP.TREATED=0.5,TESTTIME=5,
              PARMTYPE="efflogs",NULL.BETA=0,
              alpha=.05,SEED=10311){
  set.seed(SEED)
  rejectLo<-rejectHi<-pvals<-rep(0,nsim)
  CImat<-matrix(NA,nsim,8)
  for (i in 1:nsim){
    x<-createData(n=N,p=P,pi.c=PI.C,prop.treated = PROP.TREATED)
    
    CImat[i,1:2]<-bpcp2samp(x$y,x$status,x$group,TESTTIME,
                            conf.level=1-alpha,parmtype=PARMTYPE,
                            control=bpcp2sampControl(seed=NULL))$conf.int
    CImat[i,3:4]<-bpcp2samp(x$y,x$status,x$group,TESTTIME,
                            conf.level=1-alpha,parmtype=PARMTYPE,
                            midp=TRUE, control=bpcp2sampControl(seed=NULL))$conf.int
    CImat[i,5:6]<-delta2samp(x$y,x$status,x$group,TESTTIME,
                            conf.level=1-alpha,parmtype=PARMTYPE,
                            zero.one.adjustment = TRUE, 
                            method="standard")$conf.int
    CImat[i,7:8]<-delta2samp(x$y,x$status,x$group,TESTTIME,
                             conf.level=1-alpha,parmtype=PARMTYPE,
                             zero.one.adjustment = TRUE, 
                             method="adj_hybrid")$conf.int
    #print(paste0("i=",i))
    #if (CImat[i,1]<0) browser()
    
  }
  dimnames(CImat)[[2]]<-c("meld.lo","meld.hi",
                          "meld.midp.lo","meld.midp.hi",
                          "delta.zo.lo","delta.zo.hi",
                          "delta.badjhy.lo","delta.badjhy.hi")
  
  summarizeCI<-function(ci,nullBeta=NULL.BETA){
    N<- nrow(ci)
    N.lo<- sum(!is.na(ci[,1]))
    N.hi<- sum(!is.na(ci[,2]))
    prop.Res.lo<- N.lo/N
    prop.Res.hi<- N.hi/N
    reject.lo<- sum(!is.na(ci[,1]) & ci[,1]>nullBeta)/N.lo
    reject.hi<- sum(!is.na(ci[,1]) & ci[,2]<nullBeta)/N.hi
    rejectTotal.lo<- sum(!is.na(ci[,1]) & ci[,1]>nullBeta)/N
    rejectTotal.hi<- sum(!is.na(ci[,1]) & ci[,2]<nullBeta)/N
    out<-c(N=N,N.lo=N.lo,N.hi=N.hi,prop.Res.lo=prop.Res.lo,
           prop.Res.hi=prop.Res.hi,reject.lo=reject.lo,
           reject.hi=reject.hi,rejectTotal.lo=rejectTotal.lo,
           rejectTotal.hi=rejectTotal.hi)
    out
  }
  
  sout<-matrix(NA,ncol(CImat)/2,9)
  
  s1<-summarizeCI(CImat[,1:2])
  sout[1,]<- s1 
  sout[2,]<- summarizeCI(CImat[,3:4])
  sout[3,]<- summarizeCI(CImat[,5:6])
  sout[4,]<- summarizeCI(CImat[,7:8])

  dimnames(sout)<- list(
    c("meld","meld.midp","delta.zo","delta.badjhy"),
    names(s1))
  
  list(sout=sout,CImat=CImat)
}


pE <- seq(0, 0.4, 0.01)
p<-data.frame(n=rep(600,length(pE)),
              pA=rep(.5,length(pE)),
              #pE=c(.1,.2,.3,.4),
              pE=pE,
              pi.c=rep(.9,length(pE)),
              prop.treated=rep(.5,length(pE)),
              parmtype=rep("efflogs",length(pE)))


do_sim <- function(NSIM=1000, seed, output1) {
  
  out<- array(NA,c(NSIM,8,nrow(p)),
              dimnames=list(c(1:NSIM),
                            c("meld.lo","meld.hi",
                              "meld.midp.lo","meld.midp.hi",
                              "delta.zo.lo","delta.zo.hi",
                              "delta.badjhy.lo","delta.badjhy.hi"),
                            c(paste0("Scenario ",1:nrow(p))))
  )
  
  
  for (i in 1:nrow(p)){
    x<-sim(NSIM,
           N=p$n[i],P=c(p$pA[i],p$pE[i]),PI.C=p$pi.c[i],
           PROP.TREATED=p$prop.treated[i],
           PARMTYPE=p$parmtype[i],
           SEED=seed)
    out[,,i]<- x$CImat
  }
  
  saveRDS(out, output1)
  
}



# Warning, takes many hours to run the next line
#do_sim(NSIM=1000,seed=1243, output1="./simResults/efflogs2026_discrete_nsim1000.rds")