Nothing
inst/doc/DTEBOP2.R
In DTEBOP2: Bayesian Optimal Phase II Randomized Clinical Trial Design with Delayed Outcomes
## ----include = FALSE----------------------------------------------------------
#cat("<style>
#pre {
# color: black !important;
#}
#</style>")
# 全局代码块选项设置
knitr::opts_chunk$set(
collapse = TRUE, # 合并代码和输出
comment = "#>", # 输出前的注释符号
fig.width = 8, # 默认图像宽度（英寸）
fig.height = 6, # 默认图像高度（英寸）
dpi = 300, # 图像分辨率（DPI）
out.width = "100%" # 图像在 HTML 中的宽度占比
)
## ----echo=FALSE---------------------------------------------------------------
knitr::opts_chunk$set(
fig.width = 10, # 默认图像宽度
fig.height = 7, # 默认图像高度
dpi = 300, # 图像分辨率
out.width = "80%" # 图像在页面中显示为 100% 宽度
)
## ----setup,message=FALSE,warning=FALSE----------------------------------------
library(DTEBOP2)
library(invgamma)
library(truncdist)
library(doParallel)
library(survRM2adapt)
library(reconstructKM)
library(ggplot2)
data("Experimental")
data("SOC")
Experiment_NAR <- data.frame(time=seq(from=0, to=24, by=3), NAR=c(135, 68, 48, 33, 21, 15, 6, 2,0))
SOC_NAR <- data.frame(time=seq(from=0, to=24, by=3), NAR=c(137, 62, 26, 9, 6, 2,1, 0, 0))
Experiment_aug <- format_raw_tabs(raw_NAR=Experiment_NAR,
raw_surv=Experimental)
SOC_aug <- format_raw_tabs(raw_NAR=SOC_NAR,
raw_surv=SOC)
# reconstruct by calling KM_reconstruct()
Experiment_recon <- KM_reconstruct(aug_NAR=Experiment_aug$aug_NAR, aug_surv=Experiment_aug$aug_surv)
SOC_recon <- KM_reconstruct(aug_NAR=SOC_aug$aug_NAR, aug_surv=SOC_aug$aug_surv)
# put the treatment and control arms into one dataset
Experiment_IPD <- data.frame(arm=1, time=Experiment_recon$IPD_time, status=Experiment_recon$IPD_event)
SOC_IPD <- data.frame(arm=0, time=SOC_recon$IPD_time, status=SOC_recon$IPD_event)
allIPD <- rbind(Experiment_IPD, SOC_IPD)
KM_fit <- survival::survfit(survival::Surv(time, status) ~ arm, data=allIPD)
KM_plot <- survminer::ggsurvplot(
fit = KM_fit,
data = allIPD,
risk.table = TRUE,
palette = c('red', 'blue'),
legend = c(0.8, 0.9),
legend.title = '',
legend.labs = c('Docetaxel', 'Nivolumab'),
title = 'PFS',
ylab = 'Survival (%)',
xlab = 'Time (Month)',
risk.table.title = 'Number at Risk',
break.time.by = 3,
tables.height = 0.4,
risk.table.fontsize = 5
)
print(KM_plot,fig.align='center')
## ----eval=FALSE---------------------------------------------------------------
# expert_data_correct <- list(
# list(mean = 2.2, median = 2.27, sd = NULL, q25 = NULL, q975 = NULL), # expert A
# list(mean = 2.1, median = 2.3, sd = NULL, q25 = NULL, q975 = NULL), # expert B
# list(mean = 2.3, median = 2.31, sd = NULL, q25 = NULL, q975 = NULL) # expert C
# )
# param<-trunc_gamma_para(L = 2, U = 2.5, expert_data = expert_data_correct,weights = c(4,4,2,1,1) ,num_cores = 4,seed=123) # num_cores specifies the number of core to do the parallel computing.
# print(param)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=12.86,scale=0.19)),xlab="",main="")
## -----------------------------------------------------------------------------
y_prob <- c(0.54, 0.53, 0.52,0.51) #y-axis survival probability to explore
summary_fit <- summary(KM_fit)
arm_times <- list()
for (i in 1:length(KM_fit$strata)) {
time_points <- summary_fit$time[summary_fit$strata == names(KM_fit$strata)[i]]
survival_probs <- summary_fit$surv[summary_fit$strata == names(KM_fit$strata)[i]]
closest_times <- sapply(y_prob, function(tp) {
closest_index <- which.min(abs(y_prob - tp))
return(time_points[closest_index])
})
arm_times[[names(KM_fit$strata)[i]]] <- setNames(closest_times, y_prob)
}
print(arm_times)
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2.28 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior paramters of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last enrolled patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
#
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, weight=0.5,earlystop_prob =NULL,seed=123)
# #calculate the sample size for two-stage design
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_sample_size=readRDS("sample_size_1.rds")
result_sample_size=c(28,40)
print(result_sample_size)
## -----------------------------------------------------------------------------
event_fun(median.1 = 2.8,median.2 = 3.5,S_likely = 2.28,n.interim = c(28,40),rate = 6,FUP = 6,n.sim = 1000,seed=123)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2 # The most likely delayed time
# trunc.para=c(12.86,0.19) # default prior parameter of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, weight=0.5,earlystop_prob =NULL,seed=123) #calculate the sample size for two-stage design
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_sample_size_2=readRDS("sample_size_2.rds")
result_sample_size_2=c(44,63)
print(result_sample_size_2)
## -----------------------------------------------------------------------------
event_fun(median.1 = 2.8,median.2 = 3.5,S_likely = 2,n.interim = c(44,63),rate = 6,FUP = 6,n.sim = 1000,seed=123)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
set.seed(126)
n.interim = c(44,63)
nmax=63
FUP=6
rate=6
L=2
U=2.5
median.1=2.8
median.2=3.5
lambda_1=log(2)/2.8 #hazard rate before S
lambda_2=log(2)/2.8 #hazard rate after S
nobs = n.interim+1
nobs[length(nobs)] = 63
wait.t = rexp(nmax,rate = rate)
arrival.t = cumsum(wait.t)
tobs = arrival.t[nobs]
tobs[length(tobs)] = tobs[length(tobs)] + FUP
S_likely=2.1
event.t.E = generate_pe(nmax,t=S_likely,lambda1=lambda_1,lambda2=lambda_1)
event.t.C=rexp(nmax,rate=lambda_1)
####first interim
k = 1
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E) ## treatment dataset
data.C=data.frame(Time=t.event.C,Status=Ind_event_C) ## control dataset
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=63,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
k = 2
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E)
data.C=data.frame(Time=t.event.C,Status=Ind_event_C)
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=nmax,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
set.seed(126)
n.interim = c(44,63)
nmax=63
FUP=6
rate=6
L=2
U=2.5
median.1=2.8
median.2=3.5
lambda_1=log(2)/2.8 #hazard rate before S
lambda_2=log(2)/6.6 #hazard rate after S
nobs = n.interim+1
nobs[length(nobs)] = 63
wait.t = rexp(nmax,rate = rate)
arrival.t = cumsum(wait.t)
tobs = arrival.t[nobs]
tobs[length(tobs)] = tobs[length(tobs)] + FUP
S_likely=2.1
event.t.E = generate_pe(nmax,t=S_likely,lambda1=lambda_1,lambda2=lambda_2)
event.t.C=rexp(nmax,rate=lambda_1)
####first interim
k = 1
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E) ## treatment dataset
data.C=data.frame(Time=t.event.C,Status=Ind_event_C) ## control dataset
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=63,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
k = 2
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E)
data.C=data.frame(Time=t.event.C,Status=Ind_event_C)
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=nmax,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2.28 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, earlystop_prob =0.4,seed=123)
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#earlystop_sample_size_1=readRDS("early_stopping_sample_size_1.rds")
earlystop_sample_size_1=c(39,52)
print(earlystop_sample_size_1)
## -----------------------------------------------------------------------------
##S=2
n.interim=c(39,52)
getoc_2arm_piecewise(median.true=c(median.1,median.1),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = n.interim,S_likely=2.28,L=2,U=2,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)$earlystop
##S=2.5
getoc_2arm_piecewise(median.true=c(median.1,median.1),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = n.interim,S_likely=2.28,L=2.5,U=2.5,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)$earlystop
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, earlystop_prob =0.4,seed=123)
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#earlystop_sample_size_2=readRDS("early_stopping_sample_size_2.rds")
earlystop_sample_size_2=c(46,65)
print(earlystop_sample_size_2)
## ----eval=FALSE---------------------------------------------------------------
# n.interim=c(13,28,40)
# get.optimal_2arm_piecewise(median.1=median.1,median.2=median.2, L=L,U=U,S_likely=2.28,Uniform=FALSE,trunc.para=c(1,1), err1=0.1, n.interim=n.interim, rate=rate, FUP=FUP,track=TRUE,nsim = 10000,control = T,control.point = c(L,U)) #Get the optimal lambda and gamma
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal_4=readRDS("result_optimal_4.rds")
result_optimal_4=c(0.9500, 1.0000, 0.0875, 0.8678)
print(result_optimal_4)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
## ----eval=FALSE---------------------------------------------------------------
# expert_data_correct <- list(
# list(mean = 2.28, median = NULL, sd = NULL, q25 = NULL, q975 = NULL), # Expert A
# list(mean = NULL, median = 2.4, sd = NULL, q25 = NULL, q975 = NULL), # Expert B
# list(mean = 2.4, median = 2.3, sd = NULL, q25 = NULL, q975 = NULL) # Expert C
# )
#
# param_1 <- trunc_gamma_para(
# L = 2,
# U = 2.5,
# expert_data = expert_data_correct,
# weights = c(4,4,2,1,1),
# num_cores = 4,seed=123 # Number of cores for parallel computing
# )
#
# print(param_1)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=19.49,scale=0.27)),xlab="",main="")
## ----eval=FALSE---------------------------------------------------------------
# n.interim = c(28,40)
# trunc.para=param_1
# get.optimal_2arm_piecewise(median.1=median.1,median.2=median.2, L=L,U=U,S_likely=S_likely,Uniform=FALSE,trunc.para=trunc.para, err1=err1, n.interim=n.interim, rate=rate, FUP=FUP,track=TRUE,nsim = 10000,control = T,control.point = c(L,U)) #Get the optimal lambda and gamma
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal_3=readRDS("result_optimal_3.rds")
result_optimal_3=c(0.9500, 1.0000, 0.0870, 0.8628)
print(result_optimal_3)
## -----------------------------------------------------------------------------
median.1=2.8 #\bar{mu}_0
median.2=3.5 #\bar{mu}_1
L=2 # Lower bound of S
U=2.5 # Upper bound of S
S_likely=2.28 # The most likely delayed time
trunc.para=c(19.49,0.27) # prior of S
rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
FUP=6 # Follow-up time of the last patient
err1 = 0.1 # type I error
err2=0.15
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=S_likely,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=S_likely,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=1,scale=1)),xlab="",main="")
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal=readRDS("result_optimal.rds")
result_optimal=c(0.9500, 1.0000, 0.0876, 0.8694)
print(result_optimal)
Try the DTEBOP2 package in your browser
Any scripts or data that you put into this service are public.
DTEBOP2 documentation built on June 8, 2025, 1:24 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.
## ----include = FALSE----------------------------------------------------------
#cat("<style>
#pre {
# color: black !important;
#}
#</style>")
# 全局代码块选项设置
knitr::opts_chunk$set(
collapse = TRUE, # 合并代码和输出
comment = "#>", # 输出前的注释符号
fig.width = 8, # 默认图像宽度（英寸）
fig.height = 6, # 默认图像高度（英寸）
dpi = 300, # 图像分辨率（DPI）
out.width = "100%" # 图像在 HTML 中的宽度占比
)
## ----echo=FALSE---------------------------------------------------------------
knitr::opts_chunk$set(
fig.width = 10, # 默认图像宽度
fig.height = 7, # 默认图像高度
dpi = 300, # 图像分辨率
out.width = "80%" # 图像在页面中显示为 100% 宽度
)
## ----setup,message=FALSE,warning=FALSE----------------------------------------
library(DTEBOP2)
library(invgamma)
library(truncdist)
library(doParallel)
library(survRM2adapt)
library(reconstructKM)
library(ggplot2)
data("Experimental")
data("SOC")
Experiment_NAR <- data.frame(time=seq(from=0, to=24, by=3), NAR=c(135, 68, 48, 33, 21, 15, 6, 2,0))
SOC_NAR <- data.frame(time=seq(from=0, to=24, by=3), NAR=c(137, 62, 26, 9, 6, 2,1, 0, 0))
Experiment_aug <- format_raw_tabs(raw_NAR=Experiment_NAR,
raw_surv=Experimental)
SOC_aug <- format_raw_tabs(raw_NAR=SOC_NAR,
raw_surv=SOC)
# reconstruct by calling KM_reconstruct()
Experiment_recon <- KM_reconstruct(aug_NAR=Experiment_aug$aug_NAR, aug_surv=Experiment_aug$aug_surv)
SOC_recon <- KM_reconstruct(aug_NAR=SOC_aug$aug_NAR, aug_surv=SOC_aug$aug_surv)
# put the treatment and control arms into one dataset
Experiment_IPD <- data.frame(arm=1, time=Experiment_recon$IPD_time, status=Experiment_recon$IPD_event)
SOC_IPD <- data.frame(arm=0, time=SOC_recon$IPD_time, status=SOC_recon$IPD_event)
allIPD <- rbind(Experiment_IPD, SOC_IPD)
KM_fit <- survival::survfit(survival::Surv(time, status) ~ arm, data=allIPD)
KM_plot <- survminer::ggsurvplot(
fit = KM_fit,
data = allIPD,
risk.table = TRUE,
palette = c('red', 'blue'),
legend = c(0.8, 0.9),
legend.title = '',
legend.labs = c('Docetaxel', 'Nivolumab'),
title = 'PFS',
ylab = 'Survival (%)',
xlab = 'Time (Month)',
risk.table.title = 'Number at Risk',
break.time.by = 3,
tables.height = 0.4,
risk.table.fontsize = 5
)
print(KM_plot,fig.align='center')
## ----eval=FALSE---------------------------------------------------------------
# expert_data_correct <- list(
# list(mean = 2.2, median = 2.27, sd = NULL, q25 = NULL, q975 = NULL), # expert A
# list(mean = 2.1, median = 2.3, sd = NULL, q25 = NULL, q975 = NULL), # expert B
# list(mean = 2.3, median = 2.31, sd = NULL, q25 = NULL, q975 = NULL) # expert C
# )
# param<-trunc_gamma_para(L = 2, U = 2.5, expert_data = expert_data_correct,weights = c(4,4,2,1,1) ,num_cores = 4,seed=123) # num_cores specifies the number of core to do the parallel computing.
# print(param)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=12.86,scale=0.19)),xlab="",main="")
## -----------------------------------------------------------------------------
y_prob <- c(0.54, 0.53, 0.52,0.51) #y-axis survival probability to explore
summary_fit <- summary(KM_fit)
arm_times <- list()
for (i in 1:length(KM_fit$strata)) {
time_points <- summary_fit$time[summary_fit$strata == names(KM_fit$strata)[i]]
survival_probs <- summary_fit$surv[summary_fit$strata == names(KM_fit$strata)[i]]
closest_times <- sapply(y_prob, function(tp) {
closest_index <- which.min(abs(y_prob - tp))
return(time_points[closest_index])
})
arm_times[[names(KM_fit$strata)[i]]] <- setNames(closest_times, y_prob)
}
print(arm_times)
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2.28 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior paramters of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last enrolled patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
#
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, weight=0.5,earlystop_prob =NULL,seed=123)
# #calculate the sample size for two-stage design
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_sample_size=readRDS("sample_size_1.rds")
result_sample_size=c(28,40)
print(result_sample_size)
## -----------------------------------------------------------------------------
event_fun(median.1 = 2.8,median.2 = 3.5,S_likely = 2.28,n.interim = c(28,40),rate = 6,FUP = 6,n.sim = 1000,seed=123)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2 # The most likely delayed time
# trunc.para=c(12.86,0.19) # default prior parameter of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, weight=0.5,earlystop_prob =NULL,seed=123) #calculate the sample size for two-stage design
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_sample_size_2=readRDS("sample_size_2.rds")
result_sample_size_2=c(44,63)
print(result_sample_size_2)
## -----------------------------------------------------------------------------
event_fun(median.1 = 2.8,median.2 = 3.5,S_likely = 2,n.interim = c(44,63),rate = 6,FUP = 6,n.sim = 1000,seed=123)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the average type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the average power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2,U=2,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(44,63),L=2.5,U=2.5,S_likely=2.5,FUP=6,trunc.para=c(12.86,0.19),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
set.seed(126)
n.interim = c(44,63)
nmax=63
FUP=6
rate=6
L=2
U=2.5
median.1=2.8
median.2=3.5
lambda_1=log(2)/2.8 #hazard rate before S
lambda_2=log(2)/2.8 #hazard rate after S
nobs = n.interim+1
nobs[length(nobs)] = 63
wait.t = rexp(nmax,rate = rate)
arrival.t = cumsum(wait.t)
tobs = arrival.t[nobs]
tobs[length(tobs)] = tobs[length(tobs)] + FUP
S_likely=2.1
event.t.E = generate_pe(nmax,t=S_likely,lambda1=lambda_1,lambda2=lambda_1)
event.t.C=rexp(nmax,rate=lambda_1)
####first interim
k = 1
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E) ## treatment dataset
data.C=data.frame(Time=t.event.C,Status=Ind_event_C) ## control dataset
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=63,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
k = 2
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E)
data.C=data.frame(Time=t.event.C,Status=Ind_event_C)
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=nmax,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
set.seed(126)
n.interim = c(44,63)
nmax=63
FUP=6
rate=6
L=2
U=2.5
median.1=2.8
median.2=3.5
lambda_1=log(2)/2.8 #hazard rate before S
lambda_2=log(2)/6.6 #hazard rate after S
nobs = n.interim+1
nobs[length(nobs)] = 63
wait.t = rexp(nmax,rate = rate)
arrival.t = cumsum(wait.t)
tobs = arrival.t[nobs]
tobs[length(tobs)] = tobs[length(tobs)] + FUP
S_likely=2.1
event.t.E = generate_pe(nmax,t=S_likely,lambda1=lambda_1,lambda2=lambda_2)
event.t.C=rexp(nmax,rate=lambda_1)
####first interim
k = 1
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E) ## treatment dataset
data.C=data.frame(Time=t.event.C,Status=Ind_event_C) ## control dataset
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=63,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## -----------------------------------------------------------------------------
k = 2
Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0) ##event indicator for th treatment, 1: event.
Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
t.event.E = rep(0,n.interim[k])
t.event.C=rep(0,n.interim[k])
for(l in 1:length(t.event.E)){
t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])
}
data.E=data.frame(Time=t.event.E,Status=Ind_event_E)
data.C=data.frame(Time=t.event.C,Status=Ind_event_C)
print("Treatment Dataset:")
print(head(data.E))
print("Control Dataset:")
print(head(data.C))
conduct(data.E,data.C,median.1=median.1,median.2=median.2,lambda=0.95,gamma=1,nmax=nmax,S_likely=S_likely,L=L,U=U,trunc.para=c(12.86,0.19))[[4]]
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2.28 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, earlystop_prob =0.4,seed=123)
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#earlystop_sample_size_1=readRDS("early_stopping_sample_size_1.rds")
earlystop_sample_size_1=c(39,52)
print(earlystop_sample_size_1)
## -----------------------------------------------------------------------------
##S=2
n.interim=c(39,52)
getoc_2arm_piecewise(median.true=c(median.1,median.1),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = n.interim,S_likely=2.28,L=2,U=2,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)$earlystop
##S=2.5
getoc_2arm_piecewise(median.true=c(median.1,median.1),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = n.interim,S_likely=2.28,L=2.5,U=2.5,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)$earlystop
## ----eval=FALSE---------------------------------------------------------------
# median.1=2.8 #\bar{mu}_0
# median.2=3.5 #\bar{mu}_1
# L=2 # Lower bound of S
# U=2.5 # Upper bound of S
# S_likely=2 # The most likely delayed time
# trunc.para=c(12.86,0.19) # prior of S
# rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
# FUP=6 # Follow-up time of the last patient
# err1 = 0.1 # type I error
# err2=0.15 # type II error
# Two_stage_sample_size(median.1=median.1 , median.2=median.2,L=L,U=U,S_likely=S_likely,err1=err1,err2=err2,trunc.para=trunc.para, FUP=FUP, rate=rate, earlystop_prob =0.4,seed=123)
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#earlystop_sample_size_2=readRDS("early_stopping_sample_size_2.rds")
earlystop_sample_size_2=c(46,65)
print(earlystop_sample_size_2)
## ----eval=FALSE---------------------------------------------------------------
# n.interim=c(13,28,40)
# get.optimal_2arm_piecewise(median.1=median.1,median.2=median.2, L=L,U=U,S_likely=2.28,Uniform=FALSE,trunc.para=c(1,1), err1=0.1, n.interim=n.interim, rate=rate, FUP=FUP,track=TRUE,nsim = 10000,control = T,control.point = c(L,U)) #Get the optimal lambda and gamma
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal_4=readRDS("result_optimal_4.rds")
result_optimal_4=c(0.9500, 1.0000, 0.0875, 0.8678)
print(result_optimal_4)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2,U=2,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(13,28,40),L=2.5,U=2.5,S_likely=2.28,FUP=6,trunc.para=c(1,1),rate=6,nsim=10000,seed=10)
## ----eval=FALSE---------------------------------------------------------------
# expert_data_correct <- list(
# list(mean = 2.28, median = NULL, sd = NULL, q25 = NULL, q975 = NULL), # Expert A
# list(mean = NULL, median = 2.4, sd = NULL, q25 = NULL, q975 = NULL), # Expert B
# list(mean = 2.4, median = 2.3, sd = NULL, q25 = NULL, q975 = NULL) # Expert C
# )
#
# param_1 <- trunc_gamma_para(
# L = 2,
# U = 2.5,
# expert_data = expert_data_correct,
# weights = c(4,4,2,1,1),
# num_cores = 4,seed=123 # Number of cores for parallel computing
# )
#
# print(param_1)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=19.49,scale=0.27)),xlab="",main="")
## ----eval=FALSE---------------------------------------------------------------
# n.interim = c(28,40)
# trunc.para=param_1
# get.optimal_2arm_piecewise(median.1=median.1,median.2=median.2, L=L,U=U,S_likely=S_likely,Uniform=FALSE,trunc.para=trunc.para, err1=err1, n.interim=n.interim, rate=rate, FUP=FUP,track=TRUE,nsim = 10000,control = T,control.point = c(L,U)) #Get the optimal lambda and gamma
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal_3=readRDS("result_optimal_3.rds")
result_optimal_3=c(0.9500, 1.0000, 0.0870, 0.8628)
print(result_optimal_3)
## -----------------------------------------------------------------------------
median.1=2.8 #\bar{mu}_0
median.2=3.5 #\bar{mu}_1
L=2 # Lower bound of S
U=2.5 # Upper bound of S
S_likely=2.28 # The most likely delayed time
trunc.para=c(19.49,0.27) # prior of S
rate=6 # accrual rate. It assumes that the arrival time follow Poisson(rate)
FUP=6 # Follow-up time of the last patient
err1 = 0.1 # type I error
err2=0.15
##When H0 holds, the reject rate is the type I error.
getoc_2arm_piecewise(median.true=c(2.8,2.8),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=S_likely,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)
##When H1 holds, the reject rate is the power.
getoc_2arm_piecewise(median.true=c(2.8,3.5),Uniform=FALSE,lambda=0.95,gamma=1,n.interim = c(28,40),L=2.5,U=2.5,S_likely=S_likely,FUP=FUP,trunc.para=trunc.para,rate=rate,nsim=10000,seed=10)
## -----------------------------------------------------------------------------
plot(density(rtrunc_gamma(1000,L=2,U=2.5,shape=1,scale=1)),xlab="",main="")
## ----echo=FALSE---------------------------------------------------------------
#setwd("C:\\Users\\zcai13\\OneDrive - St. Jude Children's Research Hospital\\Desktop\\meeing_ht\\meeing_ht\\Jianrong Wu\\NPH\\BOP2+NPH\\R code\\cache")
#result_optimal=readRDS("result_optimal.rds")
result_optimal=c(0.9500, 1.0000, 0.0876, 0.8694)
print(result_optimal)
Try the DTEBOP2 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.