## ----include = FALSE----------------------------------------------------------
 knitr::opts_chunk$set(
   collapse = TRUE,
   comment = "#>"
 )

## ----setup--------------------------------------------------------------------
library(BayesianPlatformDesignTimeTrend)

## ----eval=FALSE---------------------------------------------------------------
#  
#  ntrials = 1000 # Number of trial replicates
#  ns = seq(120, 600, 60) # Sequence of total number of accrued patients at each interim analysis
#  null.reponse.prob = 0.4
#  alt.response.prob = 0.6
#  
#  # We investigate the type I error rate for different time trend strength
#  null.scenario = matrix(
#    c(
#      null.reponse.prob,
#      null.reponse.prob,
#      null.reponse.prob,
#      null.reponse.prob
#    ),
#    nrow = 1,
#    ncol = 4,
#    byrow = T
#  )
#  alt.scenario = matrix(
#    c(
#      null.reponse.prob,
#      alt.response.prob,
#      alt.response.prob,
#      null.reponse.prob
#    ),
#    nrow = 1,
#    ncol = 4,
#    byrow = T
#  )
#  model = "tlr" #logistic model
#  max.ar = 0.85  #limit the allocation ratio for the control group (1-max.ar < r_control < max.ar)
#  #------------Select the data generation randomisation methods-------
#  rand.type = "Urn" # Urn design
#  max.deviation = 3 # The recommended value for the tuning parameter in the Urn design
#  
#  # Require multiple cores for parallel running on HPC (Here is the number of cores i ask on Iridis 5 in University of Southampton)
#  cl = 40
#  
#  # Set the model we want to use and the time trend effect for each model used.
#  # Here the main model will be used twice for two different strength of time trend c(0,0,0,0) and c(1,1,1,1) to investigate how time trend affect the evaluation metrics in BAR setting.
#  # Then the main + stage_continuous model which is the treatment effect + stage effect model will be applied for strength equal c(1,1,1,1) to investigate how the main + stage effect model improve the evaluation metrics.
#  reg.inf = "main"
#  trend.effect = c(0,0,0,0)
#  
#  result = {
#  
#  }
#  OPC = {
#  
#  }
#  K = dim(null.scenario)[2]
#  cutoffindex = 1
#  trendindex = 1
#  
#  cutoff.information=demo_Cutoffscreening.GP (
#    ntrials = ntrials,
#    # Number of trial replicates
#    trial.fun = simulatetrial,
#    # Call the main function
#    power.type = "Conjunctive",
#    response.probs.alt = alt.scenario,
#    grid.inf = list(
#      start.length = 15,
#      grid.min = NULL,
#      grid.max = NULL,
#      confidence.level = 0.95,
#      grid.length = 101,
#      change.scale = FALSE,
#      noise = T,
#      errorrate = 0.1,
#      simulationerror = 0.01,
#      iter.max = 15,
#      plotornot = FALSE),
#    # Set up the cutoff grid for screening. The start grid has three elements. The extended grid has fifteen cutoff value under investigation
#    input.info = list(
#      response.probs = null.scenario[1,],
#      #The scenario vector in this round
#      ns = ns,
#      # Sequence of total number of accrued patients at each interim analysis
#      max.ar =  max.ar,
#      #limit the allocation ratio for the control group (1-max.ar < r_control < max.ar)
#      rand.type = rand.type,
#      # Which randomisation methods in data generation.
#      max.deviation = max.deviation,
#      # The recommended value for the tuning parameter in the Urn design
#      model.inf = list(
#        model = model,
#        #Use which model?
#        ibb.inf = list(
#          #independent beta-binomial model which can be used only for no time trend simulation
#          pi.star = 0.5,
#          # beta prior mean
#          pess = 2,
#          # beta prior effective sample size
#          betabinomialmodel = ibetabinomial.post # beta-binomial model for posterior estimation
#        ),
#        tlr.inf = list(
#          beta0_prior_mu = 0,
#          # Stan logistic model t prior location
#          beta1_prior_mu = 0,
#          # Stan logistic model t prior location
#          beta0_prior_sigma = 2.5,
#          # Stan logistic model t prior sigma
#          beta1_prior_sigma = 2.5,
#          # Stan logistic model t prior sigma
#          beta0_df = 7,
#          # Stan logistic model t prior degree of freedom
#          beta1_df = 7,
#          # Stan logistic model t prior degree of freedom
#          reg.inf =  reg.inf,
#          # The model we want to use
#          variable.inf = "Fixeffect" # Use fix effect logistic model
#        )
#      ),
#      Stop.type = "Early-Pocock",
#      # Use Pocock like early stopping boundary
#      Boundary.type = "Asymmetric",
#      # Use Symmetric boundary where cutoff value for efficacy boundary and futility boundary sum up to 1
#      Random.inf = list(
#        Fixratio = FALSE,
#        # Do not use fix ratio allocation
#        Fixratiocontrol = NA,
#        # Do not use fix ratio allocation
#        BARmethod = "Thall",
#        # Use Thall's Bayesian adaptive randomisation approach
#        Thall.tuning.inf = list(tuningparameter = "Unfixed",  fixvalue = 1) # Specified the tunning parameter value for fixed tuning parameter
#      ),
#      trend.inf = list(
#        trend.type = "linear",
#        # Linear time trend pattern
#        trend.effect = trend.effect,
#        # Stength of time trend effect
#        trend_add_or_multip = "mult" # Multiplicative time trend effect on response probability
#      )
#    ),
#    cl = 2
#  )
#  

