R/mHR2.R

Defines functions CI.boot.one do.F mHR2.LF mHR2

Documented in mHR2 mHR2.LF

#' Cox regression for a bivariate outcome
#'
#' Fits a semiparametric Cox regression model for a bivariate
#' outcome. This function computes the regression coefficients,
#' baseline hazards, and sandwich estimates of the standard
#' deviation of the regression coefficients. If desired, estimates
#' of the survival function F and marginal hazard rates Lambda11
#' can be computed using the mHR2.LF function.
#'
#' @param Y1,Y2 Vectors of event times (continuous).
#' @param Delta1,Delta2 Vectors of censoring indicators (1=event,
#' 0=censored).
#' @param X Matrix of covariates (continuous or binary).
#' @return A list containing the following elements:
#' \describe{
#' \item{Y1, Y2:}{Original vectors of event times}
#' \item{Delta1, Delta2:}{Original vectors of censoring indicators}
#' \item{X:}{Original covariate matrix}
#' \item{n10, n01:}{Total number of events for the first/second outcome}
#' \item{n11:}{Total number of double events}
#' \item{beta10, beta01, beta11:}{Regression coefficient estimates}
#' \item{lambda10, lambda01, lambda11:}{Baseline hazard estimates}
#' \item{SD.beta10, SD.beta01, SD.beta11:}{Sandwich estimates of the
#' standard deviation of the regression coefficients}
#' \item{SD.beta10.cox, SD.beta01.cox:}{Standard deviation estimates
#' for the regression coefficients based on a univariate Cox model}
#' }
#' @seealso \code{\link{mHR2.LF}}
#' @references
#' Prentice, R., Zhao, S. "The statistical analysis of multivariate
#' failure time data: A marginal modeling approach", CRC Press (2019).
#' Prentice, R., Zhao, S. "Regression models and multivariate life tables",
#' Journal of the American Statistical Association (2020) In press.
#' @useDynLib mhazard
#' @importFrom Rcpp sourceCpp
#' @importFrom survival coxph basehaz
#' @importFrom rootSolve multiroot
#' @export
#' @examples
#' x <- genClaytonReg(1000, 2, 0.5, 1, 1, log(2), log(2), log(8/3), 2, 2)
#' x.mHR2 <- mHR2(x$Y1, x$Y2, x$Delta1, x$Delta2, x$X)
mHR2 <- function(Y1, Y2, Delta1, Delta2, X){
          X <- as.matrix(X)
          n <- length(Y1)
          p <- dim(X)[2]

          T1 <- unique(sort(Y1[Delta1==1]))
          T2 <- unique(sort(Y2[Delta2==1]))
          I <- length(T1)
          J <- length(T2)

          #####counts#####
          dij_11 <- matrix(rep(0, (I+1)*(J+1)), ncol=J+1)
          dij_10 <- rep(0, I)
          dij_01 <- rep(0, J)

          id11 <- which(Delta1*Delta2==1)
          for (k in id11){
               Ik <- sum(T1<=Y1[k])
               Jk <- sum(T2<=Y2[k])
               dij_11[Ik+1, Jk+1] <- dij_11[Ik+1, Jk+1]+1
          }

          for (i in 1:I){
               dij_10[i] <- sum(Y1==T1[i] & Delta1==1)
          }
          for (j in 1:J){
               dij_01[j] <- sum(Y2==T2[j] & Delta2==1)
          }

          #####Parameter Estimates#####
          mod10 <- survival::coxph(survival::Surv(Y1, Delta1)~X, ties="breslow")
          beta10 <- mod10$coef
          mod01 <- survival::coxph(survival::Surv(Y2, Delta2)~X, ties="breslow")
          beta01 <- mod01$coef

          EE11 <- function(beta11, Y1, Y2, T1, T2, Delta1, Delta2,
                                X, dij_11){
                   return(calcEE11(beta11, Y1, Y2, T1, T2, Delta1, Delta2,
                                   X, dij_11))
          }
          beta11 <- rootSolve::multiroot(EE11, rep(0, p), Y1=Y1, Y2=Y2,
                                         T1=T1, T2=T2, Delta1=Delta1,
                                         Delta2=Delta2, X=X,
                                         dij_11=dij_11)$root


          #####variance estimates#####
          U10.k <- function(beta10, Y1, Delta1, X){
                   temp3 <- rep(NA, I)
                   X_ave <- matrix(rep(NA, p*I), nrow=p)
                   for (i in 1:I){
                        k.R <- which(Y1>=T1[i])
                        X.R <- X[k.R, , drop=FALSE]
                        temp1 <- as.vector(exp(X.R%*%beta10))
                        temp2 <- X.R*temp1
                        X_ave[, i] <- colSums(temp2)/sum(temp1)
                        temp3[i] <- dij_10[i]/sum(temp1)
                   }

                   U10.temp <- matrix(rep(0, p*n), nrow=p)
                   for (k in 1:n){
                        Ik <- sum(T1<=Y1[k])

                        temp4 <- rep(0, p)
                        if (Delta1[k]>0){
                            temp4 <- X[k, ]-X_ave[, Ik]
                        }

                        temp5 <- rep(0, p)
                        if (Ik>0){
                            temp5 <- rowSums((X[k, ]-X_ave[, 1:Ik, drop=FALSE])*temp3[1:Ik])
                        }

                        U10.temp[, k] <- temp4-temp5*as.numeric(exp(X[k, ]%*%beta10))
                   }
                   return(U10.temp)
          }

          U01.k <- function(beta01, Y2, Delta2, X){
                   temp3 <- rep(NA, J)
                   X_ave <- matrix(rep(NA, p*J), nrow=p)
                   for (j in 1:J){
                        k.R <- which(Y2>=T2[j])
                        X.R <- X[k.R, , drop=FALSE]
                        temp1 <- as.vector(exp(X.R%*%beta01))
                        temp2 <- X.R*temp1
                        X_ave[, j] <- colSums(temp2)/sum(temp1)
                        temp3[j] <- dij_01[j]/sum(temp1)
                   }

                   U01.temp <- matrix(rep(0, p*n), nrow=p)
                   for (k in 1:n){
                        Jk <- sum(T2<=Y2[k])

                        temp4 <- rep(0, p)
                        if (Delta2[k]>0){
                            temp4 <- X[k, ]-X_ave[, Jk]
                        }

                        temp5 <- rep(0, p)
                        if (Jk>0){
                            temp5 <- rowSums((X[k, ]-X_ave[, 1:Jk, drop=FALSE])*temp3[1:Jk])
                        }

                        U01.temp[, k] <- temp4-temp5*as.numeric(exp(X[k, ]%*%beta01))
                   }
                   return(U01.temp)
          }

          U11.k <- function(beta11, Y1, Y2, Delta1, Delta2, X){
                   temp3 <- matrix(rep(0, I*J), ncol=J)
                   X_ave <- Xtemp3 <- array(rep(0, p*I*J), dim=c(p, I, J))
                   for (i in 1:I){
                        for (j in 1:J){
                             if (dij_11[i+1, j+1]>0){
                                 k.R <- which(Y1>=T1[i] & Y2>=T2[j])
                                 X.R <- X[k.R, , drop=FALSE]
                                 temp1 <- as.vector(exp(X.R%*%beta11))
                                 temp2 <- X.R*temp1
                                 X_ave[, i, j] <- colSums(temp2)/sum(temp1)
                                 temp3[i, j] <- dij_11[i+1, j+1]/sum(temp1)
                                 Xtemp3[, i, j] <- X_ave[, i, j]*temp3[i, j]
                             }
                        }
                   }

                   U11.temp <- matrix(rep(0, p*n), nrow=p)
                   for (k in 1:n){
                        Ik <- sum(T1<=Y1[k])
                        Jk <- sum(T2<=Y2[k])

                        temp4 <- rep(0, p)
                        if (Delta1[k]*Delta2[k]>0){
                            temp4 <- X[k, ]-X_ave[, Ik, Jk]
                        }

                        temp5 <- temp6 <- rep(0, p)
                        if (Ik>0 & Jk>0){
                            temp5 <- X[k, ]*sum(temp3[1:Ik, 1:Jk])
                            temp6 <- apply(Xtemp3[, 1:Ik, 1:Jk, drop=FALSE], 1, sum)
                        }

                        U11.temp[, k] <- temp4-(temp5-temp6)*as.numeric(exp(X[k, ]%*%beta11))
                   }

                   return(U11.temp)
          }

          ##########Sandwich SD##########
          I10.temp <- array(rep(NA, p*p*I), dim=c(p, p, I))
          for (i in 1:I){
               t1 <- T1[i]
               k.R <- which(Y1>=t1)
               X.R <- X[k.R, , drop=FALSE]

               a <- as.vector(exp(X.R%*%beta10))
               b <- X.R*a
               temp0 <- sum(a)
               temp1 <- colSums(b)
               temp2 <- t(X.R)%*%b
               I10.temp[, , i] <- dij_10[i]*(temp0*temp2-outer(temp1, temp1))/(temp0^2)
          }
          I10 <- apply(I10.temp, c(1, 2), sum)

          I01.temp <- array(rep(NA, p*p*J), dim=c(p, p, J))
          for (j in 1:J){
               t2 <- T2[j]
               k.R <- which(Y2>=t2)
               X.R <- X[k.R, , drop=FALSE]

               a <- as.vector(exp(X.R%*%beta01))
               b <- X.R*a
               temp0 <- sum(a)
               temp1 <- colSums(b)
               temp2 <- t(X.R)%*%b
               I01.temp[, , j] <- dij_01[j]*(temp0*temp2-outer(temp1, temp1))/(temp0^2)
          }
          I01 <- apply(I01.temp, c(1, 2), sum)

          I11.temp <- array(rep(0, p*p*I*J), dim=c(p, p, I, J))
          for (i in 1:I){
               for (j in 1:J){
                    if (dij_11[i+1, j+1]>0){
                        t1 <- T1[i]
                        t2 <- T2[j]
                        k.R <- which(Y1>=t1 & Y2>=t2)

                        if (length(k.R)>0){
                            X.R <- X[k.R, , drop=FALSE]

                            a <- as.vector(exp(X.R%*%beta11))
                            b <- X.R*a
                            temp0 <- sum(a)
                            temp1 <- colSums(b)
                            temp2 <- t(X.R)%*%b
                            I11.temp[, , i, j] <- dij_11[i+1, j+1]*(temp0*temp2-outer(temp1, temp1))/(temp0^2)
                        }
                    }
               }
          }
          I11 <- apply(I11.temp, c(1, 2), sum)

          I_hat <- matrix(rep(0, 9*p*p), ncol=3*p)
          I_hat[1:p, 1:p] <- I10
          I_hat[(p+1):(2*p), (p+1):(2*p)] <- I01
          I_hat[(2*p+1):(3*p), (2*p+1):(3*p)] <- I11
          U_hat <- rbind(U10.k(beta10, Y1, Delta1, X), U01.k(beta01, Y2, Delta2, X), U11.k(beta11, Y1, Y2, Delta1, Delta2, X))
          A <- U_hat%*%t(U_hat)
          SD.sand <- sqrt(diag(solve(I_hat)%*%A%*%solve(I_hat)))

          ##########Baseline hazards##########
          cumhaz10 <- unique(survival::basehaz(mod10, centered=FALSE)$hazard)
          cumhaz01 <- unique(survival::basehaz(mod01, centered=FALSE)$hazard)
          lambda10 <- cumhaz10[2:(I+1)]-cumhaz10[1:I]
          lambda01 <- cumhaz01[2:(J+1)]-cumhaz01[1:J]
          lambda11 <- calc_lambda11(T1, T2, Y1, Y2, X, dij_11, beta11)

          return(list(Y1=Y1, Y2=Y2, Delta1=Delta1, Delta2=Delta2, X=X, n10=sum(Delta1), n01=sum(Delta2), n11=sum(Delta1*Delta2), beta10=beta10, beta01=beta01, beta11=beta11, SD.beta10=SD.sand[1:p], SD.beta01=SD.sand[(p+1):(2*p)], SD.beta11=SD.sand[(2*p+1):(3*p)], SD.beta10.cox=sqrt(mod10$var), SD.beta10.cox=sqrt(mod01$var),lambda10=lambda10, lambda01=lambda01, lambda11=lambda11))
}

