demo/Sim2026_functions.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

library(survival)
library(bpcp)
if (packageVersion("bpcp")=="1.5.1") stop("use bpcp 1.5.2 or later")

# pCENS for plot
pCENS<-0.40

## -------------------------------------------------------------------------------------------------------
# function for getting the true survival functions, censoring function,
# and proportions at risk in the two groups (and overall)
getFuncs_SurvAtRisk<-function(S1,S2,propGroup1=0.5,censor=c(0,2), propCens=pCENS){
  # S1 and S2 are the survival functions
  G<-function(t,pc=propCens,urange=censor){
     1-pc*punif(t,min=urange[1],max=urange[2])    
  }
  prop.atrisk<-function(t){ (propGroup1*S1(t)+(1-propGroup1)*S2(t))*G(t)  }
  prop1.atrisk<-function(t){ S1(t)*G(t)  }
  prop2.atrisk<-function(t){ S2(t)*G(t)  }
  out<-list(S1=S1,S2=S2,G=G,prop.atrisk=prop.atrisk,
            prop1.atrisk=prop1.atrisk,
            prop2.atrisk=prop2.atrisk)
  out
}


## -------------------------------------------------------------------------------------------------------
# Scenario 1
Sc1Rate1<-Sc1Rate2<- 0.2
S1<-function(t,rate1=Sc1Rate1){
  1-pexp(t,rate=rate1)
}
S2<-function(t,rate2=Sc1Rate2){
  1-pexp(t,rate=rate2)
}
f1<-getFuncs_SurvAtRisk(S1,S2,propGroup1=0.5,censor=c(0,2), propCens=pCENS)

# Scenario 1, failure time generating functions 
Sc1_Rtime1<-function(n,rate1=Sc1Rate1){ rexp(n,rate=rate1) }
Sc1_Rtime2<-function(n,rate2=Sc1Rate2){ rexp(n,rate=rate2) }






## -------------------------------------------------------------------------------------------------------
#Scenario 2
Sc2Rate1<- 0.1
Sc2Rate2<- 0.01
S1<-function(t,rate1=Sc2Rate1){
  1-pexp(t,rate=rate1)
}
S2<-function(t,rate2=Sc2Rate2){
  1-pexp(t,rate=rate2)
}
f2<-getFuncs_SurvAtRisk(S1,S2,propGroup1=0.5,censor=c(0,2), propCens=pCENS)

# Scenario 2, failure time generating functions 
Sc2_Rtime1<-function(n,rate1=Sc2Rate1){ rexp(n,rate=rate1) }
Sc2_Rtime2<-function(n,rate2=Sc2Rate2){ rexp(n,rate=rate2) }


## -------------------------------------------------------------------------------------------------------
# Scenario 3
Sc3Rate1<- 0.5
Sc3propB<-0.30
Sc3RateA<-0
Sc3RateB<-10

 S1<-function(t,rate1=Sc3Rate1){
  1-pexp(t,rate=rate1)
}
S2<-function(t,propB=Sc3propB,rateA=Sc3RateA,rateB=Sc3RateB){
    propB*(1-pexp(t,rate=rateB)) + (1-propB)*(1-pexp(t,rate=rateA))
}
f3<-getFuncs_SurvAtRisk(S1,S2,propGroup1=0.5,censor=c(0,2), propCens=pCENS)

# Scenario 3, failure time generating functions 
Sc3_Rtime1<-function(n,rate1=0.5){  rexp(n,rate=rate1)  }
Sc3_Rtime2<-function(n,propB=Sc3propB,rateA=Sc3RateA,rateB=Sc3RateB){
  nB<- rbinom(1,n,propB)
  nA<- n-nB
  timeB<-rexp(nB,rate=rateB)
  if (rateA==0){
    timeA<- rep(Inf,nA)
  } else {
    timeA<-rexp(nA,rate=rateA)    
  }

  time<- sample(c(timeA,timeB),replace=FALSE)
  time
}

# check
#Sc3_Rtime2(100)


## -------------------------------------------------------------------------------------------------------
# Scenario 4
Sc4propB<-0.40
Sc4RateA<-0.10
Sc4a1<-1
Sc4b1<-4
Sc4a2<-4
Sc4b2<-1



  S1<-function(t,propB=Sc4propB,rateA=Sc4RateA,a=Sc4a1,b=Sc4b1){
    propB*(1-pbeta(2*t,a,b)) + (1-propB)*(1-pexp(t,rate=rateA))
  }
  S2<-function(t,propB=Sc4propB,rateA=Sc4RateA,a=Sc4a2,b=Sc4b2){
    propB*(1-pbeta(2*(t-0.5),a,b)) + (1-propB)*(1-pexp(t,rate=rateA))
  }


