R/grpenet.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Cross Validate Group Elastic Net Regression
#'
#' @param x the model matrix
#' @param y the outcome
#' @param cv.method preferably one of "boot632" (the default), "cv", or "repeatedcv".
#' @param nfolds the number of bootstrap or cross-validation folds to use. defaults to 5.
#' @param nrep the number of repetitions for cv.method = "repeatedcv". defaults to 4.
#' @param tunlen the number of values for the unknown hyperparameter to test. defaults to 10.
#' @param crit the criterion by which to evaluate the model performance. must be one of "MAE" (the default)
#' or "MSE".
#' @param max.c the largest value of the constant for calculating lambda. defaults to 8, but
#' may be adjusted. for example, if the error metric becomes constant after a certain
#' value of C, it may be advisable to lower max.c to a smaller value to obtain
#' a more fine-grained grid over the plausible values.
#'
#' @return
#' a train object
#' @export
#'
cv_grpenet = function(x, y, idx, cv.method = "boot632", nfolds = 5, nrep = 4, tunlen = 10, crit = "MAE", max.c = 8){

  GRPENET <- list(type = "Regression",
                library = "grpreg",
                loop = NULL)

  GRPENET$parameters <- data.frame(parameter = c("C", "alpha", "base.lambda"),
                                 class = rep("numeric", 3),
                                 label = c("C", "alpha", "base.lambda"))

  huber.scale = function(y){
    MASS::hubers(y,
                 initmu =
                   MASS::hubers(y,
                          initmu = MASS::hubers(y)$mu,
                          s = sqrt(mean(c(sd(y)^2, mad(y)^2)))
                   )$mu)$s
  }

  lm.betas <- lmSolve(y ~ . , data = cbind.data.frame(y = y, x))
  model.mat <- model.matrix(y ~ . , data = cbind.data.frame(y = y, x))
  lm.pred <- as.vector(lm.betas) %*% t(model.mat)
  lm.res <- as.vector(model.frame(y ~ . , data = cbind.data.frame(y = y, x))[,1]) - lm.pred
  GRPENET$noiseSD <- huber.scale(lm.res)
  GRPENET$max.c <- max.c
  GRPENET$idx <- idx
  GRPENET$columns <- colnames(x)
  grpenetGrid <- function(x, y, idx = GRPENET$idx, max.c = GRPENET$max.c, noise.sd = GRPENET$noiseSD, len = NULL, search = "grid") {

    D = nrow(x)
    N = length(y)
    lambda0 = noise.sd * sqrt(log(D) / N)

    C <- seq(1, max.c, length.out = len)
    alphas <- seq(0.0001, 1, len = 6)

    ## use grid search:
    if(search == "grid"){
      search = "grid"
    } else {
      search = "grid"
    }

    grid <- expand.grid(C = C, alpha = alphas)
    grid$base.lambda <- rep(lambda0, nrow(grid))

    out <- grid
    return(out)
  }

  GRPENET$grid <- grpenetGrid

  grpenetFit <- function(x, y,  idx = GRPENET$idx, colms = GRPENET$columns, param, ...) {

    f = function(alpha){
      sqrt((2 - alpha)^3 / (2 - (1 - alpha)))
    }

    grpreg::grpreg(
      X = as.matrix(x),
      y = as.vector(y),
      penalty = "grLasso",
      group = idx,
      family = "gaussian",
      lambda = param$C * f(param$alpha) * param$base.lambda,
      alpha = param$alpha)
  }

  GRPENET$fit <- grpenetFit
  GRPENET$prob <- grpenetFit

  grpenetPred <- function(modelFit, newdata, preProc = NULL, submodels = NULL){
    grpreg:::predict.grpreg(modelFit, newdata, exact = TRUE)
  }

  GRPENET$predict <- grpenetPred

  postRobResamp = function(pred, obs) {

    isNA <- is.na(pred)
    pred <- pred[!isNA]
    obs <- obs[!isNA]
    if (!is.factor(obs) && is.numeric(obs)) {
      if (length(obs) + length(pred) == 0) {
        out <- rep(NA, 2)
      }
      else {

        huber.mean = function(y) {
          MASS::hubers(y, initmu =
                         MASS::hubers(y,
                                      initmu = MASS::hubers(y, s = sd(y))$mu,
                                      s = sqrt(mean(c(sd(y)^2, mad(y)^2))))$mu)$mu
        }

        robmse <- huber.mean((pred - obs)^2)
        robmae <- huber.mean(abs(pred - obs))
        out <- c(robmse, robmae)
      }
      names(out) <- c("MSE", "MAE")
    }
    else {
      if (length(obs) + length(pred) == 0) {
        out <- rep(NA, 2)
      }
      else {
        pred <- factor(pred, levels = levels(obs))
        requireNamespaceQuietStop("e1071")
        out <- unlist(e1071::classAgreement(table(obs, pred)))[c("diag",
                                                                 "kappa")]
      }
      names(out) <- c("Accuracy", "Kappa")
    }
    if (any(is.nan(out)))
      out[is.nan(out)] <- NA
    out
  }

  robustSummary = function (data, lev = NULL, model = NULL)
  {
    if (is.character(data$obs))
      data$obs <- factor(data$obs, levels = lev)
    postRobResamp(data[, "pred"], data[, "obs"])
  }


  if (cv.method == "repeatedcv") {
    fitControl <- trainControl(method = cv.method,
                               number = nfolds,
                               repeats = nrep,
                               savePredictions = "all",
                               summaryFunction = robustSummary,
                               search = "grid")
  } else {

    fitControl <- trainControl(method = cv.method,
                               number = nfolds,
                               savePredictions = "all",
                               summaryFunction = robustSummary,
                               search = "grid")
  }

  Formula <- y ~ .
  Dat <- cbind.data.frame(y = y, x)
  fitted.models <- train(form = Formula,
                         data = Dat,
                         method = GRPENET,
                         metric = crit,
                         tuneLength = tunlen,
                         maximize = FALSE,
                         preProcess = c("center", "scale"),
                         trControl = fitControl)

  f = function(alpha){
    sqrt((2 - alpha)^3 / (2 - (1 - alpha)))
  }

  lambda <- fitted.models$results$C * f(fitted.models$results$alpha) * fitted.models$results$base.lambda
  fitted.models$results <- cbind.data.frame(fitted.models$results[,1:3],
                                            lambda = lambda,
                                            fitted.models$results[,4:ncol(fitted.models$results)])

  return(fitted.models)

}