Nothing
demo/Calculation2026_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------------------------------------------------------------------------------
DOMIDP<- TRUE
#knitr::opts_chunk$set(echo = FALSE)
#'
#'
## ------------------------------------------------------------------------------------------------------
library(exact2x2)
library(bpcp)
#library(NoSleepR)
#nosleep_on()
NGRID<- 100
St<- 0:NGRID/NGRID
#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){
# create times so that length(time)=n and sum(time>testtime)=x
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))
# because of the way we create time (i.e., so that sum(time>1)=x
# within each group), Testtime must be 1
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)"
#PARMTYPE<-"efflogs"
#PARMTYPE.LABEL<-"1- log(S2(t))/log(S1(t))"
Betafunc<-function(S1,S2){ S2-S1 }
transBeta<- function(beta){ beta }
#LIMIT<-"upper"
#ALT<- "less"
LIMIT<-"lower"
ALT<-"greater"
CL<- 0.975
#N1<- 12
#N2<- 24
#N1<- 18
#N2<- 36
#N1<-30
#N2<-60
# for the lines, we number the same as for the
# uncensored plots
# 1=BPCP (b_)
# 2=mid-p BPCP (bmp_)
# 3=adj hybrid Z-O (dah_)
# 4=standard Z-O (d_)
# 5=standard (no Z-O)
# Gray scale with large width changes
LWD<-c(13,12,11,6,2)
LTY<-c(1,1,1,1,1)
COL<- gray(c(.5,.7,.3,.9,.1)-.1)
# Color changes mostly
LWD<-c(5:1)*2
LTY<-c(1,1,1,2,1)
COL<-c(palette.colors(4)[c(1,3,2,4)],palette.colors(5)[5])
#plot(1:10,1:10,type="n")
#for (i in 1:5){
# lines(c(i/2,10-i/2),c(i/2,10-i/2),lwd=LWD[i],col=COL[i])
#}
CIfunc1<-function(x1,n1,x2,n2){
CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
Midp=FALSE)
}
#CIfunc2<-function(x1,n1,x2,n2){
# CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
# Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
#}
CIfunc2a<-function(x1,n1,x2,n2){
CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
Alternative=ALT, Method="standard", Zero.one.adjustment=FALSE)
}
CIfunc2b<-function(x1,n1,x2,n2){
CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
}
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,
control=bpcp2sampControl(nmc=1e6)
)
}
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)
}
set.seed(50693)
N1<-30
N2<-30
t0<- proc.time()
# meld BPCP
d3030out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
t1<- proc.time()
t1-t0
# standard
d3030out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
t2<- proc.time()
t2-t1
#standard +Z-O
d3030out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
t3<- proc.time()
t3-t2
# meld BPCP mid-p
d3030out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
t4<- proc.time()
t4-t3
# adjusted hybrid Z-O
d3030out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t5<- proc.time()
t5-t4
LIMIT
# BPCP - standard
mean(d3030out1$CLimit-d3030out2a$CLimit)
# BPCP- standard, Z-O
mean(d3030out1$CLimit-d3030out2b$CLimit)
# BPCP - BPCP mid-p
mean(d3030out1$CLimit-d3030out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(d3030out1$CLimit-d3030out4$CLimit)
t1<-proc.time()
t1-t0
#-----------------------------------------------
# First plot only the cases when S2(t)=S1(t)
# for the difference
#----------------------------------------------
pdf("./simBPCP/simResults/Calculation2026_no_censor.pdf", width = 7, height = 5)
par(mfrow=c(1,2))
N1<-30
N2<-30
# plot at 1-S(t)
# then reverse the labels so that it is
# the labels matched the plotted values at S(t)
plot(1-St,1-diag(d3030out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
# notice the labels=rev(AT)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(d3030out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(d3030out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(d3030out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(d3030out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(d3030out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
legend("topleft",legend=c("meld BPCP","meld mid-p BPCP","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
lwd=LWD[o],lty=LTY[o],col=COL[o])
t0<- proc.time()
set.seed(991321)
N1<-30
N2<-60
LIMIT<-"upper"
ALT<- "less"
d3060out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
d3060out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
d3060out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
d3060out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
d3060out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0
LIMIT
# BPCP - standard
mean(d3060out1$CLimit-d3060out2a$CLimit)
# BPCP- standard, Z-O
mean(d3060out1$CLimit-d3060out2b$CLimit)
# BPCP - BPCP mid-p
mean(d3060out1$CLimit-d3060out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(d3060out1$CLimit-d3060out4$CLimit)
N1<-30
N2<-60
plot(1-St,1-diag(d3060out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(d3060out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(d3060out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(d3060out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(d3060out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(d3060out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
#legend("topright",legend=c("BPCP","BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
# lwd=LWD[o],lty=LTY[o],col=COL[o])
dev.off()
#dev.print(pdf,file="Calculation2026_no_censor.pdf")
#############################################
# Change Betafunc to efflogs
#
#############################################
Betafunc<-function(S1,S2){
out<-1- log(S2)/log(S1)
out[S1==S2]<- 0
# because the computer incorrectly gives
# 1-log(0)/log(1) = Inf instead of -Inf
# we need to fix this manually
out[S1==1 & S2==0]<- -Inf
out
}
# check
Betafunc(0,1)
Betafunc(0,0)
Betafunc(0.99999999,0)
Betafunc(1,0)
Betafunc(1,1)
PARMTYPE<-"efflogs"
t0<- proc.time()
set.seed(2421)
N1<-30
N2<-30
e3030out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
e3030out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
e3030out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
e3030out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
e3030out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
LIMIT
# BPCP - standard
mean(e3030out1$CLimit-e3030out2a$CLimit)
# BPCP- standard, Z-O
mean(e3030out1$CLimit-e3030out2b$CLimit)
# BPCP - BPCP mid-p
mean(e3030out1$CLimit-e3030out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(e3030out1$CLimit-e3030out4$CLimit)
t1<-proc.time()
t1-t0
pdf("./simBPCP/simResults/Calculation2026_no_censor_efflogs.pdf", width = 7, height = 5)
par(mfrow=c(1,2))
N1<-30
N2<-30
plot(1-St,1-diag(e3030out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(e3030out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(e3030out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(e3030out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(e3030out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(e3030out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
legend("topleft",legend=c("meld BPCP","meld mid-p BPCP","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
lwd=LWD[o],lty=LTY[o],col=COL[o])
t0<- proc.time()
set.seed(8393)
N1<-30
N2<-60
LIMIT<-"upper"
ALT<- "less"
e3060out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
e3060out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
e3060out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
e3060out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
e3060out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0
LIMIT
# BPCP - standard
mean(e3060out1$CLimit-e3060out2a$CLimit)
# BPCP- standard, Z-O
mean(e3060out1$CLimit-e3060out2b$CLimit)
# BPCP - BPCP mid-p
mean(e3060out1$CLimit-e3060out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(e3060out1$CLimit-e3060out4$CLimit)
N1<-30
N2<-60
plot(1-St,1-diag(e3060out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(e3060out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(e3060out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(e3060out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(e3060out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(e3060out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
#legend("topright",legend=c("BPCP","BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
# lwd=LWD[o],lty=LTY[o],col=COL[o])
dev.off()
#dev.print(pdf,file="Calculation2026_no_censor_efflogs.pdf")
#nosleep_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
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' ---
#' title: "Calculation of Type I Error Rate for No Censoring Case"
#' author: "Michael Fay"
#' date: "`r Sys.Date()`"
#' output: html_document
#' ---
#'
## ----setup, include=FALSE------------------------------------------------------------------------------
DOMIDP<- TRUE
#knitr::opts_chunk$set(echo = FALSE)
#'
#'
## ------------------------------------------------------------------------------------------------------
library(exact2x2)
library(bpcp)
#library(NoSleepR)
#nosleep_on()
NGRID<- 100
St<- 0:NGRID/NGRID
#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){
# create times so that length(time)=n and sum(time>testtime)=x
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))
# because of the way we create time (i.e., so that sum(time>1)=x
# within each group), Testtime must be 1
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)"
#PARMTYPE<-"efflogs"
#PARMTYPE.LABEL<-"1- log(S2(t))/log(S1(t))"
Betafunc<-function(S1,S2){ S2-S1 }
transBeta<- function(beta){ beta }
#LIMIT<-"upper"
#ALT<- "less"
LIMIT<-"lower"
ALT<-"greater"
CL<- 0.975
#N1<- 12
#N2<- 24
#N1<- 18
#N2<- 36
#N1<-30
#N2<-60
# for the lines, we number the same as for the
# uncensored plots
# 1=BPCP (b_)
# 2=mid-p BPCP (bmp_)
# 3=adj hybrid Z-O (dah_)
# 4=standard Z-O (d_)
# 5=standard (no Z-O)
# Gray scale with large width changes
LWD<-c(13,12,11,6,2)
LTY<-c(1,1,1,1,1)
COL<- gray(c(.5,.7,.3,.9,.1)-.1)
# Color changes mostly
LWD<-c(5:1)*2
LTY<-c(1,1,1,2,1)
COL<-c(palette.colors(4)[c(1,3,2,4)],palette.colors(5)[5])
#plot(1:10,1:10,type="n")
#for (i in 1:5){
# lines(c(i/2,10-i/2),c(i/2,10-i/2),lwd=LWD[i],col=COL[i])
#}
CIfunc1<-function(x1,n1,x2,n2){
CIfuncBPCP(x1,n1,x2,n2,Conf.level=CL,Alternative=ALT,Parmtype=PARMTYPE,
Midp=FALSE)
}
#CIfunc2<-function(x1,n1,x2,n2){
# CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
# Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
#}
CIfunc2a<-function(x1,n1,x2,n2){
CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
Alternative=ALT, Method="standard", Zero.one.adjustment=FALSE)
}
CIfunc2b<-function(x1,n1,x2,n2){
CIfuncDelta(x1,n1,x2,n2,Conf.level=CL,Parmtype=PARMTYPE,
Alternative=ALT, Method="standard", Zero.one.adjustment=TRUE)
}
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,
control=bpcp2sampControl(nmc=1e6)
)
}
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)
}
set.seed(50693)
N1<-30
N2<-30
t0<- proc.time()
# meld BPCP
d3030out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
t1<- proc.time()
t1-t0
# standard
d3030out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
t2<- proc.time()
t2-t1
#standard +Z-O
d3030out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
t3<- proc.time()
t3-t2
# meld BPCP mid-p
d3030out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
t4<- proc.time()
t4-t3
# adjusted hybrid Z-O
d3030out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t5<- proc.time()
t5-t4
LIMIT
# BPCP - standard
mean(d3030out1$CLimit-d3030out2a$CLimit)
# BPCP- standard, Z-O
mean(d3030out1$CLimit-d3030out2b$CLimit)
# BPCP - BPCP mid-p
mean(d3030out1$CLimit-d3030out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(d3030out1$CLimit-d3030out4$CLimit)
t1<-proc.time()
t1-t0
#-----------------------------------------------
# First plot only the cases when S2(t)=S1(t)
# for the difference
#----------------------------------------------
pdf("./simBPCP/simResults/Calculation2026_no_censor.pdf", width = 7, height = 5)
par(mfrow=c(1,2))
N1<-30
N2<-30
# plot at 1-S(t)
# then reverse the labels so that it is
# the labels matched the plotted values at S(t)
plot(1-St,1-diag(d3030out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
# notice the labels=rev(AT)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(d3030out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(d3030out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(d3030out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(d3030out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(d3030out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
legend("topleft",legend=c("meld BPCP","meld mid-p BPCP","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
lwd=LWD[o],lty=LTY[o],col=COL[o])
t0<- proc.time()
set.seed(991321)
N1<-30
N2<-60
LIMIT<-"upper"
ALT<- "less"
d3060out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
d3060out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
d3060out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
d3060out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
d3060out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0
LIMIT
# BPCP - standard
mean(d3060out1$CLimit-d3060out2a$CLimit)
# BPCP- standard, Z-O
mean(d3060out1$CLimit-d3060out2b$CLimit)
# BPCP - BPCP mid-p
mean(d3060out1$CLimit-d3060out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(d3060out1$CLimit-d3060out4$CLimit)
N1<-30
N2<-60
plot(1-St,1-diag(d3060out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(d3060out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(d3060out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(d3060out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(d3060out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(d3060out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
#legend("topright",legend=c("BPCP","BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
# lwd=LWD[o],lty=LTY[o],col=COL[o])
dev.off()
#dev.print(pdf,file="Calculation2026_no_censor.pdf")
#############################################
# Change Betafunc to efflogs
#
#############################################
Betafunc<-function(S1,S2){
out<-1- log(S2)/log(S1)
out[S1==S2]<- 0
# because the computer incorrectly gives
# 1-log(0)/log(1) = Inf instead of -Inf
# we need to fix this manually
out[S1==1 & S2==0]<- -Inf
out
}
# check
Betafunc(0,1)
Betafunc(0,0)
Betafunc(0.99999999,0)
Betafunc(1,0)
Betafunc(1,1)
PARMTYPE<-"efflogs"
t0<- proc.time()
set.seed(2421)
N1<-30
N2<-30
e3030out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
e3030out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
e3030out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
e3030out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
e3030out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
LIMIT
# BPCP - standard
mean(e3030out1$CLimit-e3030out2a$CLimit)
# BPCP- standard, Z-O
mean(e3030out1$CLimit-e3030out2b$CLimit)
# BPCP - BPCP mid-p
mean(e3030out1$CLimit-e3030out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(e3030out1$CLimit-e3030out4$CLimit)
t1<-proc.time()
t1-t0
pdf("./simBPCP/simResults/Calculation2026_no_censor_efflogs.pdf", width = 7, height = 5)
par(mfrow=c(1,2))
N1<-30
N2<-30
plot(1-St,1-diag(e3030out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(e3030out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(e3030out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(e3030out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(e3030out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(e3030out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
legend("topleft",legend=c("meld BPCP","meld mid-p BPCP","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
lwd=LWD[o],lty=LTY[o],col=COL[o])
t0<- proc.time()
set.seed(8393)
N1<-30
N2<-60
LIMIT<-"upper"
ALT<- "less"
e3060out1<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc1,betafunc=Betafunc,limit=LIMIT)
e3060out2a<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2a,betafunc=Betafunc,limit=LIMIT)
e3060out2b<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc2b,betafunc=Betafunc,limit=LIMIT)
e3060out3<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc3,betafunc=Betafunc,limit=LIMIT)
e3060out4<-checkCoverage(n1=N1,n2=N2,ngrid=NGRID,cifunc=CIfunc4,betafunc=Betafunc,limit=LIMIT)
t1<-proc.time()
t1-t0
LIMIT
# BPCP - standard
mean(e3060out1$CLimit-e3060out2a$CLimit)
# BPCP- standard, Z-O
mean(e3060out1$CLimit-e3060out2b$CLimit)
# BPCP - BPCP mid-p
mean(e3060out1$CLimit-e3060out3$CLimit)
# BPCP - adj hybrid, Z-O
mean(e3060out1$CLimit-e3060out4$CLimit)
N1<-30
N2<-60
plot(1-St,1-diag(e3060out1$Coverage),type="n",
xlab="S(t)",ylab="Type I Error Rate",main=paste0("n1=",N1," n2=",N2),
axes=FALSE,ylim=c(0,0.08))
AT<-c(0,.2,.4,.6,.8,1)
axis(1,at=AT,labels=rev(AT))
axis(2)
box()
# 1=BPCP
lines(1-St,1-diag(e3060out1$Coverage),lwd=LWD[1],lty=LTY[1],col=COL[1])
# 3=mid-p BPCP
lines(1-St,1-diag(e3060out3$Coverage),lwd=LWD[2],lty=LTY[2],col=COL[2])
# 4=adj hybrid
lines(1-St,1-diag(e3060out4$Coverage),lwd=LWD[3],lty=LTY[3],col=COL[3])
# 2b=standard Z-O
lines(1-St,1-diag(e3060out2b$Coverage),lwd=LWD[4],lty=LTY[4],col=COL[4])
# 2a=standard
lines(1-St,1-diag(e3060out2a$Coverage),lwd=LWD[5],lty=LTY[5],col=COL[5])
lines(c(0,1),c(0.025,0.025),lty=2,col="red")
o<-c(5:1)
#legend("topright",legend=c("BPCP","BPCP mid-p","Delta (adj hybrid,Z-O)","Delta (standard,Z-O)","Delta (standard)")[o],
# lwd=LWD[o],lty=LTY[o],col=COL[o])
dev.off()
#dev.print(pdf,file="Calculation2026_no_censor_efflogs.pdf")
#nosleep_off()
Try the bpcp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.