R/S4-BuyseTest-sensitivity.R

In BuyseTest: Generalized Pairwise Comparisons

### sensitivity.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: mar 31 2021 (14:07) 
## Version: 
## Last-Updated: Mar 20 2023 (12:08) 
##           By: Brice Ozenne
##     Update #: 316
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

## * Documentation - sensitivity
#' @docType methods
#' @name S4BuyseTest-sensitivity
#' @title  Sensitivity Analysis for the Choice of the Thresholds
#' @aliases sensitivity,S4BuyseTest-method
#' @include S4-BuyseTest.R
#' 
#' @description Evaluate the statistic of interest along various thresholds of clinical relevance.
#' @param object an \R object of class \code{\linkS4class{S4BuyseTest}}, i.e., output of \code{\link{BuyseTest}}
#' @param threshold [list] a list containing for each endpoint the thresholds to be considered.
#' @param statistic [character] the statistic summarizing the pairwise comparison:
#' \code{"netBenefit"} displays the net benefit, as described in Buyse (2010) and Peron et al. (2016)),
#' \code{"winRatio"} displays the win ratio, as described in Wang et al. (2016),
#' \code{"favorable"} displays the proportion in favor of the treatment (also called Mann-Whitney parameter), as described in Fay et al. (2018).
#' \code{"unfavorable"} displays the proportion in favor of the control.
#' Default value read from \code{BuyseTest.options()}.
#' @param null [numeric] right hand side of the null hypothesis (used for the computation of the p-value).
#' @param conf.level [numeric] confidence level for the confidence intervals.
#' Default value read from \code{BuyseTest.options()}.
#' @param alternative [character] the type of alternative hypothesis: \code{"two.sided"}, \code{"greater"}, or \code{"less"}.
#' Default value read from \code{BuyseTest.options()}.
#' @param band [logical] should simulateneous confidence intervals be computed?
#' @param adj.p.value [logical] should p-value adjusted for multiple comparisons be computed?
#' @param transformation [logical]  should the CI be computed on the logit scale / log scale for the net benefit / win ratio and backtransformed.
#' Otherwise they are computed without any transformation.
#' Default value read from \code{BuyseTest.options()}. Not relevant when using permutations or percentile bootstrap.
#' @param cpus [integer, >0] the number of CPU to use. Default value is 1.
#' @param trace [logical] Should the execution of the function be traced?
#' @param ... argument passsed to the function \code{transformCIBP} of the riskRegression package.
#'
#' @details Simulateneous confidence intervals and adjusted p-values are computed using a single-step max-test approach via the function \code{transformCIBP} of the riskRegression package.

#' @examples
#' 
#' \dontrun{
#' require(ggplot2)
#' 
#' ## simulate data
#' set.seed(10)
#' df.data <- simBuyseTest(1e2, n.strata = 2)
#'
#' ## with one endpoint
#' ff1 <- treatment ~ TTE(eventtime, status = status, threshold = 0.1)
#' BT1 <- BuyseTest(ff1, data= df.data)
#' se.BT1 <- sensitivity(BT1, threshold = seq(0,2,0.25), band = TRUE)
#' autoplot(se.BT1)
#'
#' ## with two endpoints
#' ff2 <- update(ff1, .~. + cont(score, threshold = 1))
#' BT2 <- BuyseTest(ff2, data= df.data)
#' se.BT2 <- sensitivity(BT2, threshold = list(eventtime = seq(0,2,0.25), score = 0:2),
#'                       band = TRUE)
#' autoplot(se.BT2)
#' autoplot(se.BT2, col = NA)
#' }