## -----------------------------------------------------------------------------
library(ggplot2)
# Details of grid
optimdata=optimdata_asy
# Recommend cutoff at each screening round
nextcutoff = optimdata$next.cutoff
nextcutoff$FWER=0.05
nextcutoff.predict = nextcutoff
colnames(nextcutoff.predict)=c("eff","fut","FWER")
prediction = optimdata$prediction
point.tested=optimdata$testeddata[,2:3]
tpIE=optimdata$testeddata[,1]
pow=optimdata$testeddata[,4]
point.tested=point.tested[1:sum(!is.na(tpIE)),]
tpIE=tpIE[1:sum(!is.na(tpIE))]
pow=pow[1:sum(!is.na(pow))]
cleandata=data.frame(FWER=tpIE,pow=pow,point.tested)
colnames(cleandata)[c(3,4)]=c("eff","fut")
GP.res = optimdata
xgrid.eff=optimdata$prediction$xgrid[,1]
xgrid.fut=optimdata$prediction$xgrid[,2]
grid.min=c(0.95,0)
grid.max=c(1,0.05)

library(grDevices)
library(RColorBrewer)
colormap=colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
target_line=0.1

df=data.frame(FWER=optimdata$prediction$yhat.t1E,eff=xgrid.eff,fut=xgrid.fut)
Contour.tIE<-ggplot(df,aes(eff,fut,z=FWER))+
  scale_fill_gradientn(colors = colormap)+geom_tile(aes(fill=FWER))+
  geom_contour(breaks=c(target_line, seq(min(df$FWER),max(df$FWER),by=(max(df$FWER)-min(df$FWER))/10)),color="black")+
  geom_contour(breaks=target_line,color="white",linewidth=1.1)+
  labs(title="Mean type I error rate (FWER)", x="Cutoff value for efficacy",y="Cutoff value for futility")
Contour.tIE=Contour.tIE+geom_point(data=cleandata,aes(eff,fut),color="black")+geom_point(data=nextcutoff.predict,aes(eff,fut),color="pink")

# Extract the contour data
contour_data_tIE <- ggplot_build(Contour.tIE)$data[[2]]
# Record the contour that has FWER equal to the target
contour_data_tIE_subset <- contour_data_tIE[contour_data_tIE$level == target_line, ]
# Order and split the data to ensure the plot is drawn correctly
contour_data_tIE_subset=contour_data_tIE_subset[order(contour_data_tIE_subset$piece,contour_data_tIE_subset$x),]
contour_data_tIE_subset_1=contour_data_tIE_subset[contour_data_tIE_subset$piece==1,]
contour_data_tIE_subset_2=contour_data_tIE_subset[contour_data_tIE_subset$piece==2,]
# To make sure the data frame is not empty
if (nrow(contour_data_tIE_subset_1) == 0){
  contour_data_tIE_subset_1[1,]=(rep(NA,dim(contour_data_tIE_subset_1)[[2]]))
} else if (nrow(contour_data_tIE_subset_2) == 0){
  contour_data_tIE_subset_2[1,]=(rep(NA,dim(contour_data_tIE_subset_2)[[2]]))
}


