CRAN

Author: ZXW834

.DS_Store

@@ -1,13 +1,18 @@
 Package: BayesianPlatformDesignTimeTrend
 Title: Simulate and Analyse Bayesian Platform Trial with Time Trend
-Version: 1.1.3
-Author: Ziyan Wang [aut, cre]
+Version: 1.2.0
+Author: Ziyan Wang [aut, cre],
+  David Woods [ctb]
 Maintainer: Ziyan Wang <zw7g20@soton.ac.uk>
 Authors@R: 
-	person(given = "Ziyan",
+	c(person(given = "Ziyan",
            family = "Wang",
            role = c("aut", "cre"),
-           email = "zw7g20@soton.ac.uk")
+           email = "zw7g20@soton.ac.uk"),
+  	person(given = "David",
+           family = "Woods",
+           role = c("ctb"),
+           email = "D.Woods@southampton.ac.uk"))
 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. 
@@ -20,14 +25,14 @@ Description: Simulating the sequential multi-arm multi-stage or platform trial w
     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>.
-URL: https://github.com/ZXW834/PlatFormDesignTime
+URL: https://github.com/ZXW834/BayesianPlatformDesignTimeTrend
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Biarch: true
 Depends: 
     R (>= 4.0.0),
-    rstan (>= 2.18.1)
+    rstan (>= 2.26.1)
 Imports: 
     methods,
     Rcpp (>= 0.12.0),
@@ -43,14 +48,18 @@ Imports:
     reshape (>= 0.8.8),
     BiocManager (>= 1.30.19),
     lhs (>= 1.1.6),
-    laGP (>= 1.5-9)
+    laGP (>= 1.5-9),
+    RColorBrewer (>= 1.1-3),
+    directlabels (>= 2023.8.25),
+    ggpubr (>= 0.4.0),
+    reshape2 (>= 1.4.4)
 LinkingTo: 
     BH (>= 1.66.0),
     Rcpp (>= 0.12.0),
     RcppEigen (>= 0.3.3.3.0),
     RcppParallel (>= 5.0.1),
-    StanHeaders (>= 2.18.0),
-    rstan (>= 2.18.1)
+    StanHeaders (>= 2.26.0),
+    rstan (>= 2.26.1)
 SystemRequirements: GNU make, C++17
 License: MIT + file LICENSE
 Suggests: 

.Rhistory

@@ -3,7 +3,6 @@
 export(ARmethod)
 export(AdaptiveRandomisation)
 export(Boundaryconstruction)
-export(FWER_disconjunctivepowerfunc)
 export(GP.optim)
 export(Initializetrialparameter)
 export(Meanfunc)
@@ -16,19 +15,22 @@ export(Stopboundinf)
 export(Timetrend.fun)
 export(Trial.simulation)
 export(alphaspending)
+export(conjuncativepower_or_FWER)
 export(demo_Cutoffscreening)
 export(demo_Cutoffscreening.GP)
 export(demo_multscenario)
+export(disconjunctivepowerfunc)
 export(ibetabinomial.post)
 export(intbias)
 export(modelinf.fun)
-export(perHtypeIerror_powerfunc)
+export(perHtypeIerror_marginalpowerfunc)
 export(resultrtostats)
 export(resultrtostats.rand)
 export(resultstantoRfunc)
 export(resultstantoRfunc.rand)
 export(simulatetrial)
 export(stan.logisticmodeltrans)
+export(testing_and_armdropping)
 export(trtbias)
 export(trteffect)
 export(varfunc)
@@ -44,13 +46,18 @@ importFrom(boot,logit)
 importFrom(doParallel,registerDoParallel)
 importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
+importFrom(grDevices,heat.colors)
+importFrom(graphics,image)
 importFrom(graphics,lines)
+importFrom(graphics,par)
+importFrom(graphics,points)
 importFrom(iterators,icount)
 importFrom(laGP,distance)
 importFrom(matrixStats,colVars)
 importFrom(rstan,rstan_options)
 importFrom(rstan,sampling)
 importFrom(stats,dbeta)
+importFrom(stats,dist)
 importFrom(stats,integrate)
 importFrom(stats,lm)
 importFrom(stats,model.matrix)
