MODIFIED FOR CRAN

Author: ZXW834

.DS_Store

@@ -1,25 +1,25 @@
 Package: BayesianPlatformDesignTimeTrend
-Title: Simulate and analyse Bayesian Platform Trial with time trend
-Version: 1.0.0
+Title: Simulate and Analyse Bayesian Platform Trial with Time Trend
+Version: 1.1.2
 Author: Ziyan Wang [aut, cre]
 Maintainer: Ziyan Wang <zw7g20@soton.ac.uk>
 Authors@R: 
 	person(given = "Ziyan",
            family = "Wang",
            role = c("aut", "cre"),
            email = "zw7g20@soton.ac.uk")
-Description: Simulating the multi-arm multi-stage or platform trial with Bayesian approach using the 'rstan' package, which provides the R interface for to the Stan. 
-    The randomisation approach in this package are fix ratio and Bayesian adaptive randomisation. 
-    In addition, the time trend problem of platform trial can be studied in this package.
-    There is a demo for multi-arm multi-stage trial with two different null scenario and a demo for Bayesian trial cutoff screening.
+Description: Simulating the sequential multi-arm multi-stage or platform trial with Bayesian approach using the 'rstan' package, which provides the R interface for the Stan. 
+    This package supports fixed ratio and Bayesian adaptive randomization approaches for randomization. 
+    Additionally, it allows for the study of time trend problems in platform trials. 
+    There are demos available for a multi-arm multi-stage trial with two different null scenarios, as well as for Bayesian trial cutoff screening.
     The Bayesian adaptive randomisation approaches are described in:
     Trippa et al. (2012) <doi:10.1200/JCO.2011.39.8420> and
-    Wathen et al. (2017) <doi:10.1177/1740774517692302>
-    The randomisation algorithm is described in 
-    Zhao W <doi:10.1016/j.cct.2015.06.008>
+    Wathen et al. (2017) <doi:10.1177/1740774517692302>.
+    The randomisation algorithm is described in: 
+    Zhao W <doi:10.1016/j.cct.2015.06.008>.
     The analysis methods of time trend effect in platform trial are described in:
     Saville et al. (2022) <doi:10.1177/17407745221112013> and
-    Bofill Roig et al. (2022) <doi:10.1186/s12874-022-01683-w>
+    Bofill Roig et al. (2022) <doi:10.1186/s12874-022-01683-w>.
 URL: https://github.com/ZXW834/PlatFormDesignTime
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
@@ -42,7 +42,8 @@ Imports:
     iterators (>= 1.0.13),
     reshape (>= 0.8.8),
     BiocManager (>= 1.30.19),
-    rjags
+    lhs (>= 1.1.6),
+    laGP (>= 1.5-9)
 LinkingTo: 
     BH (>= 1.66.0),
     Rcpp (>= 0.12.0),

.Rhistory

@@ -1,2 +1,2 @@
 YEAR: 2023
-COPYRIGHT HOLDER: PlatFormDesign authors
+COPYRIGHT HOLDER: BayesianPlatformDesignTimeTrend authors

DESCRIPTION

@@ -4,10 +4,9 @@ export(ARmethod)
 export(AdaptiveRandomisation)
 export(Boundaryconstruction)
 export(FWER_disconjunctivepowerfunc)
+export(GP.optim)
 export(Initializetrialparameter)
 export(Meanfunc)
-export(Mixeffect_analysis)
-export(Mixeffect_modelling)
 export(Nfunc)
 export(OutputStats.initialising)
 export(Randomisation.inf)
@@ -18,6 +17,7 @@ export(Timetrend.fun)
 export(Trial.simulation)
 export(alphaspending)
 export(demo_Cutoffscreening)
+export(demo_Cutoffscreening.GP)
 export(demo_multscenario)
 export(ibetabinomial.post)
 export(intbias)
@@ -46,17 +46,19 @@ importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
 importFrom(graphics,lines)
 importFrom(iterators,icount)
+importFrom(laGP,distance)
 importFrom(matrixStats,colVars)
 importFrom(rstan,rstan_options)
 importFrom(rstan,sampling)
 importFrom(stats,dbeta)
 importFrom(stats,integrate)
 importFrom(stats,lm)
 importFrom(stats,model.matrix)
+importFrom(stats,optimize)
 importFrom(stats,pbeta)
 importFrom(stats,pnorm)
 importFrom(stats,predict)
+importFrom(stats,qnorm)
 importFrom(stats,rbinom)
-importFrom(stats,update)
 importFrom(stats,var)
 useDynLib(BayesianPlatformDesignTimeTrend, .registration = TRUE)

LICENSE

