bpcp
Beta Product Confidence Procedure for Right Censored Data

Package index
Search the bpcp package
Vignettes
Functions
196
Source code
26
Man pages
25
Home
/
CRAN
/
bpcp
/
demo/Calculation2025_no_censoring.R

demo/Calculation2025_no_censoring.R
In bpcp: Beta Product Confidence Procedure for Right Censored Data 

#' ---
#' title: "Calculation of Type I Error Rate for No Censoring Case"
#' author: "Michael Fay"
#' date: "`r Sys.Date()`"
#' output: html_document
#' ---
#' 
## ----setup, include=FALSE------------------------------------------------------------------------------
# since MIDP takes the longest, for checking plots and code do not always run it
DOMIDP<- TRUE
#knitr::opts_chunk$set(echo = FALSE)

#' 
#' 
## ------------------------------------------------------------------------------------------------------
library(exact2x2)

NGRID<- 50

#library(bpcp)

#DIR<- "H:/My Documents/methods/bpcp2samp/R/"
#DIR<- "C:/Users/matejaam/OneDrive - National Institutes of Health/_HDrive/Projects/bpcp/R/"
#source(paste0(DIR,"delta2samp.R"))
#source(paste0(DIR,"kmciFunctions.R"))
#source(paste0(DIR,"bpcp2sample.R"))
#source(paste0(DIR,"mdiffmedian.test.R"))

#' 
#' ## Why We Don't Need a Simulation
#' 
#' For the no censoring case, the two-sample problem for comparing $S_1(t)$ and $S_2(t)$ 
#' reduces to the binomial problem for comparing two binomial proportions. 
#' So we only need to calculate the confidence interval for each possible 
#' response (i.e., $X_1=0,\ldots,n_1$ and $X_2=0,\ldots,n_2$). There are 
#' $(n_1+1)(n_2+1)$ possible responses. Then we can easily check any pair of 
#' true proportions to see if the confidence interval covers it or not. 
#' Because the methods we are using a symmetric, and our two-sided confidence 
#' intervals just use the intersection of the two one-sided intervals,
#' we only need to check one side (e.g., whether the upper limit is less than the true 
#' parameter or not).
#' 
#' 
#' # Differences
#' 
#' ## Differences Upper Limits
#' 
#' We start with testing the upper limit when the parmtype='difference'
#' 
#' 
## ------------------------------------------------------------------------------------------------------
# edited checkCoverage from supplement to Fay, et al (2015, Biometrics)
checkCoverage<-function(n1,n2,ngrid=10^2,cifunc,betafunc,limit=c("upper","lower")){
   # cifunc must give an object x with 
   # an element x$conf.int that is a vector of length 2 
    limit<-match.arg(limit)
    CLimit<-Est<-matrix(NA,n1+1,n2+1,dimnames=list(paste("X1=",0:n1,sep=""),paste("X2=",0:n2,sep="")))
    for (i in 0:n1){
        for (j in 0:n2){
            #print(paste0("i=",i," j=",j))
            # only want upper coverage, so alternative=less
            ci<- cifunc(i,n1,j,n2)$conf.int
            if (limit=="upper"){
              CLimit[i+1,j+1]<- ci[2]             
            } else if (limit=="lower"){
              CLimit[i+1,j+1]<- ci[1]                 
            }
 
            Est[i+1,j+1]<- cifunc(i,n1,j,n2)$estimate
        }
    }



    CoverageFunc<-function(p1,p2,Limit=limit,N1=n1,N2=n2,b=betafunc){
        f1<-matrix(dbinom(0:N1,N1,p1),N1+1,1)
        f2<-matrix(dbinom(0:N2,N2,p2),1,N2+1)
        f<-f1 %*% f2
        #dimnames(f)<-list(paste("X1=",0:N1,sep=""),
        #    paste("X2=",0:N2,sep=""))
        if (Limit=="upper"){
           out<-sum(f[b(p1,p2)<=CLimit], na.rm=T)         
        } else if (Limit=="lower"){
           out<-sum(f[b(p1,p2)>=CLimit], na.rm=T)  
        } else stop("limit must be either 'upper' or 'lower' ")
 
    }
    #T1<-c(0,.0001,1:100/101,.9999,1)
    #P1<-c(1/ngrid^(5:2),1:(ngrid-1)/ngrid,1-1/ngrid^(2:5))
    P1<- 0:ngrid/ngrid
    Coverage<-Beta<-matrix(NA,length(P1),length(P1),
                     dimnames=list(paste0("S1(t)=",P1),paste0("S2(t)=",P1)))
    for (i in 1:length(P1)){
        for (j in 1:length(P1)){
            Beta[i,j]<- betafunc(P1[i],P1[j])
            Coverage[i,j]<-CoverageFunc(P1[i],P1[j])
       }
    }
    list(Beta=Beta,Coverage=Coverage,CLimit=CLimit,Estimate=Est)
}




