R/sample.def1.R

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

#’ Calculate total sample size and its optimal allocation by Definition 1
#' using simulation.
#‘
#‘ @description This function calculates total sample size with and without 
#' considering desired conditional assurance probability to claim overall consistency and 
#' determines optimal sample size allocation across regions by maximizing conditional 
#' assurance probability based on Definition 1 if regional treatment effects 
#' are slightly different. (Allocate equal sample size to each region 
#' if treatment effects across regions are the same.)
#'
#' @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 beta The risk of failing to reject the null
#' hypothesis H0:r0>=0 when it is really false
#' @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
#' @param s Number of regions participating in the MRCT
#' @param u A vector presents ratios of true regional log
#' hazard ratios to true overall log hazard ratio r=u*r0
#' @param pai The given parameter in Definition 1
#' @param grid Grid interval of the grid research
#' @param n Simulation times
#' @param consistency A numeric value is the desired conditional assurance
#' probability to claim overall consistency showing only two decimal places.
#'
#' @return A list
#' @export
#'
#' @examples
#' set.seed(123)
#' Sampsize1 <- sample.def1(r0=log(0.7), alpha=0.05, beta=0.2, lamda=1, lamda_cen=1, L=2,
#' s=3, u=c(0.9,1,1.1), pai=1/3, grid=0.1, n=1000, consistency=0.8)
#'
sample.def1 <- function(r0, alpha=0.05, beta=0.2, lamda, lamda_cen, L, s, u, pai, grid=0.1, n=100000, consistency=0.8){
    if(consistency!=round(consistency,2)) stop("The value only shows two decimal places.")
    t1 <- Sys.time()
    E0 = 4*(qnorm(alpha,0,1)+qnorm(beta,0,1))^2/(r0)^2
    sampsize <- Samplesize(r0, alpha=alpha, beta=beta, E=ceiling(E0), lamda, lamda_cen, L)
    if(sum(u==rep(1,s))==s){
        f = rep(1/s,s)
        prob1 <- prob.def1(r0=r0, s=s, E0=ceiling(E0), u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
        final <- as.data.frame(t(f))
        colnames(final) <- paste0("f",1:s)
        final$con.AP <- prob1
    }else{
        final <- alloc.def1(r0=r0, alpha=alpha, beta=beta, lamda=lamda, lamda_cen=lamda_cen, L=L, s=s, u=u, pai=pai, grid=grid, n=n)$alloc_cn.AP
        f <- as.numeric(final[1,1:s])
        prob1 <- final$con.AP[1]
    }
    if(prob1 < (consistency-0.005)){
        fold <- ceiling((consistency-prob1)/0.1)+1
        prob3 <- prob.def1(r0=r0, s=s, E0=ceiling(fold*E0), u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
        while(prob3 < (consistency-0.005)){
            fold <- fold * (ceiling((consistency-prob3)/0.1)+1)
            prob3 <- prob.def1(r0=r0, s=s, E0=ceiling(fold*E0), u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
        }
        fold0 = 1
        fold1 = (fold0+fold)/2
        prob2 <- prob.def1(r0=r0, s=s, E0=ceiling(fold1*E0), u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
        while(prob2<(consistency-0.005)|prob2>=(consistency+0.005)){
            if(prob2>=(consistency+0.005)){
                fold = fold1
                prob3 = prob2
                fold1 = (fold0+fold)/2
            }
            if(prob2<(consistency-0.005)){
                fold0 = fold1
                prob1 = prob2
                fold1 = (fold0+fold)/2
            }
            prob2 <- prob.def1(r0=r0, s=s, E0=ceiling(fold1*E0), u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
        }
        for(i in 1:ceiling(fold1*E0)){
            prob <- prob.def1(r0=r0, s=s, E0=ceiling(fold1*E0)-i, u=u, f=f, pai=pai, n=n, lamda=lamda, lamda_cen=lamda_cen, L=L)$prob.cn
            if(prob < (consistency-0.005)) break
        }
        sampsize_AP <- Samplesize(r0, alpha=alpha, beta=beta, E=ceiling(fold1*E0)-i+1, lamda, lamda_cen, L)
    }else{
        sampsize_AP <- sampsize
    }
    t2 <- Sys.time()
    t2 <- Sys.time()
    full_list <- list(sampsize_AP, sampsize, f, final, t2-t1)
    names(full_list) <- c("samplesize_AP","samplesize","optimal_alloc","alloc_cn.AP","duration")
    return(full_list)
}