R/anova.ate.R

In riskRegression: Risk Regression Models and Prediction Scores for Survival Analysis with Competing Risks

### anova.ate.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: aug 19 2020 (09:18) 
## Version: 
## Last-Updated: Sep  6 2023 (09:53) 
##           By: Thomas Alexander Gerds
##     Update #: 72
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

## * anova.ate (documentation)
#' @title Risk Comparison Over Time
#' @description Comparison of risk differences or risk ratios over all timepoints.
#'
#' @param object A \code{ate} object, i.e. output of the \code{ate} function.
#' @param allContrast [matrix] contrast for which the risks should be compared.
#' Matrix with two rows, the first being the sequence of reference treatments and the second the sequence of alternative treatments. 
#' @param type [character vector] the functionnal used to compare the risks: \code{"diffRisk"} or \code{"ratioRisk"}.
#' @param estimator [character] The type of estimator relative to which the comparison should be performed. 
#' @param test [character] The type of statistic used to compare the risks over times:
#' \code{"KM"} (extremum risk), \code{"CvM"} (sum of squares of the risk), or \code{"sum"} (sum of the risks).
#' @param transform [character] Should a transformation be used, e.g. the test is performed after log-transformation of the estimate, standard error, and influence function.
#' @param alternative [character] a character string specifying the alternative hypothesis, must be one of \code{"two.sided"}, \code{"greater"} or \code{"less"}.
#' @param n.sim [integer, >0] the number of simulations used to compute the p-values.
#' @param print [logical] should the results be displayed?
#' @param ... Not used.
#'
#' @details Experimental!!!

## * confint.ate (examples)
##' @examples
##' library(survival)
##' library(data.table)
##' 
##' \dontrun{
##' ## simulate data
##' set.seed(12)
##' n <- 200
##' dtS <- sampleData(n,outcome="survival")
##' dtS$X12 <- LETTERS[as.numeric(as.factor(paste0(dtS$X1,dtS$X2)))]
##' dtS <- dtS[dtS$X12!="D"]
##'
##' ## model fit
##' fit <- cph(formula = Surv(time,event)~ X1+X6,data=dtS,y=TRUE,x=TRUE)
##' seqTime <- 1:10
##' ateFit <- ate(fit, data = dtS, treatment = "X1", contrasts = NULL,
##'               times = seqTime, B = 0, iid = TRUE, se = TRUE, verbose = TRUE, band = TRUE)
##'
##' ## display
##' autoplot(ateFit)
##' 
##' ## inference (two sided)
##' statistic <- ateFit$diffRisk$estimate/ateFit$diffRisk$se
##' confint(ateFit, p.value = TRUE, method.band = "bonferroni")$diffRisk
##' confint(ateFit, p.value = TRUE, method.band = "maxT-simulation")$diffRisk
##'
##' anova(ateFit, test = "KS")
##' anova(ateFit, test = "CvM")
##' anova(ateFit, test = "sum")
##'
##' ## manual calculation (one sided)
##' n.sim <- 1e4
##' statistic <- ateFit$diffRisk[, estimate/se]
##' iid.norm <- scale(ateFit$iid$GFORMULA[["1"]]-ateFit$iid$GFORMULA[["0"]],
##'                   scale = ateFit$diffRisk$se)
##'
##' ls.out <- lapply(1:n.sim, function(iSim){
##' iG <- rnorm(NROW(iid.norm))
##' iCurve <- t(iid.norm) %*% iG
##' data.table(max = max(iCurve), L2 = sum(iCurve^2), sum = sum(iCurve),
##' maxC = max(iCurve) - max(statistic),
##' L2C = sum(iCurve^2) - sum(statistic^2),
##' sumC = sum(iCurve) - sum(statistic),
##' sim = iSim)
##' })
##'
##' dt.out <- do.call(rbind,ls.out)
##' dt.out[,.(max = mean(.SD$maxC>=0),
##'           L2 = mean(.SD$L2C>=0),
##'           sum = mean(.SD$sumC>=0))]
##' 
##' ## permutation
##' n.sim <- 250
##' stats.perm <- vector(mode = "list", length = n.sim)
##' pb <- txtProgressBar(max = n.sim, style=3)
##' treatVar <- ateFit$variables["treatment"]
##' 
##' for(iSim in 1:n.sim){ ## iSim <- 1
##' iData <- data.table::copy(dtS)
##' iIndex <- sample.int(NROW(iData), replace = FALSE)
##' iData[, c(treatVar) := .SD[[treatVar]][iIndex]]
##'
##' iFit <- update(fit, data = iData)
##' iAteSim <- ate(iFit, data = iData, treatment = treatVar,
##'                times = seqTime, verbose = FALSE)
##' iStatistic <- iAteSim$diffRisk[,estimate/se]
##' stats.perm[[iSim]] <- cbind(iAteSim$diffRisk[,.(max = max(iStatistic),
##'                                                 L2 = sum(iStatistic^2),
##'                                                 sum = sum(iStatistic))],
##'                             sim = iSim)
##' stats.perm[[iSim]]$maxC <- stats.perm[[iSim]]$max - max(statistic)
##' stats.perm[[iSim]]$L2C <- stats.perm[[iSim]]$L2 - sum(statistic^2)
##' stats.perm[[iSim]]$sumC <- stats.perm[[iSim]]$sum - sum(statistic)
##' setTxtProgressBar(pb, iSim)
##' }
##'
##' dtstats.perm <- do.call(rbind,stats.perm)
##' dtstats.perm[,.(max = mean(.SD$maxC>=0),
##'                 L2 = mean(.SD$L2C>=0),
##'                 sum = mean(.SD$sumC>=0))]
##' }