@@ -60,5 +67,6 @@ importFrom(stats,pnorm)
 importFrom(stats,predict)
 importFrom(stats,qnorm)
 importFrom(stats,rbinom)
+importFrom(stats,runif)
 importFrom(stats,var)
 useDynLib(BayesianPlatformDesignTimeTrend, .registration = TRUE)

DESCRIPTION

@@ -1,17 +1,41 @@
+# Sept 30, 2023, BayesianPlatformDesignTimeTrend version 1.2.0
+## Summary of changes
+- This release corrects the boundary construction formula of OBF boundary; 
+- This release adds internal functions of calculating conjunctive power, disconjunctive power and marginal power.
+- This release extend the type of hypothesis testing which are Oneside test and the two side test now; 
+- This release implements method of cutoff screening for asymmetric boundary using Bayesian optimization. The cutoff screening for asymmetric boundary is released in this version (including No early stopping boundary, Pocock boundary and OBF boundary). The plot of boundary shape refers to the plot on github: https://github.com/ZXW834/BayesianPlatformDesignTimeTrend.
+- This release implements method of hyperparameter tuning for Trippa's adaptive randomisation approach. There are some flow chats in the tutorial and on the github.
+- This release adds my supervisor David Woods as the contributor to the package.
+
+## Bugs
+- In OBF boundary construction, the boundary formula was $\theta_j=\phi \left(\sqrt{\frac{J}{j}c} \right)$ which is not the same as that in the reference which is $\theta_j=\phi \left(\sqrt{\frac{J}{j}}c \right)$. Although the original boundary formula works in trial simulation, the correction is still made.
+
+## New feature
+- In previous version, the power function was only for disconjunctive power and FWER. In this version, the conjunctive power `conjuncativepower_or_FWER` and marginal power (per-hypothesis type I error)  `perHtypeIerror_marginalpowerfunc` can also be calculated for the used of design evaluation and power optimisation in asymmetric boundary cutoff screening. The function to calculate disconjunctive power is renamed as `disconjunctivepowerfunc`.
+
+- In previous version, the hypothesis testing during the interim analysis is $H_0$: $\pi_k = \pi_0$; $H_1$: $\pi_k \neq \pi_0$. Now we add one side test to only conclude for superiority: $H_0$: $\pi_k \leq \pi_0$; $H_1$: $\pi_k > \pi_0$.
+
+- In previous version, the cutoff value can only be tuned for symmetric boundary. In this version, asymmetric boundary cutoff screening is added so that user can decide whether to use a more aggressive efficacy boundary with a more aggressive futility boundary (or a more conservative efficacy boundary with more conservative futility boundary). The asymmetric boundary screening process first select the type I error rate contour equal to the target, then pick a cutoff pair (eff, fut) that maximizes the power.
+
+- In previous version, the Trippa's adaptive randomisation approach used fixed hyperparameters suggested in their paper. In this version, the hyperparameters a and b are tuned to optimise power when type I error (FWER) is controlled at target. This idea can be found on Trippa's paper (2012). We used a gaussian process model to smooth the contour plot. The pair (a, b) with highest predicted power will be exploit for hyperparameter screening until there are enough point overlapped on the contour plot. In addition, the hyperparameter point (a, b) can be pick in a crude way to save computation time via setting the maximum number of points to be looked at.
+
 # Aug 22, 2023, BayesianPlatformDesignTimeTrend version 1.1.3
 This release fix one command in fix effect model analysis of main effect that is confusing.
 
-- I fixed the bug in function `resultrtostats`. THe calculation of "stats4" was the sum of main effect and interaction effect when using the model with interaction term (main * stage_continuous / main * stage_discrete). The sum was used to calculate the posterior probability of superior or inferior. However, "stats4" only represents the mean values of main effect in the final output matrix. The wrong value was given to "stats4". This is the same for "stats5" which is the standard errors of the main effect in the model
+## Bugs
+- I fixed the bug in function `resultrtostats`. The calculation of "stats4" was the sum of main effect and interaction effect when using the model with interaction term (main * stage_continuous / main * stage_discrete). The sum was used to calculate the posterior probability of superior or inferior. However, "stats4" only represents the mean values of main effect in the final output matrix. The wrong value was given to "stats4". This is the same for "stats5" which is the standard errors of the main effect in the model
 
 # Jun 25, 2023, BayesianPlatformDesignTimeTrend version 1.1.2
 This release fix two bugs.
 