## * Method - sensitivity
#' @rdname S4BuyseTest-sensitivity
setMethod(f = "sensitivity",
          signature = "S4BuyseTest",
          definition = function(object, threshold,
                                statistic = NULL, band = FALSE, conf.level = NULL, null = NULL, transformation = NULL, alternative = NULL, adj.p.value = FALSE,
                                trace = TRUE, cpus = 1, ...){

              ## ** normalize user input
              ## band
              if(object@method.inference!="u-statistic"){
                  stop("Cannot compute confidence bands when \'method.inference\' used to obtain the object is not \"u-statistic\". \n")
              }
              
              ## endpoint
              name.endpoint <- object@endpoint
              n.endpoint <- length(name.endpoint)
              option <- BuyseTest.options()
              if(is.null(statistic)){
                  statistic <- option$statistic
              }
              if(is.null(transformation)){
                  transformation <- option$transformation
              }
              if(is.null(conf.level)){
                  conf.level <- option$conf.level
              }
              if(is.null(alternative)){
                  alternative <- option$alternative
              }
              if(is.null(null)){
                  null <- switch(statistic,
                                 "netBenefit" = 0,
                                 "winRatio" = 1,
                                 "favorable" = 1/2,
                                 "unfavorable" = 1/2)
              }else{
                  validNumeric(null, valid.length = 1,
                               min = if("statistic"=="netBenefit"){-1}else{0},
                               max = if("statistic"=="winRatio"){Inf}else{1})
              }

              ## threshold
              if(is.matrix(threshold) || is.data.frame(threshold)){
                  if(is.matrix(threshold)){
                      threshold <- as.data.frame(threshold)
                  }
                  if(NCOL(threshold)!=n.endpoint){
                      stop("When a matrix, the argument \'threshold\' should contain ",n.endpoint," columns (and not ",length(threshold),"). \n",
                           "Each column corresponds to a prioritized endpoint. \n")
                  }
                  grid.threshold <- threshold
                  if(is.null(names(grid.threshold))){
                      names(grid.threshold) <- names(name.endpoint)
                  }else{
                      if(any(duplicated(names(grid.threshold)))){
                          stop("Duplicated column names in argument \"threshold\": \"",paste0(names(grid.threshold)[duplicated(names(grid.threshold))], collapse= "\" \""),"\".\n")
                      }
                      if(any(names(grid.threshold) %in% names(name.endpoint) == FALSE)){
                          stop("Incorrect column names in argument \"threshold\": \"",paste0(names(grid.threshold)[names(grid.threshold) %in% names(name.endpoint) == FALSE], collapse= "\" \""),"\".\n",
                               "Valid names: \"",paste0(setdiff(names(name.endpoint), names(grid.threshold)), collapse= "\" \""),"\".\n")
                      }
                      grid.threshold <- grid.threshold[,names(name.endpoint)]
                  }
                  
              }else if(is.list(threshold) || is.vector(threshold)){

                  if(any(duplicated(name.endpoint))){
                      stop("Argument \'threshold\' must be a matrix when some endpoints are repeteadly used (i.e at different priorities) with different thresholds. \n")
                  }

                  if(!is.list(threshold)){
                      threshold <- list(threshold)
                  }

                  if(is.null(names(threshold))){
                      if(length(threshold)!=n.endpoint){
                          stop("When a list, the argument \'threshold\' should have length ",n.endpoint," (and not ",length(threshold),"). \n",
                               "Each element of the list corresponds to a prioritized endpoint. \n")
                      }
                  }else{
                      if(any(duplicated(names(threshold)))){
                          stop("Argument \'threshold\' must not contain duplicated names. \n",
                               "Duplicated names: \"",paste0(names(threshold)[duplicated(names(threshold))], collapse = "\" \""),"\". \n")
                      }
                      if(any(names(threshold) %in% name.endpoint == FALSE)){
                          stop("Some names used in the argument \'threshold\' does not match the existing endpoints. \n",
                               "Incorrect names: \"",paste0(names(threshold)[names(threshold) %in% name.endpoint == FALSE], collapse = "\" \""),"\". \n",
                               "Possible names: \"",paste0(setdiff(name.endpoint,names(threshold)), collapse = "\" \""),"\". \n")
                      }
                      threshold.save <- threshold
                      threshold <- setNames(vector(mode = "list", length = n.endpoint), name.endpoint)
                      threshold[names(threshold.save)] <- threshold.save
                  }

                  if(any(sapply(threshold,length)==0)){
                      threshold[sapply(threshold,length)==0] <- object@threshold[sapply(threshold,length)==0]
                  }
              
                  grid.threshold <- expand.grid(threshold)
                  colnames(grid.threshold) <- name.endpoint
              
              }else{
                  stop("Argument \'threshold\' should be a list or a matrix \n")
              }

              ## formula
              ls.args <- object@call
              if("formula" %in% names(ls.args)){
                  ls.args$formula <- NULL

                  args.tempo <- initializeFormula(object@call$formula, hierarchical = object@hierarchical)
                  ls.args[names(args.tempo)] <- args.tempo

              }
              ls.args$trace <- 0

              if (cpus == "all") { 
                  cpus <- parallel::detectCores() # this function detect the number of CPU cores 
              }

              if(band && any(object@weightObs!=1)){
                  stop("Confidence bands cannot not currently be derived with weighted observations. \n") 
              }

              ## ** run BuyseTest
              n.se <- NROW(grid.threshold)
              test.varying <- apply(grid.threshold,2,function(iX){length(unique(iX))!=1})
              if(all(test.varying==FALSE)){
                  stop("Only a single combination of thresholds. No need for a sensitivity analysis.\n")
              }
              gridRed.threshold <- grid.threshold[,which(test.varying),drop=FALSE]
                  

              if(trace>0){cat("Run ",n.se," GPC analyses: \n", sep = "")}
              
              if (cpus == 1) {
                  ls.confint <- vector(mode="list", length = n.se)
                  ls.iid <- vector(mode="list", length = n.se)

                  if(trace>0){pb <- utils::txtProgressBar(max = n.se, style = 3)}

                  for(iSe in 1:n.se){
                      if(trace>0){utils::setTxtProgressBar(pb, iSe)}
                      iLS.args <- ls.args
                      iLS.args$threshold <- as.double(grid.threshold[iSe,])
                      iBT <- do.call(BuyseTest, args = iLS.args)

                      iConfint <- confint(iBT, statistic = statistic, null = null, conf.level = conf.level, alternative = alternative, transformation = transformation)[n.endpoint,]
                      ls.confint[[iSe]] <- data.frame(c(gridRed.threshold[iSe,,drop=FALSE], iConfint))
                      if(iBT@method.inference=="u-statistic"){
                          ls.iid[[iSe]] <- getIid(iBT, statistic = statistic)
                      }
                  }

                  if(trace>0){close(pb)}

              }else{
                  if(trace>0){
                      cl <- suppressMessages(parallel::makeCluster(cpus, outfile = ""))
                      pb <- utils::txtProgressBar(max = n.se, style = 3)          
                  }else{
                      cl <- parallel::makeCluster(cpus)
                  }
                  test.lazyeval <- sapply(ls.args,function(x){inherits(x,"name")})
                  if(any(test.lazyeval)){
                      toExport <- unlist(lapply(ls.args[test.lazyeval],deparse))
                  }else{
                      toExport <- NULL
                  }

                  i <- NULL ## [:forCRANcheck:] foreach
                  ls.sensitivity <- foreach::`%dopar%`(
                                                 foreach::foreach(i=1:n.se, .export = toExport), {                                           
                                                     if(trace>0){utils::setTxtProgressBar(pb, i)}
                                                     iLS.args <- ls.args
                                                     iLS.args$threshold <- as.double(grid.threshold[i,])
                                                     iBT <- do.call(BuyseTest, args = iLS.args)

                                                     iConfint <- confint(iBT, statistic = statistic, null = null, conf.level = conf.level, alternative = alternative, transformation = transformation)[n.endpoint,]
                                                     iOut <- list(confint = data.frame(c(gridRed.threshold[i,,drop=FALSE], iConfint)))
                                                     if(iBT@method.inference=="u-statistic"){
                                                         iOut[["iid"]] <- getIid(iBT, statistic = statistic)
                                                     }
                                                     return(iOut)
                                                 })

                  parallel::stopCluster(cl)
                  if(trace>0){close(pb)}

                  ls.confint <- lapply(ls.sensitivity,"[[","confint")
                  if(object@method.inference=="u-statistic"){
                      ls.iid <- lapply(ls.sensitivity,"[[","iid")
                  }
              }
              
              df.confint <- as.data.frame(do.call(rbind,ls.confint))
              if(object@method.inference=="u-statistic"){
                  attr(df.confint, "iid") <- do.call(cbind,ls.iid)
              }
              
              ## ** compute confidence bands
              if(band){
                  requireNamespace("riskRegression")
                  A.iid <- array(NA, dim = c(NROW(attr(df.confint, "iid")), NCOL(attr(df.confint, "iid")),1))
                  A.iid[,,1] <- attr(df.confint, "iid")

                  min.value <- switch(statistic,
                                      "netBenefit" = -1,
                                      "winRatio" = 0,
                                      "favorable" = 0,
                                      "unfavorable" = 0)
                  max.value <- switch(statistic,
                                      "netBenefit" = 1,
                                      "winRatio" = Inf,
                                      "favorable" = 1,
                                      "unfavorable" = 1)

                  ## temporary fix: the next few lines should be remove when riskRegression will be updated
                  if(statistic %in% c("none","netBenefit","winRatio") || !inherits(try(riskRegression::transformCIBP(estimate = 1, se = 1, type = "atanh2", seed = NA, band = FALSE, alternative = "two.sided"),silent=TRUE),"try-error")){
                      type <- switch(statistic,
                                     "netBenefit" = "atanh",
                                     "winRatio" = "log",
                                     "favorable" = "cloglog",## note: not the same transformation as confint
                                     "unfavorable" = "cloglog", ## note: not the same transformation as confint
                                     "none" = "none") 
                  }else{
                      type <- switch(statistic,
                                     "netBenefit" = "atanh",
                                     "winRatio" = "log",
                                     "favorable" = "atanh2",
                                     "unfavorable" = "atanh2",
                                     "none" = "none") 
                  }
                  
                  dots <- list(...)
                  if("seed" %in% names(dots) == FALSE){
                      dots$seed <- NA
                  }
                  if("method.band" %in% names(dots) == FALSE){
                      dots$method.band <- "maxT-integration"
                  }

                  iBand <- do.call(riskRegression::transformCIBP,
                                   args = c(list(estimate = rbind(df.confint$estimate[df.confint$se>0]),
                                                 se = rbind(df.confint$se[df.confint$se>0]),
                                                 iid = A.iid[,df.confint$se>0,,drop=FALSE],
                                                 null = null,
                                                 conf.level = conf.level,
                                                 alternative = alternative,
                                                 ci = TRUE, type = type, min.value = min.value, max.value = max.value,
                                                 band = TRUE, p.value = adj.p.value),
                                            dots))
                  attr(df.confint,"quantileBand") <- iBand$quantile
                  df.confint$lower.band <- rep(0,length(df.confint$se))
                  df.confint$lower.band[df.confint$se>0] <- iBand$lowerBand[1,]
                  df.confint$upper.band <- rep(0,length(df.confint$se))
                  df.confint$upper.band[df.confint$se>0] <- iBand$upperBand[1,]
                  if(adj.p.value==TRUE){
                      df.confint$adj.p.value <- rep(1,length(df.confint$se))
                      df.confint$adj.p.value[df.confint$se>0] <- iBand$adj.p.value[1,]
                  }
         
                  ## iBand <- do.call(riskRegression::transformCIBP,
                  ##                  args = c(list(estimate = rbind(df.confint$estimate),
                  ##                                se = rbind(df.confint$se),
                  ##                                iid = A.iid,
                  ##                                null = null,
                  ##                                conf.level = conf.level,
                  ##                                alternative = alternative,
                  ##                                ci = TRUE, type = type, min.value = min.value, max.value = max.value,
                  ##                                band = TRUE, p.value = adj.p.value),
                  ##                           dots))

                  ## attr(df.confint,"quantileBand") <- iBand$quantile
                  ## df.confint$lower.band <- iBand$lowerBand[1,]
                  ## df.confint$upper.band <- iBand$upperBand[1,]
                  ## if(adj.p.value==TRUE){
                  ##     df.confint$adj.p.value <- iBand$adj.p.value[1,]
                  ## }
              }

              ## ** export
              attr(df.confint,"statistic") <- statistic
              attr(df.confint,"grid") <- grid.threshold
              attr(df.confint,"gridRed") <- gridRed.threshold
              class(df.confint) <- append("sensitivity",class(df.confint))
              return(df.confint)

          })


##----------------------------------------------------------------------
### sensitivity.R ends here