## * anova.ate (code)
#' @rdname anova.ate
#' @method anova ate
#' @export
anova.ate <- function(object,
                      allContrast = NULL, type = "diff", estimator = object$estimator[1],
                      test = "CvM", transform = NULL,  alternative = "two.sided", n.sim = 1e4,
                      print = TRUE, ...){

    ## ** initialize and check arguments
    if(is.null(allContrast)){
        contrasts <- object$contrasts
        allContrasts <- utils::combn(contrasts, m = 2)
    }
    n.allContrasts <- NCOL(allContrasts)
    type <- match.arg(type, choices = c("diff", "ratio"))
    null <- switch(type, "diff" = 0, "ratio" = 1)
    estimator <- match.arg(estimator, choices = object$estimator)
    if(object$inference$se==FALSE){
        stop("Cannot perform the test without the standard errors \n",
             "Set argument \'se\' to TRUE when calling the ate function \n")
    }
    if(is.null(object$iid)){
        stop("Cannot perform the test without the standard errors \n",
             "Set argument \'iid\' to TRUE when calling the ate function \n")
    }
    if(is.null(transform)){
        if(!is.null(object$transform[[paste0(type,"Risk")]])){
            transform <- "none"
        }else{
            transform <- object$transform[[paste0(type,"Risk")]]
        }
    }
    n.sample <- NROW(object$iid[[1]][[1]])
    n.time <- NCOL(object$iid[[1]][[1]])
    test <- match.arg(test, c("KS","CvM","sum"))
    alternative <- match.arg(alternative, c("two.sided","greater","less"))
    
    ## ** prepare arguments for cpp routine
    beta <- matrix(NA, nrow = n.time, ncol = n.allContrasts)
    iid2cpp <- array(NA, dim = c(n.time, n.sample, n.allContrasts))
    statistic2cpp <- matrix(NA, nrow = n.time, ncol = n.allContrasts)

    typeRisk <- object[[paste0(type,"Risk")]]
    
    for(iC in 1:n.allContrasts){## iC <- 1
        ## *** gather information
        iC.A <- allContrasts[1,iC]
        iC.B <- allContrasts[2,iC]

        iRowIndex <- which((typeRisk$estimator==estimator)*(typeRisk$A==iC.A)*(typeRisk$B==iC.B)==1)
        beta[,iC] <- typeRisk[iRowIndex, .SD$estimate]
        beta.se <- typeRisk[iRowIndex, .SD$se]

        iid.A <- object$iid[[estimator]][[iC.A]]
        iid.B <- object$iid[[estimator]][[iC.B]]
        
        if(type=="diff"){
            iid.AB <- iid.B-iid.A
        }else if(type=="ratio"){
            iRowIndex <- which((object$meanRisk$estimator==estimator)*(object$meanRisk$treatment==iC.A)==1)
            risk.A <- object$meanRisk[iRowIndex, .SD$estimate]
            iRowIndex <- which((object$meanRisk$estimator==estimator)*(object$meanRisk$treatment==iC.B)==1)
            risk.B <- object$meanRisk[iRowIndex, .SD$estimate]
            
            iid.AB <- rowScale_cpp(iid.B, scale = risk.A)-rowScale_cpp(iid.A, scale = risk.B/risk.A^2)
        }
        ## *** transformation
        beta.se <- transformSE(estimate = beta[,iC], se = beta.se, type = transform)
        statistic2cpp[,iC] <- transformT(estimate = beta[,iC], se = beta.se, null = null, type = transform)
        iid2cpp[,,iC] <- t(rowScale_cpp(transformIID(estimate = beta[,iC], iid = iid.AB, type = transform), scale = beta.se))
    }

    ## ** run simulation
    resCpp <- sampleMaxProcess_cpp(nSample = n.sample,
                                   nContrast = n.allContrasts,
                                   nSim = n.sim,
                                   value = statistic2cpp,
                                   iid =  iid2cpp,
                                   global = FALSE,
                                   alternative = switch(alternative,
                                                        "two.sided" = 3,
                                                        "greater" = 2,
                                                        "less" = 1),
                                   type = switch(test, "KS"=1, "CvM"=2, "sum"=3)
                                   )

    ## ** process results
    out <- data.frame(matrix(NA, nrow = n.allContrasts, ncol = 4,
                             dimnames = list(NULL,c("treatment.A","treatment.B","statistic","p.value"))))
    out$treatment.A <- allContrasts[1,]
    out$treatment.B <- allContrasts[2,]
    if(test == "KS"){
        if(alternative == "less"){
            out$statistic <- apply(statistic2cpp,2,min)
        }else if(alternative == "greater"){
            out$statistic <- apply(statistic2cpp,2,max)
        }else if(alternative == "two.sided"){
            out$statistic <- apply(abs(statistic2cpp),2,max)
        }
    }else if(test == "CvM"){
        out$statistic <- colSums(statistic2cpp^2)
    }else if(test == "sum"){
        out$statistic <- colSums(statistic2cpp)
    }
    out$p.value <- colMeans(resCpp>0)

    ## ** display
    if(print){
        txt.alternative <- switch(alternative,
                                  "two.sided" = "unequal mean risk between treatment groups",
                                  "greater" = "greater mean risk with treatment B compared to treatment A",
                                  "less" = "lower mean risk with treatment B compared to treatment A")
        txt.test <- switch(test,
                           "KS" = "maximum",
                           "CvM" = "sum of squares",
                           "sum" = "sum")
        txt.type <- switch(type,
                           "diff" = "differences",
                           "ratio" = "ratios")

        cat("    Comparison of the Average Treatment Effect over all times \n\n")
        cat(" - test statistic        : ",txt.test," of the risk ",txt.type," \n",sep="")
        cat(" - alternative hypothesis: ",txt.alternative," \n",sep="")
        print(out)
    }
    
    ## ** export
    return(invisible(out))
        
}

######################################################################
### anova.ate.R ends here