f4<-getFuncs_SurvAtRisk(S1,S2,propGroup1=0.5,censor=c(0,2), propCens=pCENS)

# Scenario 4, failure time generating functions 
Sc4_Rtime1<-function(n,propB=Sc4propB,rateA=Sc4RateA,a=Sc4a1,b=Sc4b1){
  nB<- rbinom(1,n,propB)
  nA<- n-nB
  if (rateA==0){
    timeA<- rep(Inf,nA)
  } else {
    timeA<-rexp(nA,rate=rateA)    
  }
  timeB<- rbeta(nB,a,b)/2
  time<- sample(c(timeA,timeB),replace=FALSE)
  time
}
Sc4_Rtime2<-function(n,propB=Sc4propB,rateA=Sc4RateA,a=Sc4a2,b=Sc4b2){
  nB<- rbinom(1,n,propB)
  nA<- n-nB
  if (rateA==0){
    timeA<- rep(Inf,nA)
  } else {
    timeA<-rexp(nA,rate=rateA)    
  }
  timeB<- 0.5+rbeta(nB,a,b)/2
  time<- sample(c(timeA,timeB),replace=FALSE)
  time
}


## -------------------------------------------------------------------------------------------------------
plotSimGenerating<- function(f,TITLE=""){
  par(mar=c(4.1,4.1,2.1,4.1))

  plot(c(0,2),c(-0.4,1),type="n",xlab="t",ylab="",
       axes=FALSE,
       main=TITLE,xaxs="i",yaxs="i",las=2)
  axis(1)
  AT<- 0:5/5
  axis(2,at=c(0,.2,.4,.6,.8,1),labels=AT,cex=0.5,las=2)
  axis(2,at=c(-0.2,-0.4),labels=c("Pct","at risk"),las=2,cex=0.5,tick=FALSE)
  axis(2,at=0.5,labels="S(t)",line=1.6,las=2,cex=0.5,tick=FALSE)
  axis(4,at=c(0,-0.2,-0.4),labels=c("100%","50%","0%"),las=2)
  
  tt<- seq(from=0,to=2,length.out=100)
  lines(tt,f$S1(tt),lwd=4,col="gray")
  lines(tt,f$S2(tt),lwd=2,lty=2)

  # draw lines for comparison landmark times
  X0<-c(0.8,1,1.2,1.4,1.8)
  Y0<- rep(0,length(X0))
  
  segments(x0=X0,y0=Y0,x1=X0,y1=Y0+1,lty=2)
  
  # create box for percent of sample at risk
  lines(c(-10,2),c(0,0))
  lines(c(-10,2),c(1,1))
  lines(c(2,2),c(1,0))
  polygon(c(0,2,2,0,0),c(-0.4,-0.4,0,0,-0.4),col=gray(.9))
  #polygon(c(0,tt,2,2,0),-0.4+c(0,0.4*f$prop.atrisk(tt),0.4*f$S.atrisk(2),0,0),col="blue")
  lines(tt,-0.4+0.4*f$prop1.atrisk(tt),lwd=4,col=gray(.7))
  lines(tt,-0.4+0.4*f$prop2.atrisk(tt),lwd=2,lty=2,col="black")
  
  
  
}


## -------------------------------------------------------------------------------------------------------
par(mfrow=c(2,2))
plotSimGenerating(f1,"Scenario 1")
plotSimGenerating(f2,"Scenario 2")
plotSimGenerating(f3,"Scenario 3")
plotSimGenerating(f4,"Scenario 4")


# create pdf for paper
# dev.print(pdf,file="Sim2026_4Scenarios.pdf")


## -------------------------------------------------------------------------------------------------------

########################################################
#
#  Functions for simulations
#
###########################################################
# get statistics from output of one method for simulation table
getstats<-function(out){
  c(est=unname(out$estimate),
    lo=unname(out$conf.int[1]),
    hi=unname(out$conf.int[2]),
    p=unname(out$p.value))
}

