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.0.1
 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 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)
@@ -41,8 +41,7 @@ Imports:
     foreach (>= 1.5.1),
     iterators (>= 1.0.13),
     reshape (>= 0.8.8),
-    BiocManager (>= 1.30.19),
-    rjags
+    BiocManager (>= 1.30.19)
 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

@@ -6,8 +6,6 @@ export(Boundaryconstruction)
 export(FWER_disconjunctivepowerfunc)
 export(Initializetrialparameter)
 export(Meanfunc)
-export(Mixeffect_analysis)
-export(Mixeffect_modelling)
 export(Nfunc)
 export(OutputStats.initialising)
 export(Randomisation.inf)
@@ -57,6 +55,5 @@ importFrom(stats,pbeta)
 importFrom(stats,pnorm)
 importFrom(stats,predict)
 importFrom(stats,rbinom)
-importFrom(stats,update)
 importFrom(stats,var)
 useDynLib(BayesianPlatformDesignTimeTrend, .registration = TRUE)

LICENSE

@@ -1,2 +1,48 @@
-# BayesianPlatformDesignTimeTrend 1.0.0
-* Added a `NEWS.md` file to track changes to the package.
+# 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
+
+# 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)

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,
@@ -138,7 +139,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,
@@ -178,7 +179,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,

NEWS.md

