R/prob.def4.R

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

#’ Calculate assurance probability by Definition 4 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 4.
#'
#' @param r0 True overall log hazard ratio
#' @param alpha The risk of rejecting the null hypothesis H0:r0>=0
#' when it is really true
#' @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 eps Significance level of not rejecting H0 in Definition 4
#' @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
#' alpha=0.05
#' r0=log(0.7)
#' E0=Eventnum(r=r0, alpha=0.05, beta=0.2)
#' set.seed(123)
#' AP4 <- prob.def4(r0=r0, alpha=alpha, s=3, E0=E0, u=c(1,1,1), f=c(1/3,1/3,1/3),
#' eps=0.1, n=1000, lamda=1, lamda_cen=1, L=2)
#'
prob.def4 <- function(r0, alpha, s, E0, u, f, eps=0.1, n=100000, lamda, lamda_cen, L){
    library(mvtnorm)
    E = f *E0
    r_reg = r0*u
    r_glb = sum(r_reg*f)
    delta_reg = 1-exp(r_reg)
    delta_glb = 1-exp(r_glb)
    
    covar3 = matrix(0,s,s)
    diag(covar3) = 4*exp(2*r_reg)/E
    sigma_delta_glb = sqrt(exp(2*r_glb)*(4/E0))
    
    n_cn = 0
    reject.un = 0
    reject.cn = 0
    z.eps = qchisq(1-eps, df=s)
    z.alpha = qnorm(1 - alpha)
    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))
        Z = trt - trt_glb
        Q = t(Z)%*%solve(covar3)%*%Z
        
        if(Q>z.eps) reject.un = reject.un+1
        
        if(trt_glb>z.alpha*sigma_delta_glb){
            n_cn = n_cn+1
            if(Q>z.eps) reject.cn = reject.cn+1
        }
    }
    AP.un = 1-reject.un/n
    AP.cn = 1-reject.cn/n_cn
    return( list(prob.un=AP.un,prob.cn=AP.cn) )
}