# simulation function to create simulation table
sim<-function(nsim,TESTTIME,n1,n2,pCens=0.30,simData,funcs=f,domidp=TRUE){
  parm<- funcs$S2(TESTTIME) - funcs$S1(TESTTIME)
  outTable<-matrix(rep(NA,NSIM*4*5),NSIM,4*5)
  for (i in 1:nsim){
    d<-simData(n1,n2,pCens)
    bout<- bpcp2samp(d$time,d$status,d$grp,testtime=TESTTIME,
                     control=bpcp2sampControl(seed=NULL),conf.level=0.95,
                     parmtype=PARMTYPE)
    if (domidp){
      bout_midp<-bpcp2samp(d$time,d$status,d$grp,testtime=TESTTIME,
                           control=bpcp2sampControl(seed=NULL),
                              conf.level=0.95,midp=TRUE,
                           parmtype=PARMTYPE)
    } else bout_midp<-list(estimate=0,conf.int=c(-Inf,-Inf),p.value=1)
    dout<-  delta2samp(d$time,d$status,d$grp,testtime=TESTTIME,conf.level=0.95, 
                       zero.one.adjustment=TRUE,method="standard",
                       parmtype=PARMTYPE)
    dout_rh<-  delta2samp(d$time,d$status,d$grp,testtime=TESTTIME,conf.level=0.95, 
                          zero.one.adjustment=TRUE,method="reg_hybrid",
                          parmtype=PARMTYPE)
    dout_ah<-  delta2samp(d$time,d$status,d$grp,testtime=TESTTIME,conf.level=0.95, 
                          zero.one.adjustment=TRUE,method="adj_hybrid",
                          parmtype=PARMTYPE)
    
    
    
    tableRow<-c(getstats(bout),getstats(bout_midp),getstats(dout),
                getstats(dout_rh),getstats(dout_ah))
    outTable[i,]<- tableRow
    #if (is.na(tableRow["d_p"])) browser()
  }
  
  dimnames(outTable)[[2]]<-paste0(rep(c("b_","bmp_","d_","drh_","dah_"),each=4),
                                  names(tableRow))
  outTable
}


# transform beta so that beta is on the -1 to 1 scale 
# and the the CI width masks more sense 
if (PARMTYPE=="difference"){
  beta_trans<-function(beta){ beta }
  inv_bt<-function(trans){ trans }
} else if (PARMTYPE=="efflogs"){
  beta_trans<-function(beta){
    out<-beta/(2-beta)
    out[beta==-Inf]<- -1
    out
  }
  inv_bt<-function(trans){
    beta<- 2*trans/(1+trans)
    beta[trans==-1]<- -Inf
    beta
  }
}
# CHECK TRANSFORMATION FOR efflogs
#beta_trans(c(-Inf,-4,0,.5,1))



# function for summarizing the simulation table
simSummary<-function(tab,beta=0,alpha=0.05,betaEqual=BETAEQUAL){
  nsim<- nrow(tab)
  # if KM1(t)=KM2(t)=1 then d_p=NaN
  nNaN_d_p<- sum(is.na(tab[,"d_p"]))
  c(meanTransEst=mean(beta_trans(tab[,"b_est"])),
    invMTE=inv_bt(mean(beta_trans(tab[,"b_est"]))),
    propNA= nNaN_d_p/nsim,
    b_err_lo=sum(tab[,"b_lo"]>beta)/nsim,
    b_err_hi=sum(tab[,"b_hi"]<beta)/nsim,
    b_width=mean(beta_trans(tab[,"b_hi"])-beta_trans(tab[,"b_lo"])),
    b_power_gr=sum(tab[,"b_lo"]>betaEqual)/nsim,
    b_power_less=sum(tab[,"b_hi"]<betaEqual)/nsim,
    bmp_err_lo=sum(tab[,"bmp_lo"]>beta)/nsim,
    bmp_err_hi=sum(tab[,"bmp_hi"]<beta)/nsim,
    bmp_width=mean(beta_trans(tab[,"bmp_hi"])-beta_trans(tab[,"bmp_lo"])),
    bmp_power_gr=sum(tab[,"bmp_lo"]>betaEqual)/nsim,
    bmp_power_less=sum(tab[,"bmp_hi"]<betaEqual)/nsim,
    d_err_lo=sum(tab[,"d_lo"]>beta)/nsim,
    d_err_hi=sum(tab[,"d_hi"]<beta)/nsim,
    d_width=mean(beta_trans(tab[,"d_hi"])-beta_trans(tab[,"d_lo"])),
    d_power_gr=sum(tab[,"d_lo"]>betaEqual,na.rm=TRUE)/(nsim-nNaN_d_p),
    d_power_less=sum(tab[,"d_hi"]<betaEqual,na.rm=TRUE)/(nsim-nNaN_d_p),
    drh_err_lo=sum(tab[,"drh_lo"]>beta)/nsim,
    drh_err_hi=sum(tab[,"drh_hi"]<beta)/nsim,
    drh_width=mean(beta_trans(tab[,"drh_hi"])-beta_trans(tab[,"drh_lo"])),
    drh_power_gr=sum(tab[,"drh_lo"]>betaEqual,na.rm=TRUE)/(nsim-nNaN_d_p),
    drh_power_less=sum(tab[,"drh_hi"]<betaEqual,na.rm=TRUE)/(nsim-nNaN_d_p),
    dah_err_lo=sum(tab[,"dah_lo"]>beta)/nsim,
    dah_err_hi=sum(tab[,"dah_hi"]<beta)/nsim,
    dah_width=mean(beta_trans(tab[,"dah_hi"])-beta_trans(tab[,"dah_lo"])),
    dah_power_gr=sum(tab[,"dah_lo"]>betaEqual)/nsim,
    dah_power_less=sum(tab[,"dah_hi"]<betaEqual)/nsim
    )
}

