R/prob.def3.R

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

#’ Calculate assurance probability by Definition 3 using asymptotic distributions derived.
#‘
#‘ @description This function uses derived distributions to calculate the unconditional
#' and conditional probabilities of claiming consistency of treatments effects in
#' an MRCT with a survival endpoint by Definition 3.
#'
#' @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 eps Significance level of not rejecting H0 in Definition 3
#' @param pai The preserved proportion of overall treatment effect
#' @param alpha The risk of rejecting the null hypothesis H0:r0>=0
#' when it is really true
#'
#' @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)
#' AP3 <- prob.def3(r0=r0, s=3, E0=E0, u=c(1,1,1), f=c(1/3,1/3,1/3),
#' eps=0.3,pai=1/3, alpha=alpha)
#'
prob.def3 <- function(r0, s, E0, u, f, eps, pai, alpha){
    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)
    
    mean3 = delta_reg-pai*delta_glb
    
    v1 = (pai^2)*exp(2*r_glb)*(4/E0)
    v2 = 4*pai*exp(r_reg+r_glb)/E0
    v3 = 4/E*exp(2*r_reg)
    e = rep(1,s)
    tmp_covar3 = -e%*%t(v2)-v2%*%t(e)+v1*e%*%t(e)
    covar3 = tmp_covar3
    diag(covar3) = v3-2*v2+v1
    
    z.eps = qnorm(1-eps)
    lower = z.eps*sqrt(v3-2*v2+v1)
    upper = rep(Inf,s)
    
    prob.un <- pmvnorm(lower,upper,mean=mean3,sigma = covar3)
    
    mean3_1 = c(mean3,delta_glb)
    
    v4 = 4*exp(r_reg+r_glb)/E0-pai*exp(2*r_glb)*(4/E0)
    covar3_1 = matrix(0,(s+1),(s+1))
    covar3_1[1:s,1:s] = covar3
    covar3_1[(s+1),1:s] = v4
    covar3_1[1:s,(s+1)] = v4
    covar3_1[(s+1),(s+1)] = exp(2*r_glb)*(4/E0)
    
    z_alpha = qnorm(1-alpha) 
    lower_1 = c(lower,z_alpha*sqrt(covar3_1[(s+1),(s+1)]))
    upper_1 = rep(Inf,(s+1))
    prob.cn <- pmvnorm(lower=lower_1,upper=upper_1,mean=mean3_1,
    sigma=covar3_1)/ (1-pnorm(z_alpha*sqrt(covar3_1[(s+1),(s+1)]),
    mean=delta_glb,sd=sqrt(covar3_1[(s+1),(s+1)])))
    
    return(list(prob.un=prob.un[1],prob.cn=prob.cn[1]))
}