## ----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)