#simSummary(outTable)


## -------------------------------------------------------------------------------------------------------
# function for generating simulated data
# all Scenarios have the same independent 
# censoring, so for each Scenario we just need to change the
# rtime1 and rtime2 
# (functions to create random failure times in the two groups)
createSimData<-function(n1,n2,rtime1,rtime2,censor=c(0,2), 
                    propCens=pCENS, maxTime=2){
  # maxTime+0.001 = censor time for all those not "censored"
  grp<-c(rep(1,n1),rep(2,n2))
  y<- c(rtime1(n1),rtime2(n2))
  ncens<-rbinom(1,n1+n2,propCens)
  # create cens times, then end of study times
  cens<-c(runif(ncens,min=censor[1],max=censor[2]),
          rep(maxTime+ 0.001,n1+n2-ncens))
  # randomly choose the cens or those that reach the end of the study
  cens<- sample(cens,replace=FALSE)        
  # observed failure time or censor time
  time<- pmin(y,cens)
  # status=1 is failure, status=0 is censored
  status<- ifelse(time==y,1,0)
  out<-data.frame(time=time,status=status,grp=grp)
  out
}


## -------------------------------------------------------------------------------------------------------
## Simulation for Scenario 1
simScen1<-function(N1,N2,propCens=pCENS,RTIME1=Sc1_Rtime1,RTIME2=Sc1_Rtime2){ 
  createSimData(N1,N2,RTIME1,RTIME2,propCens=pCENS) 
}



# test by plotting Kaplan-Meier with a large sample size
d1<-simScen1(1e4,1e4,0.30)
s<-survfit(Surv(time,status)~grp,data=d1)
plot(s)


## -------------------------------------------------------------------------------------------------------
## Simulation for Scenario 2

simScen2<-function(N1,N2,propCens=pCENS,RTIME1=Sc2_Rtime1,RTIME2=Sc2_Rtime2){ 
  createSimData(N1,N2,RTIME1,RTIME2,propCens=pCENS) 
}


# test by plotting Kaplan-Meier with a large sample size
d<-simScen2(1e4,1e4,0.30)
s<-survfit(Surv(time,status)~grp,data=d)
plot(s)



## -------------------------------------------------------------------------------------------------------
## Simulation for Scenario 3

simScen3<-function(N1,N2,propCens=pCENS,RTIME1=Sc3_Rtime1,RTIME2=Sc3_Rtime2){ 
  createSimData(N1,N2,RTIME1,RTIME2,propCens=pCENS) 
}


# test by plotting Kaplan-Meier with a large sample size
d<-simScen3(1e4,1e4,0.30)
s<-survfit(Surv(time,status)~grp,data=d)
plot(s)



## -------------------------------------------------------------------------------------------------------
## Simulation for Scenario 4

simScen4<-function(N1,N2,propCens=pCENS,RTIME1=Sc4_Rtime1,RTIME2=Sc4_Rtime2){ 
  createSimData(N1,N2,RTIME1,RTIME2,propCens=pCENS) 
}


# test by plotting Kaplan-Meier with a large sample size
d<-simScen4(1e4,1e4,0.30)
s<-survfit(Surv(time,status)~grp,data=d)
plot(s)