df=data.frame(sd=optimdata$prediction$sd.t1E,eff=xgrid.eff,fut=xgrid.fut)
Contour.sd<-ggplot(df,aes(eff,fut,z=sd))+
  scale_fill_gradientn(colors = colormap)+geom_tile(aes(fill=sd))+
  geom_contour(breaks=seq(min(df$sd),max(df$sd),by=(max(df$sd)-min(df$sd))/10),color="black")+labs(title="sd of contour", x="Cutoff value for efficacy",y="Cutoff value for futility")
Contour.sd=Contour.sd+
  geom_path(data = contour_data_tIE_subset_1, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_path(data = contour_data_tIE_subset_2, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_point(data=cleandata,aes(eff,fut,z=NA),color="black")+geom_point(data=nextcutoff.predict,aes(eff,fut,z=NA),color="pink")

df=data.frame(Power=optimdata$prediction$yhat.pow,eff=xgrid.eff,fut=xgrid.fut)
Contour.pow<-ggplot(df,aes(eff,fut,z=Power))+
  scale_fill_gradientn(colors = colormap)+geom_tile(aes(fill=Power))+
  geom_contour(breaks=seq(min(df$Power),max(df$Power),by=(max(df$Power)-min(df$Power))/10),color="black")+labs(title="Mean power", x="Cutoff value for efficacy",y="Cutoff value for futility")
Contour.pow=Contour.pow+
  geom_path(data = contour_data_tIE_subset_1, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_path(data = contour_data_tIE_subset_2, aes(x,y,z=NA),color="white",linewidth=1.1)+
  geom_point(data=cleandata,aes(eff,fut,z=NA),color="black")+geom_point(data=nextcutoff.predict,aes(eff,fut,z=NA),color="pink")

df=data.frame(NullESS=optimdata$prediction$yhat.ESS.null,eff=xgrid.eff,fut=xgrid.fut)
Contour.nullESS<-ggplot(df,aes(eff,fut,z=NullESS))+
  scale_fill_gradientn(colors = colormap)+geom_tile(aes(fill=NullESS))+
  geom_contour(breaks=seq(min(df$NullESS),max(df$NullESS),by=(max(df$NullESS)-min(df$NullESS))/10),color="black")+labs(title="Mean ESS under null", x="Cutoff value for efficacy",y="Cutoff value for futility")
Contour.nullESS=Contour.nullESS+
  geom_path(data = contour_data_tIE_subset_1, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_path(data = contour_data_tIE_subset_2, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_point(data=cleandata,aes(eff,fut,z=NA),color="black")+geom_point(data=nextcutoff.predict,aes(eff,fut,z=NA),color="pink")

df=data.frame(AltESS=optimdata$prediction$yhat.ESS.alt,eff=xgrid.eff,fut=xgrid.fut)
Contour.altESS<-ggplot(df,aes(eff,fut,z=AltESS))+
  scale_fill_gradientn(colors = colormap)+geom_tile(aes(fill=AltESS))+
  geom_contour(breaks=seq(min(df$AltESS),max(df$AltESS),by=(max(df$AltESS)-min(df$AltESS))/10),color="black")+labs(title="Mean ESS under alternative", x="Cutoff value for efficacy",y="Cutoff value for futility")
Contour.altESS=Contour.altESS+
  geom_path(data = contour_data_tIE_subset_1, aes(x,y,z=NA),color="white",linewidth=1.1)+geom_path(data = contour_data_tIE_subset_2, aes(x,y,z=NA),color="white",linewidth=1.1)+
  geom_point(data=cleandata,aes(eff,fut,z=NA),color="black")+geom_point(data=nextcutoff.predict,aes(eff,fut,z=NA),color="pink")

## ----fig.align='center',fig.height=9,fig.width=7,warning=FALSE----------------
library(ggpubr)
ggarrange(Contour.tIE,Contour.pow,Contour.nullESS,Contour.altESS,Contour.sd,ncol = 2,nrow=3)