R/truncAIPW_c2.R

In truncAIPW: Doubly Robust Estimation under Covariate-Induced Dependent Left Truncation

### function that compute the truncAIPW estimator under "c2" type of right censoring - censoring always after left truncation

#' @title Doubly Robust Estimation under Covariate-induced Dependent Left Truncation and Noninformative Right Censoring where Censoring is always after Left Truncation
#' @description Doubly robust estimation of the mean of an arbitrarily transformed survival time under covariate-induced dependent left truncation and noninformative right censoring where censoring is always after left truncation. Inverse probability of censoring weighting is used to handle the right censoring.
#' @references Wang, Y., Ying, A., Xu, R. (2022) "Doubly robust estimation under covariate-induced dependent left truncation" <arXiv:2208.06836>.
#' @param dat data frame that contains the data for constructing the estimating equation.
#' @param nu transformation that defines the parameter of interest.
#' @param Fuz.mx matrix for the estimated conditional CDF of the event time given covariates. Each row corresponds to a subject, and each column corresponds to a time point. The column names of the matrix are the time points. See \code{\link{F_est}} for an example of computing this input matrix for the conditional CDF.
#' @param Gvz.mx matrix for the estimated conditional CDF of the truncation time given covariates. Each row corresponds to a subject, and each column corresponds to a time point. The column names of the matrix are the time points. See \code{\link{G_est}} for an example of computing this input matrix for the conditional CDF.
#' @param wd vector for the inverse probability of residual censoring weights \eqn{\Delta/\hat S_D(X-Q)}.
#' @param X.name name of the censored event time variable X = min(T, C).
#' @param Q.name name of the left truncation time variable.
#' @param status.name name of the event time indicator.
#' @param trim constant that is used to bound from below for the denominators involved in the computation.
#' @return \code{truncAIPW_cen2()} returns a list of estimators (`dr', `IPW.Q', `Reg.T1', `Reg.T2').
#' \item{dr}{doubly robust estimator 'dr'.}
#' \item{IPW.Q}{inverse probability of truncation weighted estimator 'IPW.Q'.}
#' \item{Reg.T1}{regression based estimator 'Reg.T1'.}
#' \item{Reg.T2}{regression based estimator 'Reg.T2'.}
#' @export
#' @import survival stats
#' @seealso See also \code{\link{truncAIPW}} for estimation under no censoring, and \code{\link{truncAIPW_cen1}} for estimation under another type of noninformative right censoring. See also \code{\link{F_est}}, \code{\link{G_est}} as examples for computing the input matrices of the conditional CDF's.
#' @examples
#' library(survival)
#' data("simu_c2")
#' nu <- function(t){ return(as.numeric(t>3)) }
#'
#' u = c(min(simu_c2$X)-1e-10, sort(simu_c2$X), max(simu_c2$X)+1e-10)
#' v = c(min(simu_c2$Q)-1e-10, sort(simu_c2$Q), max(simu_c2$Q)+1e-10)
#'
#' kmfit.D = survfit(Surv(X-Q, 1-delta)~1, data = simu_c2, type = "kaplan-meier")
#' Sd = stepfun(kmfit.D$time,  c(1, kmfit.D$surv))
#' wd = rep(0, nrow(simu_c2))
#' wd[which(simu_c2$delta == 1)] = 1/Sd(simu_c2$X - simu_c2$Q)[which(simu_c2$delta == 1)]
#'
#' simu_c2$wd = wd
#' simu_c2.1 = simu_c2[simu_c2$delta==1,]
#' wd_1 = simu_c2.1$wd
#'
#' Fuz.mx = F_est(simu_c2, simu_c2, u, "Cox", "X", "Q", "delta", c("Z1","Z2"))
#' Gvz.mx = G_est(simu_c2.1, simu_c2, v, "Cox", "X", "Q", "delta", c("Z1","Z2"), weights = wd_1)
#'
#' est = truncAIPW_cen2(simu_c2, nu, Fuz.mx, Gvz.mx, wd, "X", "Q", "delta", trim = 1e-7)
#' est
truncAIPW_cen2 <- function(dat, nu, Fuz.mx, Gvz.mx, wd, X.name, Q.name, status.name, trim = 1e-7){

    efs = truncAIPW_transMean_EF(dat, nu, Fuz.mx, Gvz.mx, X.name, Q.name, trim)

    Num_AIPW = efs$Num_AIPW
    Den_AIPW = efs$Den_AIPW

    Num_IPW.Q = efs$Num_IPW.Q
    Den_IPW.Q = efs$Den_IPW.Q

    Num_Reg.T1 = efs$Num_Reg.T1
    Num_Reg.T2 = efs$Num_Reg.T2
    Den_Reg = efs$Den_Reg


    est_dr = mean(wd*Num_AIPW)/mean(wd*Den_AIPW)
    est_IPWQ = mean(wd*Num_IPW.Q)/mean(wd*Den_IPW.Q)
    est_RegT1 = mean(wd*Num_Reg.T1)/mean(wd*Den_Reg)
    est_RegT2 = mean(wd*Num_Reg.T2)/mean(wd*Den_Reg)


    return(list(dr = est_dr,
                IPW.Q = est_IPWQ,
                Reg.T1 = est_RegT1,
                Reg.T2 = est_RegT2))
}