R/final.selection.R

In HCTR: Higher Criticism Tuned Regression

#' @title Final Selection
#' @description Returns the index of final selected variables in the final chosen model.
#'
#' @param cv.fit An object of either class "cv.glmnet" from glmnet::cv.glmnet() or class "cv.ncvreg" from ncvreg::cv.ncvreg(), which is a list generated by a cross-validation fit.
#' @param pihat eatimated proprtion from HCTR::est.prop().
#' @param p Total number of variables, except for covariates.
#' @param cov.num Number of covariates in model, default is 0. Covariate matrix, W, is assumed on the left side of variable matrix, X. The column index of covariates are before those of variables.
#'
#' @return A sequence of index of final selected variables in the final chosen model.
#'
#' @export
#' @import glmnet
#' @import ncvreg
#' @importFrom stats coef
#'
#' @examples
#' set.seed(10)
#' X <- matrix(rnorm(20000), nrow = 100)
#' beta <- rep(0, 200)
#' beta[1:100] <- 5
#' Y <- MASS::mvrnorm(n = 1, mu = X%*%beta, Sigma = diag(100))
#' fit <- glmnet::cv.glmnet(x = X, y = Y)
#' pihat <- 0.01
#' result <- est.lambda(cv.fit = fit, pihat = pihat, p = ncol(X))
#' lambda.seq <- seq(from = result$lambda.min, to = result$lambda.max, length.out = 100)
#' # Note: The lambda sequences in glmnet and ncvreg are diffrent.
#' fit2 <- glmnet::cv.glmnet(x = X, y = Y, lambda = lambda.seq)
#' result2 <- final.selection(cv.fit = fit2, pihat = 0.01, p = ncol(X))
final.selection <- function(cv.fit, pihat, p, cov.num = 0) {
  upper <- ceiling(pihat*p)
  lower <- floor(pihat*p)
  if (cov.num > 0) {
    upper <- (upper + cov.num)
    lower <- (lower + cov.num)
  }
  if (class(cv.fit) == "cv.glmnet") {
    lambda.index <- which((cv.fit$nzero == lower)|(cv.fit$nzero == upper))
    lambda_hat <- cv.fit$lambda[lambda.index[which.min(cv.fit$cvm[lambda.index])]]
    beta.est <- coef(cv.fit, s = lambda_hat)
    selected.index <- beta.est@i[-1]
  } else if (class(cv.fit) == "cv.ncvreg") {
    nzero <- (colSums(cv.fit$fit$beta != 0) - 1)
    lambda.index <- which((nzero == lower)|(nzero == upper))
    lambda_hat <- cv.fit$lambda[lambda.index[which.min(cv.fit$cve[lambda.index])]]
    beta.est <- coef(cv.fit, lambda=lambda_hat)
    selected.index <- which(beta.est!=0,arr.ind = T)[-1]-1
  } else {
    stop("Invalid fitted model!")
  }

  return (selected.index)
}