+## 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.
 
+## Bugs
 - 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.
 

NAMESPACE

@@ -2,13 +2,13 @@
 #' @description This function reads in the output matrix of a number of trial replicates to calculate the error rate
 #' @param res A list of output matrix of a number of trial replicates
 #'
-#' @return Type I error rate or power of each treatment - control comparison
+#' @return "per-hypothesis" Type I error rate or power (marginal power) for each treatment - control comparison
 #' @export
 #'
 #' @examples
 #' \dontrun{perHtypeIerror_powerfunc(res)}
 #' @author Ziyan Wang
-perHtypeIerror_powerfunc = function(res) {
+perHtypeIerror_marginalpowerfunc = function(res) {
   colMeans(t(sapply(res, function(x) {
     resname = colnames(x)
     K = sum(stringr::str_detect(colnames(x), "H")) + 1
@@ -25,16 +25,16 @@ perHtypeIerror_powerfunc = function(res) {
     }
   })))
 }
-#' @title FWER_disconjunctivepowerfunc
+#' @title disconjunctivepowerfunc
 #' @description This function reads in the output matrix of a number of trial replicates to calculate the Family wise error rate or disconjunctive power
 #' @param res A list of output matrix of a number of trial replicates
 #'
-#' @return Family wise error rate or disconjunctive power
+#' @return Disconjunctive power
 #' @export
 #'
 #' @examples
-#' \dontrun{FWER_disconjunctivepowerfunc(res)}
-FWER_disconjunctivepowerfunc = function(res) {
+#' \dontrun{disconjunctivepowerfunc(res)}
+disconjunctivepowerfunc = function(res) {
   mean(sapply(res, function(x) {
     resname = colnames(x)
     K = sum(stringr::str_detect(colnames(x), "H")) + 1
@@ -47,6 +47,113 @@ FWER_disconjunctivepowerfunc = function(res) {
   }))
 }
 
+#' @title conjuncativepower_or_FWER
+#'
+#' @description This function reads in the output matrix of a number of trial replicates to calculate the Family wise error rate or Conjunctive power
+#'
+#' @param res A list of output matrix of a number of trial replicates
+#' @param scenario The true scenario used to generate the res list
+#' @param test.type The indicator of whether using one side or two side testing.
+#'     Please make sure that the input test.type does not conflicts to the data. Otherwise the conjunctive power calculation is wrong
+#'
+#' @return Family wise error rate or Conjunctive power
+#' @export
+#'
+#' @examples
+#' \dontrun{conjuncativepower_or_FWER(res)}
+conjuncativepower_or_FWER=function(res,scenario,test.type){
+  K=mean(sapply(res,function(x){
+    K=sum(stringr::str_detect(colnames(x),"H"))+1
+    return(K)
+  }))
+  hypres=mean(sapply(res,function(x){
+    stage=dim(x)[1]
+    resname=colnames(x)
+    K=sum(stringr::str_detect(colnames(x),"H"))+1
+
+    # True scenario construction
+    if (sum(scenario-min(scenario))>0 & test.type == "Twoside"){
+      sce = (scenario-scenario[1])[-1]!=0
+    }
+    else if (sum(scenario-min(scenario))>0 & test.type == "Oneside"){
+      sce = (scenario-scenario[1])[-1]>0
+    }
+    else{
+      sce = rep(0,K-1)
+    }
+
+    # # Check whether the data output conflicts the test type. This will lead to error if the hypothesis test data input are all zero under the alternative
+    # if (test.type == "Twoside"){
+    #   test.data = x[,(K-1+2*K+1):(K-1+2*K+K-1)]
+    #   superiorarm.loc = which((scenario-scenario[1])[-1]>0)
+    #   inferiorarm.loc = which((scenario-scenario[1])[-1]<0)
+    #   for (i in 1:length(inferiorarm.loc)){
+    #     if (test.data[min(which(is.na(test.data[,inferiorarm.loc[i]])))-1,inferiorarm.loc[i]] == 0){
+    #       stop("The data type used a Two side test while the input variable is one side.")
+    #     }
+    #   }
+    #   for (i in 1:length(superiorarm.loc)){
+    #     if (test.data[min(which(is.na(test.data[,superiorarm.loc[i]])))-1,superiorarm.loc[i]] == 0){
+    #       stop("The data type used a Two side test while the input variable is one side.")
+    #     }
+    #   }
+    # }
+    # else if (test.type == "Oneside"){
+    #   inferiorarm.loc = which((scenario-scenario[1])[-1]<0)
+    #   test.data = x[,(K-1+2*K+1):(K-1+2*K+K-1)]
+    #   # check if any inferior arm is identified to be inferior for one side test
+    #   for (i in 1:length(inferiorarm.loc)){
+    #     if (test.data[min(which(is.na(test.data[,inferiorarm.loc[i]])))-1,inferiorarm.loc[i]] == 1){
+    #     stop("The data type used a Two side test while the input variable is one side.")
+    #     }
+    #     }
+    # }
+    # else{
+    #   Stop("Input test.type is invalid.")
+    # }
+    #Indentify which hypothesis is rejected
+    reject=which(matrix(x[,(K-1+2*K+1):(K-1+2*K+K-1)] %in% 1,ncol=K-1),1)[,2]
+    drop.at=which(matrix(x[,(K-1+2*K+1):(K-1+2*K+K-1)] %in% 1,ncol=K-1),1)[,1]
+    drop.at.all=rep(stage,K-1)
+    drop.at.all[reject]=drop.at
+    treatmentindex=seq(1,K-1)
+    trtmean.loc=cbind(drop.at.all,treatmentindex)
+
+    res=matrix(x[,(K-1+2*K+1):(K-1+2*K+K-1)] %in% 1,ncol=K-1)
+    result=rep(NA,K-1)
+    for (i in 1:(K-1)){
+      result[i] = res[trtmean.loc[i,1],trtmean.loc[i,2]]
+    }
+    if (sum(result == sce) == (K-1)){
+      return(1)
+    }
+    else{
+      return(0)
+    }
+  }))
+  if (test.type == "Twoside"){
+    if (sum(scenario-min(scenario))>0){
+    current_scenario = "Alt"
+    return(hypres)
+  }
+  else{
+    current_scenario = "Null"
+    return(1-hypres)
+  }
+  }
+  else{
+    if ((max(scenario)-scenario[1]) <= 0){
+      current_scenario = "Null"
+      return(1-hypres)
+    }
+    else{
+      current_scenario = "Alt"
+      return(hypres)
+    }
+  }
+
+}
+
 #' @title Meanfunc
 #' @description This function reads in the output matrix of a number of trial replicates to calculate mean treatment effect estimate.
 #' @param res A list of output matrix of a number of trial replicates

NEWS.md

@@ -33,6 +33,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
                                   max.ar = 0.75,
                                   rand.algo = "Urn",
                                   max.deviation = 3,
+                                  test.type = "Twoside",
                                   model.inf = list(
                                     model = "tlr",
                                     ibb.inf = list(
@@ -104,7 +105,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
     )
 
     # perHtypeIerror=mean(perHtypeIerrorfunc(result))
-    FWER = FWER_disconjunctivepowerfunc(result)
+    FWER = conjuncativepower_or_FWER(result,input.info$response.probs,test.type = input.info$test.type)
     startgrid[j, 1] = FWER
   }
 
@@ -151,7 +152,7 @@ demo_Cutoffscreening = function(ntrials = 1000,
       Random.inf = input.info$Random.inf,
       trend.inf = input.info$trend.inf
     )
-    FWER = FWER_disconjunctivepowerfunc(restlr090five)
+    FWER = conjuncativepower_or_FWER(restlr090five, input.info$response.probs,test.type = input.info$test.type)
     extendgrid[cutoffindex, 1] = FWER
     extendgrid$cutoff2 <- extendgrid$cutoff ^ 2
     quadratic.model <-