R/HMT.R

In MLGL: Multi-Layer Group-Lasso

#' @title Hierarchical Multiple Testing procedure
#'
#' @description Apply Hierarchical Multiple Testing procedure on a \code{\link{MLGL}} object
#'
#' @param res \code{\link{MLGL}} object
#' @param X matrix of size n*p
#' @param y vector of size n.
#' @param control either "FDR" or "FWER"
#' @param alpha control level for testing procedure
#' @param test test used in the testing procedure. Default is \link{partialFtest}
#' @param addRoot If TRUE, add a common root containing all the groups
#' @param Shaffer If TRUE, a Shaffer correction is performed (only if control = "FWER")
#' @param ... extra parameters for \link{selFDR}
#'
#' @return a list containing:
#' \describe{
#'   \item{lambdaOpt}{lambda values maximizing the number of rejects}
#'   \item{var}{A vector containing the index of selected variables for the first \code{lambdaOpt} value}
#'   \item{group}{A vector containing the values index of selected groups for the first \code{lambdaOpt} value}
#'   \item{selectedGroups}{Selected groups for the first \code{lambdaOpt} value}
#'   \item{indLambdaOpt}{indices associated with optimal lambdas}
#'   \item{reject}{Selected groups for all lambda values}
#'   \item{alpha}{Control level}
#'   \item{test}{Test used in the testing procedure}
#'   \item{control}{"FDR" or "FWER"}
#'   \item{time}{Elapsed time}
#'   \item{hierTest}{list containing the output of the testing function for each lambda. Each element can be used 
#'   with the \link{selFWER} or \link{selFDR} functions.}
#'   \item{lambda}{lambda path}
#'   \item{nGroup}{Number of groups before testing}
#'   \item{nSelectedGroup}{Numer of groups after testing}
#' }
#'
#' @examples
#' set.seed(42)
#' X <- simuBlockGaussian(50, 12, 5, 0.7)
#' y <- X[, c(2, 7, 12)] %*% c(2, 2, -2) + rnorm(50, 0, 0.5)
#' res <- MLGL(X, y)
#'
#' # perform hierarchical testing with FWER control
#' out <- HMT(res, X, y, alpha = 0.05)
#'
#' # test a new value of alpha for a specific lambda
#' selFWER(out$hierTest[[60]], alpha = 0.1)
#' @seealso \link{hierarchicalFWER} \link{hierarchicalFDR} \link{selFWER} \link{selFDR}
#'
#' @export
HMT <- function(res, X, y, control = c("FWER", "FDR"), alpha = 0.05, test = partialFtest, 
                addRoot = FALSE, Shaffer = FALSE, ...) {
  control <- match.arg(control)

  t1 <- proc.time()
  # choose the right function
  hierTestFunction <- hierarchicalFWER
  selFunction <- selFWER
  extraParams <- list(addRoot = addRoot, Shaffer = Shaffer)
  if (control == "FDR") {
    hierTestFunction <- hierarchicalFDR
    selFunction <- selFDR
    extraParams <- list(addRoot = addRoot)
  }

  # testing procedure for each lambda
  REJECT <- TEST <- list()
  nbReject <- rep(0, length(res$lambda))
  prevSelGroup <- selGroup <- c()
  for (i in seq_along(res$lambda))
  {
    # if no groups are selected we do nothing
    if (length(res$group[[i]]) > 0) {
      selGroup <- unique(res$group[[i]])

      # if the selected groups have not changed compared with the last iteration, we copy the result
      if (setequal(prevSelGroup, selGroup)) {
        TEST[[i]] <- TEST[[i - 1]]
        REJECT[[i]] <- REJECT[[i - 1]]
        nbReject[i] <- nbReject[i - 1]
      }
      else {
        # hierarchical testing and selection
        resTest <- do.call(hierTestFunction, c(list(X, y, res$group[[i]], res$var[[i]], test), extraParams))
        
        resSel <- selFunction(resTest, alpha, ...)

        # keep outerNode (need for FDR outer = FALSE, do not change in other cases)
        groupSel <- outerNode(resSel$toSel, resTest$hierMatrix)
        # Id of rejected groups
        REJECT[[i]] <- (resSel$groupId[resSel$toSel])[groupSel]
        TEST[[i]] <- resTest

        # number of rejects for the lambda value
        nbReject[i] <- length(REJECT[[i]])
      }
    } # end if no selection

    prevSelGroup <- selGroup
  } # end for lambda
  t2 <- proc.time()

  hierTestTime <- (t2 - t1)[3]
  time <- c(res$time, hierTestTime)
  names(time) <- c(names(res$time), "test")

  # indices of optimal lambda: the one with the greatest number of reject
  indLambdaOpt <- which(nbReject == max(nbReject))



  # groups selected for the lambda optimal
  indGroupSel <- res$group[[indLambdaOpt[1]]] %in% REJECT[[indLambdaOpt[1]]]
  group <- res$group[[indLambdaOpt[1]]][indGroupSel]
  var <- res$var[[indLambdaOpt[1]]][indGroupSel]


  out <- list(
    lambdaOpt = res$lambda[indLambdaOpt], selectedGroups = REJECT[[indLambdaOpt[1]]], indLambdaOpt = indLambdaOpt, 
    lambda = res$lambda, nGroup = res$nGroup, nSelectedGroup = nbReject, group = group, var = var, test = test, 
    alpha = alpha, reject = REJECT, control = control, time = time, hierTest = TEST
  )

  class(out) <- "HMT"

  return(out)
}