@@ -1,2 +1,76 @@
-# BayesianPlatformDesignTimeTrend 1.0.0
-* Added a `NEWS.md` file to track changes to the package.
+# Jun 25, 2023, BayesianPlatformDesignTimeTrend version 1.1.2
+This release fix two bugs.
+
+- I fixed the bug in function `GP.optim` where the formula of information weighed randomisation is wrong.
+- I fixed the bug in function `demo_Cutoffscreening` where the nextcutoff vector for sample may have only one element. This will lead to the error in function `sample` when you only want to sample one value greater than 1. The argument 'ntrials' in each example should be large (> 100) instead 2 to make the example more like an actual simulation example. I use ntrials = 2 in the example to speed up the check process.
+
+# Jun 11, 2023, BayesianPlatformDesignTimeTrend version 1.1.1
+This release fix one bug reported by the cran team.
+
+- I fixed the bug in function `GP.optim` where the nextcutoff vector for sample may have only one element. This will lead to the error in function `sample` when you only want to sample one value greater than 1.
+- The argument 'ntrials' in each example should be large (> 100) instead 2 to make the example more like an actual simulation example. I use ntrials = 2 in the example to speed up the check process.
+
+# May 2, 2023, BayesianPlatformDesignTimeTrend version 1.1.0
+This release implements additional method of cutoff screening using Bayesian optimization.
+
+## New feature
+- Add a new Demo function called `demo_Cutoffscreening_GP`. The function indicates how to use Gaussian process-based Bayesian optimisation algorithm to calibrate the critical value for stopping boundary under the null scenario to control the Type I error rate or FWER before trial simulation. 
+
+- Add a function `GP_optim` which returns a cutoff value for the next evaluation during the cutoff screening process. The function is used in demo `demo_Cutoffscreening_GP`.
+
+## Bugs
+- I found that the output data matrix can not be generated  when there are four arm, the 'variable.inf' = "Mixeffect" and 'reg.inf' is "main"/"main + stage_continuous"/"main * stage_discrete".
+It only work for 'variable.inf'="Mixeffect" and 'reg.inf = main + stage_discrete'. The reason why this happen is that the length of stats6 and stats7 variable are not generated correctly at stage 1 since the fixed effect model is used for 'variable.inf'="Mixeffect" at stage 1 to speed up the simulation.
+
+## Debug
+- I add one ifelse command to enforce the generation of stats6 and stats7 when 'variable.inf'="Mixeffect". It does not depend on the value of 'reg.inf' now.
+
+# May 1, 2023, BayesianPlatformDesignTimeTrend version 1.0.1
+## Major changes
+- I fixed a bug where the fix ratio allocation method code is not suitable for the multi-arm trial
+
+- Old code: rpk[-1] = rep((armleft - 1) / (armleft - 1 + Fixratiocontrol), armleft - 1)
+- New code: rpk[-1] = rep(1 / (armleft - 1 + Fixratiocontrol), armleft - 1)
+
+# Apr 25, 2023, BayesianPlatformDesignTimeTrend version 1.0.1
+## Major changes
+- I added a check command for max.ar which is the upper boundary for randomisation ratio for each arm. The command is added in the MainFunction.R. Details are shown below.
+-   #-max.ar check
+  if (1 - max.ar > 1/K){
+    stop("Error: The lower allocation ratio should be at least 1/K. Please check the number of arm at the beginning and the max.ar")
+  }
+  
+- Another change is the randomisation ratio adjustment in the Simulation_AdaptiveRandomisationmethodRatioCalc.R.
+
+- I modified the randomisation ratio adjustment command for Thall's approach.
+
+- ##----Original code------ (This code only protects the control arm's allocation ratio)
+- rpk = matrix(rep(0, armleft), ncol = armleft)
+- randomprob = matrix(rep(0, K), ncol = K)
+- colnames(rpk) = c(1, treatmentindex + 1)
+- colnames(randomprob) = seq(1, K)
+- rpk[1] = alloc.prob.best[1]
+- rpk[-1] = alloc.prob.best[-1][treatmentindex]
+- rpk = rpk / sum(rpk)
+- lower = ifelse(rpk[1] < (1 - max.ar), 1 - max.ar, rpk[1])
+- rpk[1] = lower
+- upper = ifelse(rpk[1] > max.ar, max.ar, rpk[1])
+- rpk[1] = upper
+- rpk[-1] = (1 - rpk[1]) * rpk[-1] / sum(rpk[-1])
+- randomprob[as.numeric(colnames(rpk))] = randomprob[as.numeric(colnames(rpk))] + rpk
+
+- ##----New code------ (This code protects all arms' allocation ratio. Also, this code work together with max.ar check code)
+- rpk = matrix(rep(0,armleft),ncol = armleft)
+- randomprob = matrix(rep(0,K),ncol = K)
+- colnames(rpk) = c(1,treatmentindex+1)
+- colnames(randomprob) = seq(1,K)
+- rpk[1] = alloc.prob.best[1]
+- rpk[-1] = alloc.prob.best[-1][treatmentindex]
+- rpk = rpk/sum(rpk)
+- lower = ifelse(rpk<(1-max.ar),1-max.ar,rpk)
+- rpk = lower
+- upper = ifelse(rpk>max.ar,max.ar,rpk)
+- rpk = upper
+- rpk[!(rpk==(1-max.ar))]=(1-sum(rpk[rpk==1-max.ar]))*(rpk[!(rpk==1-max.ar)]/sum(rpk[!(rpk==1-max.ar)]))
+- randomprob[as.numeric(colnames(rpk))] = randomprob[as.numeric(colnames(rpk))]+rpk
+