#' Bivariate regression survival function and marginal hazards estimation
#'
#' Estimates the survival function F and the marginal hazards Lambda11
#' for a bivariate Cox regression model. F and Lambda11 are estimated
#' at two specified values of the covariates. If desired, (bootstrap)
#' confidence intervals or confidence bounds for F and Lambda11 may also
#' be computed.
#'
#' @param mHR2.obj Output from the mHR2 function.
#' @param X0_out,X1_out Two possible sets of values for the covariates.
#' F and Lambda will be estimated at X=X0_out and X=X1_out.
#' @param T1_out,T2_out Vector of time points at which F and Lambda11
#' should be estimated. If confidence="CB", then both vectors must
#' have length 3.
#' @param confidence Type of confidence estimate to be computed.
#' Possible values include "none", "CI" (to compute confidence
#' intervals), and "CB" (to compute confidence bands). Defaults to
#' "none".
#' @param n.boot Number of bootstrap iterations for computing the
#' confidence intervals/bands. Defaults to 100. Ignored if
#' confidence="none".
#' @section Details:
#' If confidence="CI" or confidence="CB", then 95% bootstrap confidence
#' bounds are computed by estimating the standard errors of F/Lambda11
#' based on n.boot bootstrap iterations. Currently confidence bounds
#' can only be computed at three specified T1out/T2out combinations
#' (meaning that T1out and T2out must both have length 3 if
#' confidence="CB"). No confidence measures will be returned if
#' confidence="none".
#' @return A list containing the following elements:
#' \describe{
#' \item{n10, n01:}{Total number of events for the first/second outcome}
#' \item{n11:}{Total number of double events}
#' \item{beta10, beta01, beta11:}{Regression coefficient estimates}
#' \item{lambda10, lambda01, lambda11:}{Baseline hazard estimates}
#' \item{Lambda11_out_Z0, Lambda11_out_Z1:}{Estimates of Lambda11 at
#' T1_out, T2_out for X=X0_out and X=X1_out}
#' \item{F_out_X0, F_out_X1:}{Estimates of F at T1_out, T2_out for
#' X=X0_out and X=X1_out}
#' \item{CI_Lambda11_X0.lb, CI_Lambda11_X0.ub:}{Lower and upper bounds
#' for Lambda11 at X=X0_out}
#' \item{CI_Lambda11_X1.lb, CI_Lambda11_X1.ub:}{Lower and upper bounds
#' for Lambda11 at X=X1_out}
#' \item{CI_F_X0.lb, CI_F_X0.ub:}{Lower and upper bounds for F at
#' X=X0_out}
#' \item{CI_F_X1.lb, CI_F_X1.ub:}{Lower and upper bounds for F at
#' X=X1_out}
#' \item{CB1_Lambda11_X0.lb, CB1_Lambda11_X0.ub, CB2_Lambda11_X0.lb,
#' CB2_Lambda11_X0.ub, CB3_Lambda11_X0.lb, CB3_Lambda11_X0.ub:}{Lower
#' and upper bounds for Lambda11 at X=X0_out, at three T1_out, T2_out
#' combinations}
#' \item{CB1_Lambda11_X1.lb, CB1_Lambda11_X1.ub, CB2_Lambda11_X1.lb,
#' CB2_Lambda11_X1.ub, CB3_Lambda11_X1.lb, CB3_Lambda11_X1.ub:}{Lower
#' and upper bounds for Lambda11 at X=X1_out, at three T1_out, T2_out
#' combinations}
#' \item{CB1_F_X0.lb, CB1_F_X0.ub, CB2_F_X0.lb, CB2_F_X0.ub, CB3_F_X0.lb,
#' CB3_F_X0.ub:}{Lower and upper bounds for F at X=X0_out, at three
#' T1_out, T2_out combinations}
#' \item{CB1_F_X1.lb, CB1_F_X1.ub, CB2_F_X1.lb, CB2_F_X1.ub, CB3_F_X1.lb,
#' CB3_F_X1.ub:}{Lower and upper bounds for F at X=X1_out, at three
#' T1_out, T2_out combinations}
#' }
#' @seealso \code{\link{mHR2}}
#' @references
#' Prentice, R., Zhao, S. "The statistical analysis of multivariate
#' failure time data: A marginal modeling approach", CRC Press (2019).
#' Prentice, R., Zhao, S. "Regression models and multivariate life tables",
#' Journal of the American Statistical Association (2020) In press.
#' @useDynLib mhazard
#' @importFrom Rcpp sourceCpp
#' @importFrom survival coxph basehaz Surv
#' @importFrom rootSolve multiroot
#' @importFrom stats quantile sd
#' @export
#' @examples
#' x <- genClaytonReg(1000, 2, 0.5, 1, 1, log(2), log(2), log(8/3), 2, 2)
#' x.mHR2 <- mHR2(x$Y1, x$Y2, x$Delta1, x$Delta2, x$X)
#' x.LF <- mHR2.LF(x.mHR2, 0, 1, c(0.25, 0.5, 1), c(0.25, 0.5, 1))
#' \donttest{x.LF.CI <- mHR2.LF(x.mHR2, 0, 1, c(0.25, 0.5, 1),
#' c(0.25, 0.5, 1), confidence="CI")}
#' \donttest{x.LF.CB <- mHR2.LF(x.mHR2, 0, 1, c(0.25, 0.5, 1),
#' c(0.25, 0.5, 1), confidence="CB")}
mHR2.LF <- function(mHR2.obj, X0_out, X1_out, T1_out, T2_out,
                    confidence=c("none", "CI", "CB"), n.boot=100){
    Y1 <- mHR2.obj$Y1
    Y2 <- mHR2.obj$Y2
    Delta1 <- mHR2.obj$Delta1
    Delta2 <- mHR2.obj$Delta2
    X <- as.matrix(mHR2.obj$X)
    n.sub <- length(Y1)

    confidence <- match.arg(confidence)
    if (confidence=="CI") {
        result <- CI.boot.one(Y1, Y2, Delta1, Delta2, X, X0_out, X1_out, T1_out, T2_out)
        Lambda11.boot_Z0 <- Lambda11.boot_Z1 <- F.boot_Z0 <- F.boot_Z1 <- array(rep(NA, length(T1_out)*length(T2_out)*n.boot), dim=c(length(T1_out), length(T2_out), n.boot))

        for (i.boot in 1:n.boot){
            index <- sample(1:n.sub, replace=TRUE)
            Y1.boot <- Y1[index]
            Y2.boot <- Y2[index]
            Delta1.boot <- Delta1[index]
            Delta2.boot <- Delta2[index]
            X.boot <- as.matrix(X[index, ])
            result.boot <- CI.boot.one(Y1.boot, Y2.boot, Delta1.boot, Delta2.boot, X.boot, X0_out, X1_out, T1_out, T2_out)
            Lambda11.boot_Z0[, , i.boot] <- result.boot$Lambda11_out_Z0
            Lambda11.boot_Z1[, , i.boot] <- result.boot$Lambda11_out_Z1
            F.boot_Z0[, , i.boot] <- result.boot$F_out_Z0
            F.boot_Z1[, , i.boot] <- result.boot$F_out_Z1
        }

        SE_Lambda11_Z0 <- apply(Lambda11.boot_Z0, c(1, 2), function(x){sd(x, na.rm=TRUE)})
        SE_Lambda11_Z1 <- apply(Lambda11.boot_Z1, c(1, 2), function(x){sd(x, na.rm=TRUE)})
        SE_F_Z0 <- apply(F.boot_Z0, c(1, 2), function(x){sd(x, na.rm=TRUE)})
        SE_F_Z1 <- apply(F.boot_Z1, c(1, 2), function(x){sd(x, na.rm=TRUE)})

        CI_Lambda11_Z0.lb <- result$Lambda11_out_Z0-1.96*SE_Lambda11_Z0
        CI_Lambda11_Z0.ub <- result$Lambda11_out_Z0+1.96*SE_Lambda11_Z0
        CI_Lambda11_Z1.lb <- result$Lambda11_out_Z1-1.96*SE_Lambda11_Z1
        CI_Lambda11_Z1.ub <- result$Lambda11_out_Z1+1.96*SE_Lambda11_Z1

        CI_F_Z0.lb <- result$F_out_Z0-1.96*SE_F_Z0
        CI_F_Z0.ub <- result$F_out_Z0+1.96*SE_F_Z0
        CI_F_Z1.lb <- result$F_out_Z1-1.96*SE_F_Z1
        CI_F_Z1.ub <- result$F_out_Z1+1.96*SE_F_Z1

        return(list(n10=result$n10, n01=result$n01, n11=result$n11,
                    beta10=result$beta10, beta01=result$beta01, beta11=result$beta11,
                    lambda10=result$lambda10, lambda01=result$lambda01, lambda11=result$lambda11,
                    Lambda11_out_X0=result$Lambda11_out_Z0, Lambda11_out_X1=result$Lambda11_out_Z1,
                    CI_Lambda11_X0.lb=CI_Lambda11_Z0.lb, CI_Lambda11_X0.ub=CI_Lambda11_Z0.ub,
                    CI_Lambda11_X1.lb=CI_Lambda11_Z1.lb, CI_Lambda11_X1.ub=CI_Lambda11_Z1.ub,
                    F_out_X0=result$F_out_Z0, F_out_X1=result$F_out_Z1,
                    CI_F_X0.lb=CI_F_Z0.lb, CI_F_X0.ub=CI_F_Z0.ub,
                    CI_F_X1.lb=CI_F_Z1.lb, CI_F_X1.ub=CI_F_Z1.ub))
    }
    else if (confidence=="CB") {
        if (length(T1_out)!=3 | length(T2_out)!=3) {
            stop("T1_out, T2_out must have length 3 for confidence='CB'")
        }
        T1 <- sort(unique(Y1[Delta1==1]))
        T2 <- sort(unique(Y2[Delta2==1]))
        result <- CI.boot.one(Y1, Y2, Delta1, Delta2, X, X0_out, X1_out, T1, T2)

        Lambda11_Z0_hat <- result$Lambda11_out_Z0
        Lambda11_Z1_hat <- result$Lambda11_out_Z1
        F_Z0_hat <- result$F_out_Z0
        F_Z1_hat <- result$F_out_Z1

        ind.T1 <- unlist(lapply(T1_out, function(x){sum(x>=T1)}))
        ind.T2 <- unlist(lapply(T2_out, function(x){sum(x>=T2)}))

        mW_Lambda11.boot_Z0 <- mW_Lambda11.boot_Z1 <- mW_F.boot_Z0 <- mW_F.boot_Z1 <- matrix(rep(NA, 3*n.boot), ncol=3)

        for (i.boot in 1:n.boot){
            index <- sample(1:n.sub, replace=TRUE)
            Y1.boot <- Y1[index]
            Y2.boot <- Y2[index]
            Delta1.boot <- Delta1[index]
            Delta2.boot <- Delta2[index]
            X.boot <- as.matrix(X[index, ])
            result.boot <- CI.boot.one(Y1.boot, Y2.boot, Delta1.boot, Delta2.boot, X.boot, X0_out, X1_out, T1, T2, calc_rij=TRUE)


            W_Lambda11.boot_Z0 <- sqrt(n.sub)*(result.boot$Lambda11_out_Z0-Lambda11_Z0_hat)*I(result.boot$rij_out>0)
            W_Lambda11.boot_Z1 <- sqrt(n.sub)*(result.boot$Lambda11_out_Z1-Lambda11_Z1_hat)*I(result.boot$rij_out>0)
            W_F.boot_Z0 <- sqrt(n.sub)*(result.boot$F_out_Z0-F_Z0_hat)*I(result.boot$rij_out>0)
            W_F.boot_Z1 <- sqrt(n.sub)*(result.boot$F_out_Z1-F_Z1_hat)*I(result.boot$rij_out>0)

            r1.boot.Z0 <- sum((Y1.boot[X.boot==0]>=T1_out[1]) & (Y2.boot[X.boot==0]>T2_out[1]))
            r2.boot.Z0 <- sum((Y1.boot[X.boot==0]>=T1_out[2]) & (Y2.boot[X.boot==0]>T2_out[2]))
            r3.boot.Z0 <- sum((Y1.boot[X.boot==0]>=T1_out[3]) & (Y2.boot[X.boot==0]>T2_out[3]))
            r1.boot.Z1 <- sum((Y1.boot[X.boot==1]>=T1_out[1]) & (Y2.boot[X.boot==1]>T2_out[1]))
            r2.boot.Z1 <- sum((Y1.boot[X.boot==1]>=T1_out[2]) & (Y2.boot[X.boot==1]>T2_out[2]))
            r3.boot.Z1 <- sum((Y1.boot[X.boot==1]>=T1_out[3]) & (Y2.boot[X.boot==1]>T2_out[3]))

            if (r1.boot.Z0>0){
                mW_Lambda11.boot_Z0[i.boot, 1] <- max(abs(W_Lambda11.boot_Z0[1:ind.T1[1], 1:ind.T2[1]]))
                mW_F.boot_Z0[i.boot, 1] <- max(abs(W_F.boot_Z0[1:ind.T1[1], 1:ind.T2[1]]))
            }
            if (r2.boot.Z0>0){
                mW_Lambda11.boot_Z0[i.boot, 2] <- max(abs(W_Lambda11.boot_Z0[1:ind.T1[2], 1:ind.T2[2]]))
                mW_F.boot_Z0[i.boot, 2] <- max(abs(W_F.boot_Z0[1:ind.T1[2], 1:ind.T2[2]]))
            }
            if (r3.boot.Z0>0){
                mW_Lambda11.boot_Z0[i.boot, 3] <- max(abs(W_Lambda11.boot_Z0[1:ind.T1[3], 1:ind.T2[3]]))
                mW_F.boot_Z0[i.boot, 3] <- max(abs(W_F.boot_Z0[1:ind.T1[3], 1:ind.T2[3]]))
            }

            if (r1.boot.Z1>0){
                mW_Lambda11.boot_Z1[i.boot, 1] <- max(abs(W_Lambda11.boot_Z1[1:ind.T1[1], 1:ind.T2[1]]))
                mW_F.boot_Z1[i.boot, 1] <- max(abs(W_F.boot_Z1[1:ind.T1[1], 1:ind.T2[1]]))
            }
            if (r2.boot.Z1>0){
                mW_Lambda11.boot_Z1[i.boot, 2] <- max(abs(W_Lambda11.boot_Z1[1:ind.T1[2], 1:ind.T2[2]]))
                mW_F.boot_Z1[i.boot, 2] <- max(abs(W_F.boot_Z1[1:ind.T1[2], 1:ind.T2[2]]))
            }
            if (r3.boot.Z1>0){
                mW_Lambda11.boot_Z1[i.boot, 3] <- max(abs(W_Lambda11.boot_Z1[1:ind.T1[3], 1:ind.T2[3]]))
                mW_F.boot_Z1[i.boot, 3] <- max(abs(W_F.boot_Z1[1:ind.T1[3], 1:ind.T2[3]]))
            }
        }

        C.Lambda11_Z0 <- apply(mW_Lambda11.boot_Z0, 2, function(x){quantile(x, 0.95, na.rm=TRUE)})
        C.Lambda11_Z1 <- apply(mW_Lambda11.boot_Z1, 2, function(x){quantile(x, 0.95, na.rm=TRUE)})
        C.F_Z0 <- apply(mW_F.boot_Z0, 2, function(x){quantile(x, 0.95, na.rm=TRUE)})
        C.F_Z1 <- apply(mW_F.boot_Z1, 2, function(x){quantile(x, 0.95, na.rm=TRUE)})

        CB1_Lambda11_Z0.lb <- Lambda11_Z0_hat[1:ind.T1[1], 1:ind.T2[1]]-C.Lambda11_Z0[1]/sqrt(n.sub)
        CB1_Lambda11_Z0.ub <- Lambda11_Z0_hat[1:ind.T1[1], 1:ind.T2[1]]+C.Lambda11_Z0[1]/sqrt(n.sub)
        CB1_Lambda11_Z1.lb <- Lambda11_Z1_hat[1:ind.T1[1], 1:ind.T2[1]]-C.Lambda11_Z1[1]/sqrt(n.sub)
        CB1_Lambda11_Z1.ub <- Lambda11_Z1_hat[1:ind.T1[1], 1:ind.T2[1]]+C.Lambda11_Z1[1]/sqrt(n.sub)

        CB1_F_Z0.lb <- F_Z0_hat[1:ind.T1[1], 1:ind.T2[1]]-C.F_Z0[1]/sqrt(n.sub)
        CB1_F_Z0.ub <- F_Z0_hat[1:ind.T1[1], 1:ind.T2[1]]+C.F_Z0[1]/sqrt(n.sub)
        CB1_F_Z1.lb <- F_Z1_hat[1:ind.T1[1], 1:ind.T2[1]]-C.F_Z1[1]/sqrt(n.sub)
        CB1_F_Z1.ub <- F_Z1_hat[1:ind.T1[1], 1:ind.T2[1]]+C.F_Z1[1]/sqrt(n.sub)

        CB2_Lambda11_Z0.lb <- Lambda11_Z0_hat[1:ind.T1[2], 1:ind.T2[2]]-C.Lambda11_Z0[2]/sqrt(n.sub)
        CB2_Lambda11_Z0.ub <- Lambda11_Z0_hat[1:ind.T1[2], 1:ind.T2[2]]+C.Lambda11_Z0[2]/sqrt(n.sub)
        CB2_Lambda11_Z1.lb <- Lambda11_Z1_hat[1:ind.T1[2], 1:ind.T2[2]]-C.Lambda11_Z1[2]/sqrt(n.sub)
        CB2_Lambda11_Z1.ub <- Lambda11_Z1_hat[1:ind.T1[2], 1:ind.T2[2]]+C.Lambda11_Z1[2]/sqrt(n.sub)

        CB2_F_Z0.lb <- F_Z0_hat[1:ind.T1[2], 1:ind.T2[2]]-C.F_Z0[2]/sqrt(n.sub)
        CB2_F_Z0.ub <- F_Z0_hat[1:ind.T1[2], 1:ind.T2[2]]+C.F_Z0[2]/sqrt(n.sub)
        CB2_F_Z1.lb <- F_Z1_hat[1:ind.T1[2], 1:ind.T2[2]]-C.F_Z1[2]/sqrt(n.sub)
        CB2_F_Z1.ub <- F_Z1_hat[1:ind.T1[2], 1:ind.T2[2]]+C.F_Z1[2]/sqrt(n.sub)

        CB3_Lambda11_Z0.lb <- Lambda11_Z0_hat[1:ind.T1[3], 1:ind.T2[3]]-C.Lambda11_Z0[3]/sqrt(n.sub)
        CB3_Lambda11_Z0.ub <- Lambda11_Z0_hat[1:ind.T1[3], 1:ind.T2[3]]+C.Lambda11_Z0[3]/sqrt(n.sub)
        CB3_Lambda11_Z1.lb <- Lambda11_Z1_hat[1:ind.T1[3], 1:ind.T2[3]]-C.Lambda11_Z1[3]/sqrt(n.sub)
        CB3_Lambda11_Z1.ub <- Lambda11_Z1_hat[1:ind.T1[3], 1:ind.T2[3]]+C.Lambda11_Z1[3]/sqrt(n.sub)

        CB3_F_Z0.lb <- F_Z0_hat[1:ind.T1[3], 1:ind.T2[3]]-C.F_Z0[3]/sqrt(n.sub)
        CB3_F_Z0.ub <- F_Z0_hat[1:ind.T1[3], 1:ind.T2[3]]+C.F_Z0[3]/sqrt(n.sub)
        CB3_F_Z1.lb <- F_Z1_hat[1:ind.T1[3], 1:ind.T2[3]]-C.F_Z1[3]/sqrt(n.sub)
        CB3_F_Z1.ub <- F_Z1_hat[1:ind.T1[3], 1:ind.T2[3]]+C.F_Z1[3]/sqrt(n.sub)

        return(list(n=n.sub, n10=result$n10, n01=result$n01, n11=result$n11, T1=T2, T2=T2,
                    beta10=result$beta10, beta01=result$beta01, beta11=result$beta11, lambda10=result$lambda10, lambda01=result$lambda01, lambda11=result$lambda11,
                    Lambda11_out_X0=result$Lambda11_out_Z0, Lambda11_out_X1=result$Lambda11_out_Z1,
                    CB1_Lambda11_X0.lb=CB1_Lambda11_Z0.lb, CB1_Lambda11_X0.ub=CB1_Lambda11_Z0.ub,
                    CB1_Lambda11_X1.lb=CB1_Lambda11_Z1.lb, CB1_Lambda11_X1.ub=CB1_Lambda11_Z1.ub,
                    CB2_Lambda11_X0.lb=CB2_Lambda11_Z0.lb, CB2_Lambda11_X0.ub=CB2_Lambda11_Z0.ub,
                    CB2_Lambda11_X1.lb=CB2_Lambda11_Z1.lb, CB2_Lambda11_X1.ub=CB2_Lambda11_Z1.ub,
                    CB3_Lambda11_X0.lb=CB3_Lambda11_Z0.lb, CB3_Lambda11_X0.ub=CB3_Lambda11_Z0.ub,
                    CB3_Lambda11_X1.lb=CB3_Lambda11_Z1.lb, CB3_Lambda11_X1.ub=CB3_Lambda11_Z1.ub,
                    F_out_X0=result$F_out_Z0, F_out_X1=result$F_out_Z1,
                    CB1_F_X0.lb=CB1_F_Z0.lb, CB1_F_X0.ub=CB1_F_Z0.ub,
                    CB1_F_X1.lb=CB1_F_Z1.lb, CB1_F_X1.ub=CB1_F_Z1.ub,
                    CB2_F_X0.lb=CB2_F_Z0.lb, CB2_F_X0.ub=CB2_F_Z0.ub,
                    CB2_F_X1.lb=CB2_F_Z1.lb, CB2_F_X1.ub=CB2_F_Z1.ub,
                    CB3_F_X0.lb=CB3_F_Z0.lb, CB3_F_X0.ub=CB3_F_Z0.ub,
                    CB3_F_X1.lb=CB3_F_Z1.lb, CB3_F_X1.ub=CB3_F_Z1.ub))
    }
    else {
        result <- CI.boot.one(Y1, Y2, Delta1, Delta2, X, X0_out, X1_out, T1_out, T2_out)
        return(list(n10=result$n10, n01=result$n01, n11=result$n11,
                    beta10=result$beta10, beta01=result$beta01, beta11=result$beta11,
                    lambda10=result$lambda10, lambda01=result$lambda01, lambda11=result$lambda11,
                    Lambda11_out_X0=result$Lambda11_out_Z0, Lambda11_out_X1=result$Lambda11_out_Z1,
                    F_out_X0=result$F_out_Z0, F_out_X1=result$F_out_Z1))
    }
}