createTime<-function(x,n,testtime=1){
  if (x==n){
    time<- testtime + 1:n
  } else if (x==0){
    time<- testtime* (c(1:n)/(n+1))
  } else {
    time<- c( testtime*c((x+1):n)/(n+1), 
              testtime+ 1:x)
  }
  time
}


CIfuncDelta<-function(x1,n1,x2,n2,Conf.level=0.975,Parmtype="difference",
                      Alternative="less", Method="standard", Zero.one.adjustment=TRUE){
  # make survival data with no censoring such that 
  # Kaplan-Meier survival estimates at t, 
  # say S1(t) and S2(t) are
  #
  # S1(t) = x1/n1 and S2(t)=x2/n2
  #
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  Testtime<- 1 
  if (Testtime!=1) stop("testtime must be 1")

  delta2samp(time=Time,
            status=Status,
            group=Group, 
            testtime=Testtime,
            parmtype=Parmtype,
            alternative=Alternative, 
            conf.level=Conf.level, 
            method=Method, 
            zero.one.adjustment = Zero.one.adjustment)
}


CIfuncBPCP<-function(x1,n1,x2,n2,Conf.level=0.975,Midp=FALSE, 
                     Parmtype="difference",Alternative="less"){
  # make survival data with no censoring such that 
  # Kaplan-Meier survival estimates at t, 
  # say S1(t) and S2(t) are
  #
  # S1(t) = x1/n1 and S2(t)=x2/n2
  #
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  Testtime<- 1 
  if (Testtime!=1) stop("testtime must be 1")

  bpcp2samp(time=Time,
            status=Status,
            group=Group, 
            testtime=Testtime,
            parmtype=Parmtype,
            alternative=Alternative, 
            conf.level=Conf.level, 
            midp=Midp)
}

#' 
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"difference"
PARMTYPE.LABEL<- "S2(t)-S1(t)"
Betafunc<-function(S1,S2){ S2-S1 }
transBeta<- function(beta){ beta }
LIMIT<-"upper"
ALT<- "less"
CL<- 0.975
N1<- 12
N2<- 24


CIfunc1<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
  CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=FALSE)
}
CIfunc3<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
  #               midp=TRUE)
    CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method="adj_hybrid", Zero.one.adjustment=TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
