R/corrTime.R

Defines functions corrTime

Documented in corrTime

#' Correlation Between Two Logrank Test Statistics Over Time With Staggered Entry
#' 
#' This function calculates the correlation of two logrank test statistics on overlapping populations A and B
#' over time calculated from first subject in. 
#' 
#' @param T  Analysis time calculated from first subject randomization date.
#' @param r  Proportion of experimental subjects in both population A and B, 
#' A not B, B not A. Default, 1/2 for 1:1 randomization stratified by A and B
#' @param h0 Hazard function of control arm for subjects in both population A and B, 
#' A not B, B not A. h0(t) = log(2)/m0 means T~exponential distribution with 
#' median m0. For study design without considering heterogeneous effect in 
#' strata for control arm, then specify the same h0(t) function across strata.
#' @param S0 Survival function of control arm for subjects in both population A and B, 
#' A not B, B not A. h0(t) = log(2)/m0 means T~exponential distribution with 
#' median m0. For study design without considering heterogeneous effect in 
#' strata for control arm, then specify the same S0(t) function across strata.
#' The density function f0(t) = h0(t) * S0(t).
#' @param h1 Hazard function of experimental arm for subjects in both population A and B, 
#' A not B, B not A. For study design without considering heterogeneous effect in 
#' strata for the experimental arm, then specify the same h1(t) function across strata.
#' @param S1 Survival function of experimental arm for subjects in both population A and B, 
#' A not B, B not A. For study design without considering heterogeneous effect in 
#' strata for the experimental arm, then specify the same h1(t) function across strata.
#' @param F.entry Distribution function of enrollment. For uniform enrollment, 
#' F.entry(t) = (t/A) where A is the enrollment period, i.e., F.entry(t) = t/A for 0<=t<=A, and 
#' F.entry(t) = 1 when t > A. For more general non-uniform enrollment with weight psi, 
#' F.entry(t) = (t/A)^psi*I(0<=t<=A) + I(t>A). Default F.entry is uniform distribution function.
#' @param G.ltfu Distribution function of lost-to-follow-up censoring process. The observed
#' survival time is min(survival time, lost-to-follow-up time). For 3\% drop-off 
#' every year as a constant rate, then G.ltfu = 1-exp(-0.03/12*t), i.e., drop-off~exp distribution.
#' Default G.ltfu = 0 (no lost-to-followup)
#' @param n Total sample size for two arms for subjects in both population A and B, 
#' A not B, B not A. Default is NULL. 
#' @param strat.ana Stratified Analysis by A and B "Y" or "N". Default "Y"
#'
#' @return An object with dataframes below.
#'  \describe{
#'  \item{n.events}{Expected number of events}
#'       \itemize{
#'       \item n.events0: number of events for control group
#'       \item n.events1: number of events for experimental group
#'       \item n.events.total: total number of events for two groups
#'       \item subgroup  Label of subgroups
#'       }
#'  \item{corr}{Correlation of the two test statistics for populations A and B at time T}
#'  }
#'  
#' @examples
#' ################################
#' # Example 1: PD-L1+ subgroup 
#' # and overall population tests.
#' ################################
#' #Control arm has exponential
#' #distribution with median 12 months in all strata. 
#' #1:1 randomization stratified by PD-L1+ status. 300 subjects in PD-L1+ and 
#' #450 total subjects in overall population. 3\% drop-off every year. Enrollment
#' #period is 18 months and weight 1.5.
#' 
#' #################
#' # h(t) and S(t)
#' #################
#' 
#' h0 = function(t){log(2)/12}
#' S0 = function(t){exp(-log(2)/12 * t)}
#' h0.5 = function(t){log(2)/12*0.5}
#' h0.65 = function(t){log(2)/12*0.65}
#' h0.8 = function(t){log(2)/12*0.8}
#' S0.5 = function(t){exp(-log(2)/12 * 0.5 * t)}
#' S0.65 = function(t){exp(-log(2)/12 * 0.65 * t)}
#' S0.8 = function(t){exp(-log(2)/12 * 0.8 * t)}
#' 
#' #Entry distribution: enrollment period 18 mo, acceleration weight 1.5.
#' Fe = function(t){(t/18)^1.5*as.numeric(t <= 18) + as.numeric(t > 18)}
#' 
#' #Drop-off distribution: 3\% drop-off every year.
#' G = function(t){1-exp(-0.03/12*t)}
#' 
#' #Plot correlation over time after enrollment complete
#' t = seq(18, 50, 1) #Analysis time, must be greater than enrollment period 18.
#' omega = matrix(NA, nrow=6, ncol=length(t))
#' 
#' for (i in 1:length(t)){
#'   #(1) Homogeneous regardless of PD-L1 status: HR = 0.65
#'   omega[1,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.65), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.65),
#'       F.entry = Fe, G.ltfu = G, strat.ana="Y")$corr
#'       
#'   omega[2,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.65), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.65),
#'       F.entry = Fe, G.ltfu = G, strat.ana="N")$corr
#'       
#'   #(2) Stronger effect in PD-L1+: HR = 0.65; PD-L1-: HR = 0.85
#'   omega[3,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.8), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.8),
#'       F.entry = Fe, G.ltfu = G, strat.ana="Y")$corr
#'       
#'   omega[4,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.8), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.8),
#'       F.entry = Fe, G.ltfu = G, strat.ana="N")$corr
#'       
#'   #(3) Weaker effect in PD-L1+: HR = 0.85; PD-L1-: HR = 0.65
#'   omega[5,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.5), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.5),
#'       F.entry = Fe, G.ltfu = G, strat.ana="Y")$corr
#'       
#'   omega[6,i]=corrTime(T = t[i], n = list(AandB = 300, AnotB=0, BnotA=450), 
#'       r = list(AandB=1/2, AnotB =0, BnotA = 1/2, A=1/2, B=1/2), 
#'       h0=list(AandB=h0, AnotB=h0, BnotA=h0), 
#'       S0=list(AandB=S0, AnotB=S0, BnotA=S0),
#'       h1=list(AandB=h0.65, AnotB=NULL, BnotA=h0.5), 
#'       S1=list(AandB=S0.65, AnotB=NULL, BnotA=S0.5),
#'       F.entry = Fe, G.ltfu = G, strat.ana="N")$corr
#' }
#' 
#' #Plot the correlations vs time
#' plot(t, omega[1,], type="n", ylim=range(omega),
#' xlab="Analysis Time", ylab="Correlation")
#' lines(t, omega[1,], lty=1, col=1, lwd=2)
#' lines(t, omega[2,], lty=2, col=2, lwd=2)
#' lines(t, omega[3,], lty=1, col=3, lwd=2)
#' lines(t, omega[4,], lty=2, col=3, lwd=2)
#' lines(t, omega[5,], lty=1, col=4, lwd=2)
#' lines(t, omega[6,], lty=2, col=4, lwd=2)
#' legend(18,0.642, c("PD-L1+/- HR 0.65/0.65: S", 
#' "PD-L1+/- HR 0.65/0.65: U", "PD-L1+/- HR 0.65/0.8: S"),
#' lwd=rep(2,3), col=c(1,2,3), lty=c(1,2,1), bty="n", cex=0.6) 
#' 
#' legend(18,0.632, c("PD-L1+/- HR 0.65/0.8: U","PD-L1+/- HR 0.65/0.5: S","PD-L1+/- HR 0.65/0.5: U"),
#' lwd=rep(2,3), col=c(3,4,4), lty=c(2,1,2), bty="n", cex=0.6) 
#'
#' @export
#' 
corrTime = function(T = 24, n = list(AandB = 300, AnotB=0, BnotA=450), 
    r = list(AandB=1/2, AnotB =0, BnotA = 1/2), rA=1/2, rB=1/2, 
    h0=list(AandB=function(t){log(2)/12}, AnotB=function(t){log(2)/12}, 
            BnotA=function(t){log(2)/12}), 
    S0=list(AandB=function(t){exp(-log(2)/12*t)}, AnotB=function(t){exp(-log(2)/12*t)},
            BnotA=function(t){exp(-log(2)/12*t)}),
    h1=list(AandB=function(t){log(2)/12*0.70},    AnotB=function(t){log(2)/12*0.70},
            BnotA=function(t){log(2)/12*0.70}), 
    S1=list(AandB=function(t){exp(-log(2)/12 * 0.7 * t)},AnotB=function(t){exp(-log(2)/12 * 0.7 * t)},
            BnotA=function(t){exp(-log(2)/12 * 0.7 * t)}),
    F.entry = function(t){(t/18)^1.5*as.numeric(t <= 18) + as.numeric(t > 18)}, 
    G.ltfu = function(t){0}, strat.ana="Y"){

  #Expected events at time T
  e = corrEvents(T = T, n = n, r = r,h0=h0,S0=S0,h1=h1,S1=S1, F.entry = F.entry, G.ltfu = G.ltfu)
  eAandB = e$n.events.total[e$subgroup=="AandB"]
  eAnotB = e$n.events.total[e$subgroup=="AnotB"]
  eBnotA = e$n.events.total[e$subgroup=="BnotA"]
  g = n$AandB / (n$AandB+n$AnotB+n$BnotA)
    
  #Call corrBounds function for convenience to calculate the correlation
  #alpha and w are actually not needed here. Just set some values.
  bd=corrBounds(sf=list(sfuA=gsDesign::sfLDOF, sfuB=gsDesign::sfLDOF), 
                eAandB = eAandB, eAnotB = eAnotB, eBnotA = eBnotA,
                r = r, rA=rA, rB=rB, gamma = g,
                strat.ana=strat.ana, alpha=0.025, w=c(1/2,1/2), 
                method=c("Balanced Allocation"))
  corr = bd$corr[1,2]

  o=list()
  o$n.events = e
  o$corr = corr
  
  return(o)
}
phe2189/corrTests documentation built on Oct. 7, 2022, 11:13 a.m.