R/prob.def2.R

In carolinewei/apsurvival: Total sample size and its allocation for the MRCT

#’ Calculate assurance probability by Definition 2 using simulation.
#‘
#‘ @description This function uses simulation to calculate the unconditional and
#' conditional probabilities of claiming consistency of treatments effects in an
#' MRCT with a survival endpoint by Definition 2. It is 
#' used in the article and other functions.
#'
#' @param r0 True overall log hazard ratio
#' @param s Number of regions participating in the MRCT
#' @param E0 Event number of all regions with two groups combined
#' @param u A vector presents ratios of true regional log
#' hazard ratios to true overall log hazard ratio r=u*r0
#' @param f A Vector presents proportions of total event number 
#' assigned to each region
#' @param bound The given parameter in Definition 2
#' @param n Simulation times
#' @param lamda The event hazard rate for placebo
#' @param lamda_cen The discontinuation hazard rate
#' @param L The whole study duration of fixed study duration 
#' design
#'
#' @return A list
#' @export
#' 
#' @examples
#' r0=log(0.7)
#' E0=Eventnum(r=r0, alpha=0.05, beta=0.2)
#' set.seed(123)
#' AP2 <- prob.def2(r0=r0, s=3, E0=E0, u=c(1,1,1), f=c(1/3,1/3,1/3),
#' bound=0.1, n=1000, lamda=1, lamda_cen=1, L=2)
#'

prob.def2 <- function(r0, s, E0, u, f, bound, n=100000, lamda, lamda_cen, L){
    
  E = f*E0
  r_reg = r0*u
  
  n.cn = 0
  count = 0
  count.cn = 0
  j = 0
  
  while(j < n){
      gamma <- rep(NA,s)
      for(k in 1:s){
          gamma[k] <- estimatelogHR(r=r_reg[k], lamda=lamda, lamda_cen=lamda_cen, L=L, E=E[k])
      }
      if(sum(is.nan(gamma))>0){
          j = j
          next
      }
      else{
          j = j+1
      }
      
      trt <- 1-exp(gamma)
      trt_glb <- 1-exp(sum(f*gamma))
      
      temp = sum( trt > bound )
      
      if(temp==s) count=count+1
      
      if( (temp==s) && (trt_glb>0) ) count.cn=count.cn+1
      
      if(trt_glb>0) n.cn=n.cn+1
      
  }
  
  prob.un = count/n
  prob.cn = count.cn/n.cn
  
  return( list(prob.un=prob.un,prob.cn=prob.cn) )
}