plotCoverage<-function(ptype){
  par(mfrow=c(2,2))
  
  if (ptype=="difference"){
    tr<- function(x){ x }
    axesFunc<-function(){
      axis(1)
      axis(2, at=c(0.8, 0.85, 0.9, 0.95, 1), las=2)
      box()
    }
  }  else if (ptype=="ratio"){
    tr<- transBeta
    axesFunc<-function(){
      AT<- c(0.01,0.1,.25,.5,1,2,4,10,100)
      axis(1, at=transBeta(AT),labels=AT)
      axis(2, las=2)
      box()
    }
  } else if (ptype=="efflogs"){
    tr<- transBeta
    axesFunc<-function(){
      AT<- c(-Inf,-4,-1,0,0.5,0.8,1)
      axis(1, at=transBeta(AT),labels=AT)
      axis(2, las=2)
      box()
    }
  }
  
  YLAB<- paste0("Coverage (",LIMIT," CL)")
  plot(tr(as.vector(out1$Beta)),as.vector(out1$Coverage),ylim=c(0.8,1),
       xlab=PARMTYPE.LABEL,ylab=YLAB,
       main="Melding",
       axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
  axesFunc()
  lines(c(-1,1),c(0.975,0.975),lty=2,col="red")
  
  
  plot(tr(as.vector(out2$Beta)),as.vector(out2$Coverage),ylim=c(0.8,1),
       xlab=PARMTYPE.LABEL,ylab=YLAB,
       main="Delta Method (Greenwood)",
       axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
  axesFunc()
  lines(c(-1,1),c(0.975,0.975),lty=2,col="red")
  
  
  plot(tr(as.vector(out3$Beta)),as.vector(out3$Coverage),ylim=c(0.8,1),
       xlab=PARMTYPE.LABEL,ylab=YLAB,
       main="Melding (mid-p version)",
       axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
  axesFunc()
  lines(c(-1,1),c(0.975,0.975),lty=2,col="red")
  
  
  plot(tr(as.vector(out4$Beta)),as.vector(out4$Coverage),ylim=c(0.8,1),
       xlab=PARMTYPE.LABEL,ylab=YLAB,
       main="Delta Method (Borkowf Adjusted Hybrid)",
       axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
  axesFunc()
  lines(c(-1,1),c(0.975,0.975),lty=2,col="red")
}


pdf("DiffNoCensUp_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
plotEstCL<-function(ptype){
  par(mfrow=c(2,2))   
  if (ptype=="difference"){
    tr<- function(x){ x }
    axesFunc<-function(){
      axis(1)
      axis(2, las=2)
      box()
    }
  }  else if (ptype=="ratio"){
    tr<- transBeta
    axesFunc<-function(){
      AT<- c(0.01,0.1,.25,.5,1,2,4,10,100)
      axis(1, at=transBeta(AT),labels=AT)
      axis(2, at=transBeta(AT),labels=AT, las=2)
      box()
    }
  } else if (ptype=="efflogs"){
    tr<- transBeta
    axesFunc<-function(){
      AT<- c(-Inf,-4,-1,0,0.5,0.8,1)
      axis(1, at=transBeta(AT),labels=AT)
      axis(2, at=transBeta(AT),labels=AT, las=2)
      box()
    }
  }
  
  YLAB<- paste0(round(100*CL,1),"% ",LIMIT," Limit") 
  XLAB<- paste0("Estimate of ",PARMTYPE.LABEL)

plot(tr(as.vector(out1$Estimate)),tr(as.vector(out1$CLimit)),xlim=c(-1,1), ylim=c(-1,1),
     xlab=XLAB,ylab=YLAB,
     main="Melding",
     axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
axesFunc()
lines(c(-1,1),c(-1,1),lty=2)

plot(tr(as.vector(out2$Estimate)),tr(as.vector(out2$CLimit)),xlim=c(-1,1), ylim=c(-1,1),
     xlab=XLAB,ylab=YLAB,
     main="Delta Method (Standard)",
     axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
axesFunc()
lines(c(-1,1),c(-1,1),lty=2)

plot(tr(as.vector(out3$Estimate)),tr(as.vector(out3$CLimit)),xlim=c(-1,1), ylim=c(-1,1),
     xlab=XLAB,ylab=YLAB,
     main="Melding (mid-p version)",
     axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
axesFunc()
lines(c(-1,1),c(-1,1),lty=2)

plot(tr(as.vector(out4$Estimate)),tr(as.vector(out4$CLimit)),xlim=c(-1,1), ylim=c(-1,1),
     xlab=XLAB,ylab=YLAB,
     main="Delta Method (Borkowf Adjusted Hybrid)",
     axes=FALSE,col=adjustcolor("black",alpha.f=0.1),pch=16)
axesFunc()
lines(c(-1,1),c(-1,1),lty=2)
}


pdf("DiffNoCensUp_Limits.pdf", width = 12, height = 8)
plotEstCL("difference")
dev.off()

#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' ## Differences Lower Limits 
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"difference"
Betafunc<-function(S1,S2){ S2-S1 }
transBeta<- function(beta){ beta }
LIMIT<-"lower"
ALT<- "greater"
CL<- 0.975
N2<- 12
N1<- 24


CIfunc1<-function(x1,n1,x2,n2){
  binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
}
CIfunc3<-function(x1,n1,x2,n2){
  binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
                 midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method="adj_hybrid", Zero.one.adjustment=TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
pdf("DiffNoCensLo_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
pdf("DiffNoCensLo_Limits.pdf", width = 12, height = 8)
plotEstCL(PARMTYPE)
dev.off()

#' 
#' 
#' 
#' 
#' 
#' 
#' # Ratios
#' 
#' ## Ratios Upper Limits
#' 
#' We start with testing the upper limit when the parmtype='ratio'
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"ratio"
PARMTYPE.LABEL<-"S2(t)/S1(t)"
Betafunc<-function(S1,S2){ S2/S1 }
transBeta<- function(beta){ (beta-1)/(beta+1) }
LIMIT<-"upper"
ALT<- "less"
CL<- 0.975
N1<- 12
N2<- 24



CIfunc1<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
  CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=FALSE)
}
CIfunc3<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
  #               midp=TRUE)
    CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method = "standard", Zero.one.adjustment = TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method = "adj_hybrid", Zero.one.adjustment = TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("RatioNoCensUp_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("RatioNoCensUp_Limits.pdf", width = 12, height = 8)
plotEstCL(PARMTYPE)
dev.off()

#' 
#' 
#' ## Ratios Lower Limits 
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"ratio"
PARMTYPE.LABEL<-"S2(t)/S1(t)"
Betafunc<-function(S1,S2){ S2/S1 }
transBeta<-function(beta){ (beta-1)/(beta+1) }
LIMIT<-"lower"
ALT<- "greater"
CL<- 0.975
N2<- 12
N1<- 24


CIfunc1<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
  CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=FALSE)
}
CIfunc3<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
  #               midp=TRUE)
    CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method = "standard", Zero.one.adjustment = TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method = "adj_hybrid", Zero.one.adjustment = TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
#' 
#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
pdf("RatioNoCensLo_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("RatioNoCensLo_Limits.pdf", width = 12, height = 8)
plotEstCL(PARMTYPE)
dev.off()

#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' # Efficacy using Log Survival Ratios
#' 
#' ## Efficacy Log(S) Ratios Upper Limits
#' 
#' We start with testing the upper limit when the parmtype='efflogs'
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"efflogs"
PARMTYPE.LABEL<-"1-log(S2(t))/log(S1(t))"
Betafunc<-function(S1,S2){ 1-log(S2)/log(S1) }
transBeta<- function(beta){ beta/(2-beta) }
LIMIT<-"upper"
ALT<- "less"
CL<- 0.975
N1<- 12
N2<- 24



CIfunc1<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
  CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=FALSE)
}
CIfunc3<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
  #               midp=TRUE)
    CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method = "standard", Zero.one.adjustment = TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method = "adj_hybrid", Zero.one.adjustment = TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("EfflogsNoCensUp_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("EfflogsNoCensUp_Limits.pdf", width = 12, height = 8)
plotEstCL(PARMTYPE)
dev.off()

#' 
#' 
#' ## Efficacy Log(S) Ratios Lower Limits
#' 
## ------------------------------------------------------------------------------------------------------
PARMTYPE<-"efflogs"
PARMTYPE.LABEL<-"1-log(S2(t))/log(S1(t))"
Betafunc<-function(S1,S2){ 1-log(S2)/log(S1) }
transBeta<- function(beta){ beta/(2-beta) }
LIMIT<-"lower"
ALT<- "greater"
CL<- 0.975
N1<- 12
N2<- 24


CIfunc1<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,alternative=ALT,parmtype=PARMTYPE)
  CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=FALSE)
}
CIfunc3<-function(x1,n1,x2,n2){
  #binomMeld.test(x1,n1,x2,n2,conf.level=CL,parmtype=PARMTYPE,alternative=ALT,
  #               midp=TRUE)
    CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
             Midp=TRUE)
}
CIfunc2<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
              Alternative=ALT, Method = "standard", Zero.one.adjustment = TRUE)
}
CIfunc4<-function(x1,n1,x2,n2){
  CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,
              Parmtype=PARMTYPE,Alternative=ALT, 
              Method = "adj_hybrid", Zero.one.adjustment = TRUE)
}

t0<- proc.time()
out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
out2<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2,betafunc=Betafunc,limit=LIMIT)
if (DOMIDP){
  out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
} else out3<-out1
out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0

#' 
## ----echo=FALSE, fig.height=8--------------------------------------------------------------------------
pdf("EfflogsNoCensLo_Coverage.pdf", width = 12, height = 8)
plotCoverage(PARMTYPE)
dev.off()

#' 
## ----fig.height=8--------------------------------------------------------------------------------------
pdf("EfflogsNoCensLo_Limits.pdf", width = 12, height = 8)
plotEstCL(PARMTYPE)
dev.off()

#' 
#'

Try the bpcp package in your browser

Any scripts or data that you put into this service are public.

bpcp documentation built on July 21, 2026, 5:08 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close