@@ -66,7 +66,7 @@ demo_multscenario = function(ntrials = 1000,
           response.probs = scenario[i,],
           ns = ns[[z]],
           max.ar = 0.75,
-          rand.type = "Urn",
+          rand.algo = "Urn",
           max.deviation = 3,
           model.inf = list(
             model = "tlr",

R/Demo_CutoffScreening.R

@@ -3,8 +3,8 @@
 #' @param ii Meaning less parameter but required for foreach function in doParallel package
 #' @param response.probs A vector of true response probability for each arm. Default response.probs = c(0.4, 0.4).
 #' @param ns  A vector of accumulated number of patient at each stage. Default is ns = c(30, 60, 90, 120, 150).
-#' @param max.ar The upper boundary for randomisation ratio for each arm. Default is 0.75 for a two arm trial.
-#' @param rand.type The method of applying patient allocation with a given randomisation probability vector. Default is "Urn".
+#' @param max.ar The upper boundary for randomisation ratio for each arm. Default is 0.75 for a two arm trial. The minimum value depends on K where 1 - max.ar <= 1/K
+#' @param rand.algo The method of applying patient allocation with a given randomisation probability vector. Default is "Urn".
 #' @param max.deviation The tuning parameter for Urn randomisation method. Default is 3.
 #' @param Stopbound.inf The list of stop boundary information for more see \code{\link{Stopboundinf}}
 #' @param Random.inf The list of Adaptive randomisation information for more see \code{\link{Randomisation.inf}}
@@ -19,7 +19,7 @@
 #' simulatetrial(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",
@@ -60,7 +60,7 @@ simulatetrial <- function(ii,
                            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",
@@ -113,9 +113,14 @@ simulatetrial <- function(ii,
   randomprob = initialised.par$randomprob
   z = initialised.par$z
   y = initialised.par$y
-
   group_indicator = initialised.par$group_indicator
   post.prob.best.mat = initialised.par$post.prob.best.mat
+
+  #-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")
+  }
+
   #Initialize the output data frame: stats
   stats = OutputStats.initialising(model.inf$tlr.inf$variable.inf,
                                    model.inf$tlr.inf$reg.inf,
@@ -139,10 +144,10 @@ simulatetrial <- function(ii,
 
     # Patient randomisation based on given randomization ratio
     # Parameter checking
-    if (rand.type %in% c("Coin", "Urn")) {
+    if (rand.algo %in% c("Coin", "Urn")) {
       random.output = AdaptiveRandomisation(
         Fixratio,
-        rand.type,
+        rand.algo,
         K,
         n.new,
         randomprob,
@@ -159,7 +164,7 @@ simulatetrial <- function(ii,
       )
     }
     else {
-      stop("Error randomisation type wrong")
+      stop("Error: randomisation type wrong")
     }
 
     # Extract the output list
@@ -491,66 +496,6 @@ simulatetrial <- function(ii,
           treatmentindex = treatmentindex[is.na(match(treatmentindex, treatmentdrop))]
         }
       }
-      else if (model.inf$tlr.inf$variable.inf == "Mixeffect.jag") {
-        jagmodel<-"  model {
-            for(i in 1:armleft){
-              for(j in 1:stage){
-                Y[i,j] ~ dbin(p[i,j],N[i,j])
-                logit(p[i,j]) = beta0 + alpha[stage-(j-1)] + beta1[i]
-              }
-            }
-           alpha[1] = 0
-           alpha[2] ~ dnorm(0, tau2)
-           for(k in 3:stage ) {
-             alpha[k] ~ dnorm(2*alpha[k-1] - alpha[k-2],tau2)
-        }
-           beta1[1] <- 0
-           for(i in 2:armleft){
-             beta1[i] ~ dnorm(0,0.31)
-           }
-              tau2 ~ dgamma(0.1, 0.01)
-           beta0 ~ dnorm(0, 0.31)
-          }"
-        postsamp.list = Mixeffect_modelling(ytemp=ytemp, treatmentindex=treatmentindex, group=group, ntemp=ntemp, armleft=armleft, jagmodel=jagmodel)
-        analysis = Mixeffect_analysis(postsamp.list=postsamp.list, group=group, treatmentindex=treatmentindex, ns = ns, K = K)
-        stats1 = analysis$stats1
-        statsbeta0 = analysis$statsbeta0
-        stats4 = analysis$stats4
-        stats5 = analysis$stats5
-        stats6 = analysis$stats6
-        stats7 = analysis$stats7
-        sampoutcome = analysis$sampoutcome
-        sampefftotal = analysis$sampefftotal
-        post.prob.btcontrol = analysis$post.prob.btcontrol
-        #-Calculating posterior probability of each arm (including control) to be the best arm-
-        # post.prob.best: The posterior probability of each arm (including control) to be the best arm
-        # This is required for Thall's randomisation approach
-        for (q in 1:armleft) {
-          post.prob.best.mat[group, zlevel[q]] = (sum(max.col(sampefftotal) == q)) /
-            2500
-        }
-        post.prob.best = post.prob.best.mat[group,]
-        #Normalizing in case any value equals zero
-        post.prob.best = post.prob.best + 1e-7
-        post.prob.best = post.prob.best / sum(post.prob.best)
-        #----Justify if type I error was made for each arm----
-        # post.prob.btcontrol>cutoffeff[group]: Efficacy boundary is hit at this stage
-        # post.prob.btcontrol<cutoffful[group]: Fultility boundary is hit at this stage
-        Justify = post.prob.btcontrol > cutoffeff[group] |
-          post.prob.btcontrol < cutoffful[group]
-        #Identify which active arm should be dropped at current stage
-        treatmentdrop = treatmentindex[post.prob.btcontrol > cutoffeff[group] |
-                                         post.prob.btcontrol < cutoffful[group]]
-        stats3 = rep(NA, K - 1)
-        names(stats3) = seq(1, K - 1)
-        stats3[treatmentindex] = Justify
-        if (sum(Justify) > 0) {
-          armleft = armleft - sum(Justify)
-          #Debugged for K arm by Ziyan Wang on 12:00 26/07/2022 for three arm. Used to be treatmentindex = treatmentindex[-treatmentdrop]
-          #Debugged for K arm by Ziyan Wang on 18:58 26/07/2022 for more than 3 arm. Used to be treatmentindex = treatmentindex[!(treatmentindex==treatmentdrop)]
-          treatmentindex = treatmentindex[is.na(match(treatmentindex, treatmentdrop))]
-        }
-      }
 
       #-Adjust the posterior randomisation ratio-
       randomprob = ARmethod(

R/Demo_multiplescenariotrialsimulation.R

@@ -52,6 +52,13 @@ ARmethod = function(Fixratio,
                     max.ar,
                     armleft,
                     treatmentindex) {
+  # Validate inputs
+  if (!is.logical(Fixratio)) stop("Error: Fixratio should be a logical value (TRUE/FALSE)")
+  if (!is.character(BARmethod)) stop("Error: BARmethod should be a character value")
+  if (!is.numeric(group) || group <= 0) stop("Error: group should be a positive numeric value")
+  if (!is.numeric(max.ar) || max.ar <= 0 || max.ar >= 1) stop("Error: max.ar should be a numeric value between 0 and 1")
+  if (K < 2) stop("Error: K should be an integer value greater than or equal to 2")
+
   #---------------------Trippa's approach---------------------
   if (Fixratio == F & BARmethod == "Trippa") {
     ##Tuning the paprameter using method mentioned in Trippa's paper (2014)
@@ -97,20 +104,19 @@ ARmethod = function(Fixratio,
     }
     #-------------------Allocation bounds restriction (K arm: restriction on control)---------------
     else if (K > 2) {
-      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 = 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
+      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
     }
 
   }