NAMESPACE

@@ -31,7 +31,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
                                   response.probs = c(0.4, 0.4),
                                   ns = c(30, 60, 90, 120, 150),
                                   max.ar = 0.75,
-                                  rand.type = "Urn",
+                                  rand.algo = "Urn",
                                   max.deviation = 3,
                                   model.inf = list(
                                     model = "tlr",
@@ -67,6 +67,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
                                 ),
                                 cl = 2) {
   old <- options()# code line i
+  on.exit(options(old))
   #Set start grid of screening
   startgrid <-
     data.frame(tpIE = rep(NA, length(grid.inf$start)), cutoff = grid.inf$start)
@@ -94,7 +95,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
       response.probs = input.info$response.probs,
       ns = input.info$ns,
       max.ar = input.info$max.ar,
-      rand.type = input.info$rand.type,
+      rand.algo = input.info$rand.algo,
       max.deviation = input.info$max.deviation,
       model.inf = input.info$model.inf,
       Stopbound.inf = Stopbound.inf,
@@ -121,7 +122,12 @@ demo_Cutoffscreening = function(ntrials = 1000,
   randomprobability = (1 / (abs(predictedtpIE[abs(predictedtpIE - 0.05) <=
                                                 0.0025] - 0.05) + e)) / sum(1 / (abs(predictedtpIE[abs(predictedtpIE - 0.05) <=
                                                                                                      0.0025] - 0.05) + e))
-  nextcutoff = sample(potentialcutoff, 1, replace = T, prob = randomprobability)
+  if (length(potentialcutoff) > 1){
+    nextcutoff = sample(potentialcutoff, 1, replace = T, prob = randomprobability)
+  }
+  else {
+    nextcutoff = cutoffgrid[which.min(abs(predictedtpIE - 0.05))]
+  }
   extendgrid[1, 2] = nextcutoff
   recommand = {
 
@@ -138,7 +144,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
       response.probs = input.info$response.probs,
       ns = input.info$ns,
       max.ar = input.info$max.ar,
-      rand.type = input.info$rand.type,
+      rand.algo = input.info$rand.algo,
       max.deviation = input.info$max.deviation,
       model.inf = input.info$model.inf,
       Stopbound.inf = Stopbound.inf,
@@ -164,9 +170,16 @@ demo_Cutoffscreening = function(ntrials = 1000,
       randomprobability = 1
       potentialcutoff = extendgrid[cutoffindex, 2]
     }
-    extendgrid[cutoffindex + 1, 2] = sample(potentialcutoff, 1, replace = T, prob = randomprobability)
-    recommand = c(recommand, cutoffgrid[as.numeric(names(which.max(randomprobability)))])
-    message(paste("Finished extend grid screening round", cutoffindex))
+    if (length(potentialcutoff) > 1){
+      extendgrid[cutoffindex + 1, 2] = sample(potentialcutoff, 1, replace = T, prob = randomprobability)
+      recommand = c(recommand, cutoffgrid[as.numeric(names(which.max(randomprobability)))])
+      message(paste("Finished extend grid screening round", cutoffindex))
+    }
+    else {
+      extendgrid[cutoffindex + 1, 2] = cutoffgrid[which.min(abs(predictedtpIE - 0.05))]
+      recommand = c(recommand, cutoffgrid[which.min(abs(predictedtpIE - 0.05))])
+      message(paste("Finished extend grid screening round", cutoffindex))
+    }
   }
   message("Output data recording")
   dataloginformd = data.frame(rbind(startgrid, extendgrid))
@@ -178,7 +191,6 @@ demo_Cutoffscreening = function(ntrials = 1000,
     predict(quadratic.model,
             list(cutoff = cutoffgrid, cutoff2 = cutoffgrid ^ 2))
   doParallel::stopImplicitCluster()
-  on.exit(options(old))
   return(
     list(
       detailsforgrid = dataloginformd,