do.F <- function(T1, T2, beta10, beta01, beta11, lambda10, lambda01, lambda11, x){
        I <- length(T1)
        J <- length(T2)

        lambda10.x <- rep(NA, I+1)
        lambda01.x <- rep(NA, J+1)
        lambda11.x <- matrix(rep(NA, (I+1)*(J+1)), ncol=(J+1))
        lambda10.x[1] <- 0
        lambda01.x[1] <- 0
        lambda10.x[2:(I+1)] <- lambda10*exp(as.numeric(x%*%beta10))
        lambda01.x[2:(J+1)] <- lambda01*exp(as.numeric(x%*%beta01))
        lambda11.x <- lambda11*exp(as.numeric(x%*%beta11))

        F.x <- calcF(lambda10.x, lambda01.x, lambda11.x)
        return(F.x)
}

CI.boot.one <- function(Y1, Y2, Delta1, Delta2, X, X0_out, X1_out, T1_out, T2_out, calc_rij=FALSE){
               X <- as.matrix(X)
               n <- length(Y1)
               p <- dim(X)[2]

               T1 <- unique(sort(Y1[Delta1==1]))
               T2 <- unique(sort(Y2[Delta2==1]))
               I <- length(T1)
               J <- length(T2)

               #####counts#####
               if (calc_rij) {
                   junk <- calc_drij(Y1, Y2, c(0, T1), c(0, T2), Delta1, Delta2)
                   dij_11 <- junk[,,2]
                   rij <- junk[,,1]
               }
               else {
                   dij_11 <- matrix(rep(0, (I+1)*(J+1)), ncol=J+1)
                   id11 <- which(Delta1*Delta2==1)
                   for (k in id11){
                       Ik <- sum(T1<=Y1[k])
                       Jk <- sum(T2<=Y2[k])
                       dij_11[Ik+1, Jk+1] <- dij_11[Ik+1, Jk+1]+1
                   }
               }

               #####Parameter Estimates#####
               mod10 <- survival::coxph(survival::Surv(Y1, Delta1)~X, ties="breslow")
               beta10 <- mod10$coef
               mod01 <- survival::coxph(survival::Surv(Y2, Delta2)~X, ties="breslow")
               beta01 <- mod01$coef

               EE11 <- function(beta11, Y1, Y2, T1, T2, Delta1, Delta2,
                                X, dij_11){
                   return(calcEE11(beta11, Y1, Y2, T1, T2, Delta1, Delta2,
                                   X, dij_11))
               }
               beta11 <- rootSolve::multiroot(EE11, rep(0, p), Y1=Y1, Y2=Y2,
                                              T1=T1, T2=T2, Delta1=Delta1,
                                              Delta2=Delta2, X=X,
                                              dij_11=dij_11)$root

               cumhaz10 <- unique(survival::basehaz(mod10,
                                                    centered=FALSE)$hazard)
               cumhaz01 <- unique(survival::basehaz(mod01,
                                                    centered=FALSE)$hazard)
               lambda10 <- cumhaz10[2:(I+1)]-cumhaz10[1:I]
               lambda01 <- cumhaz01[2:(J+1)]-cumhaz01[1:J]
               lambda11 <- calc_lambda11(T1, T2, Y1, Y2, X, dij_11, beta11)

               lambda11_Z0 <- lambda11*exp(as.numeric(X0_out%*%beta11))
               lambda10_Z0 <- lambda10*exp(as.numeric(X0_out%*%beta10))
               lambda01_Z0 <- lambda01*exp(as.numeric(X0_out%*%beta01))

               lambda11_Z1 <- lambda11*exp(as.numeric(X1_out%*%beta11))
               lambda10_Z1 <- lambda10*exp(as.numeric(X1_out%*%beta10))
               lambda01_Z1 <- lambda01*exp(as.numeric(X1_out%*%beta01))

               Lambda11_Z0 <- t(apply(apply(lambda11_Z0, 2, cumsum), 1, cumsum))
               Lambda11_Z1 <- t(apply(apply(lambda11_Z1, 2, cumsum), 1, cumsum))

               F_Z0 <- do.F(T1, T2, beta10, beta01, beta11, lambda10, lambda01, lambda11, X0_out)
               F_Z1 <- do.F(T1, T2, beta10, beta01, beta11, lambda10, lambda01, lambda11, X1_out)

               index.i <- unlist(lapply(T1_out, function(x){sum(T1<=x)}))
               index.j <- unlist(lapply(T2_out, function(x){sum(T2<=x)}))

               if (calc_rij) {
                   rij_out <- rij[index.i+1, index.j+1]
               }
               Lambda11_out_Z0 <- Lambda11_Z0[index.i+1, index.j+1]
               Lambda11_out_Z1 <- Lambda11_Z1[index.i+1, index.j+1]
               F_out_Z0 <- F_Z0[index.i+1, index.j+1]
               F_out_Z1 <- F_Z1[index.i+1, index.j+1]

               if (calc_rij) {
                   return(list(n10=sum(Delta1), n01=sum(Delta2), n11=sum(Delta1*Delta2),
                               T1=T1, T2=T2,
                               beta10=beta10, beta01=beta01, beta11=beta11,
                               lambda10=lambda10, lambda01=lambda01, lambda11=lambda11,
                               rij_out=rij_out,
                               Lambda11_out_Z0=Lambda11_out_Z0, Lambda11_out_Z1=Lambda11_out_Z1,
                               F_out_Z0=F_out_Z0, F_out_Z1=F_out_Z1))
               }
               else {
                   return(list(n10=sum(Delta1), n01=sum(Delta2), n11=sum(Delta1*Delta2),
                               T1=T1, T2=T2,
                               beta10=beta10, beta01=beta01, beta11=beta11,
                               lambda10=lambda10, lambda01=lambda01, lambda11=lambda11,
                               Lambda11_out_Z0=Lambda11_out_Z0, Lambda11_out_Z1=Lambda11_out_Z1,
                               F_out_Z0=F_out_Z0, F_out_Z1=F_out_Z1))
               }
}

Try the mhazard package in your browser

Any scripts or data that you put into this service are public.

mhazard documentation built on May 30, 2022, 1:07 a.m.