#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.