R/estimatelogHR.R

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

#’ Estimate log hazard ratio based on logrank test statisitic and
#' simulating survival data under model hazard_1(t)=hazard_0(t)*e^r.
#'
#' @param r True log hazard ratio
#' @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 E Event number of two groups combined
#'
#' @return A single value.
#' @export
#'
#' @examples
#' r0=log(0.7)
#' E0=Eventnum(alpha=0.05, beta=0.2, r=r0)
#' set.seed(123)
#' estimatelogHR(r=r0, lamda=1, lamda_cen=1, L=2, E=E0)
#'
estimatelogHR <- function(r, lamda=1, lamda_cen=1, L=2, E){
    lamda1 = lamda*exp(r)
    H0 = lamda/(lamda+lamda_cen)*(1-exp(-L*(lamda+lamda_cen)))
    H1 = lamda1/(lamda1+lamda_cen)*(1-exp(-L*(lamda1+lamda_cen)))
    N = round(2*E/(H0+H1))
    if (N %% 2 == 1) N = N+1
    ratio = 0.5
    n = N*ratio
    
    T0 <- round(rexp(n,lamda),2)
    T1 <- round(rexp(n,lamda*exp(r)),2)
    T_sur <- c(T0,T1)
    cen <- round(rexp(N,lamda_cen),2)
    T_cen <- pmin(cen, L)
    y <- pmin(T_sur, T_cen)
    status <- as.numeric(T_sur > T_cen)
    ifdeath <- 1-status
    eventnum <- N-sum(status)
    
    g0 <- as.data.frame( cbind(y[1:n],status[1:n],ifdeath[1:n]) )
    colnames(g0) <- c("survTime", "ifcensor", "ifdeath")
    g1 <- as.data.frame( cbind(y[(n+1):N],status[(n+1):N],ifdeath[(n+1):N]) )
    colnames(g1) <- c("survTime", "ifcensor", "ifdeath")
    combine <- rbind(g0,g1)
    
    event <- combine[which(combine[,2]==0),1]
    event <- sort(event)
    event <- unique(event)
    
    o <- rep(NA,length(event))
    e <- rep(NA,length(event))
    v <- rep(NA,length(event))
    
    for(i in 1: length(event)){
        tao = event[i]
        d0 = sum(g0$ifdeath[which(g0$survTime==tao)])
        d1 = sum(g1$ifdeath[which(g1$survTime==tao)])
        y0 = n-sum(g0$ifdeath[which(g0$survTime<tao)])-sum(g0$ifcensor[which(g0$survTime<tao)])
        y1 = n-sum(g1$ifdeath[which(g1$survTime<tao)])-sum(g1$ifcensor[which(g1$survTime<tao)])
        
        o[i] = d1
        e[i] = y1/(y0+y1)*(d0+d1)
        v[i] = y0*y1*(d0+d1)*(y0+y1-d1-d0)/(y0+y1)^2/(y0+y1-1)
    }
    
    v <- v[!is.nan(v)]
    z <- sum(o-e)/sqrt(sum(v))
    loghr <- 2*z/sqrt(eventnum)
    return(loghr)
}