R/predict.SIS.R

In statcodes/SIS: Sure Independence Screening

#' Model prediction based on a fitted SIS object.
#'
#' Similar to the usual predict methods, this function returns predictions from
#' a fitted \code{'SIS'} object.
#'
#' @export
#' @param object Fitted \code{'SIS'} model object.
#' @param newx Matrix of new values for \code{x} at which predictions are to be
#' made, without the intercept term.
#' @param lambda Penalty parameter \code{lambda} of the final fitted
#' model by (I)SIS at which predictions are required. By default, only the
#' lambda minimizing the criterion \code{tune} is returned.
#' @param which Indices of the penalty parameter \code{lambda} of the final fitted
#' model by (I)SIS at which predictions are required. If supplied, will overwrite
#' the default \code{lambda} value.
#' @param type Type of prediction required. Type \code{'response'} gives the
#' fitted values for \code{'gaussian'}, fitted probabilities for
#' \code{'binomial'}, fitted mean for \code{'poisson'}, and the fitted relative
#' risk for \code{'cox'}. Type \code{'link'} returns the linear predictors for
#' \code{'binomial'}, \code{'poisson'} and \code{'cox'} models; for
#' \code{'gaussian'} models it is equivalent to type \code{'response'}. Type
#' \code{'class'} applies only to \code{'binomial'} models, and produces the
#' class label corresponding to the maximum probability (0-1 labels).
#' @param \dots Not used. Other arguments to predict.
#' @return The object returned depends on type.
#' @author Jianqing Fan, Yang Feng, Diego Franco Saldana, Richard Samworth, and
#' Yichao Wu
#' @seealso \code{\link{SIS}}
#' @references
#' Diego Franco Saldana and Yang Feng (2018) SIS: An R package for Sure Independence Screening in
#' Ultrahigh Dimensional Statistical Models, \emph{Journal of Statistical Software}, \bold{83}, 2, 1-25.
#'
#' Jianqing Fan and Jinchi Lv (2008) Sure Independence Screening
#' for Ultrahigh Dimensional Feature Space (with discussion). \emph{Journal of
#' Royal Statistical Society B}, \bold{70}, 849-911.
#'
#' Jianqing Fan and Rui Song (2010) Sure Independence Screening in Generalized
#' Linear Models with NP-Dimensionality.  \emph{The Annals of Statistics},
#' \bold{38}, 3567-3604.
#'
#' Jianqing Fan, Richard Samworth, and Yichao Wu (2009) Ultrahigh Dimensional
#' Feature Selection: Beyond the Linear Model. \emph{Journal of Machine
#' Learning Research}, \bold{10}, 2013-2038.
#'
#' Jianqing Fan, Yang Feng, and Yichao Wu (2010) High-dimensional Variable
#' Selection for Cox Proportional Hazards Model. \emph{IMS Collections},
#' \bold{6}, 70-86.
#'
#' Jianqing Fan, Yang Feng, and Rui Song (2011) Nonparametric Independence
#' Screening in Sparse Ultrahigh Dimensional Additive Models. \emph{Journal of
#' the American Statistical Association}, \bold{106}, 544-557.
#'
#' Diego Franco Saldana and Yang Feng (2014) SIS: An R package for Sure Independence Screening in
#' Ultrahigh Dimensional Statistical Models, \emph{Journal of Statistical Software}.
#' @keywords models
#' @examples
#'
#'
#' set.seed(0)
#' n <- 400
#' p <- 50
#' rho <- 0.5
#' corrmat <- diag(rep(1 - rho, p)) + matrix(rho, p, p)
#' corrmat[, 4] <- sqrt(rho)
#' corrmat[4, ] <- sqrt(rho)
#' corrmat[4, 4] <- 1
#' corrmat[, 5] <- 0
#' corrmat[5, ] <- 0
#' corrmat[5, 5] <- 1
#' cholmat <- chol(corrmat)
#' x <- matrix(rnorm(n * p, mean = 0, sd = 1), n, p)
#' x <- x %*% cholmat
#' colnames(x) <- unlist(lapply(seq(1:dim(x)[2]), function(y) paste0('V',y)))
#' testX <- matrix(rnorm(10 * p, mean = 0, sd = 1), nrow = 10, ncol = p)
#'
#' # gaussian response
#' set.seed(1)
#' b <- c(4, 4, 4, -6 * sqrt(2), 4 / 3)
#' y <- x[, 1:5] %*% b + rnorm(n)
#' model1 <- SIS(x, y, family = "gaussian", tune = "bic", varISIS = "aggr", seed = 11)
#'
#' predict(model1, testX, type = "response")
#' predict(model1, testX, which = 1:10, type = "response")
#' \dontrun{
#' # binary response
#' set.seed(2)
#' feta <- x[, 1:5] %*% b
#' fprob <- exp(feta) / (1 + exp(feta))
#' y <- rbinom(n, 1, fprob)
#' model2 <- SIS(x, y, family = "binomial", tune = "bic", varISIS = "aggr", seed = 21)
#'
#' predict(model2, testX, type = "response")
#' predict(model2, testX, type = "link")
#' predict(model2, testX, type = "class")
#'
#' predict(model2, testX, which = 1:10, type = "response")
#' predict(model2, testX, which = 1:10, type = "link")
#' predict(model2, testX, which = 1:10, type = "class")
#'
#' # poisson response
#' set.seed(3)
#' b <- c(0.6, 0.6, 0.6, -0.9 * sqrt(2))
#' myrates <- exp(x[, 1:4] %*% b)
#' y <- rpois(n, myrates)
#' model3 <- SIS(x, y, 
#' family = "poisson", penalty = "lasso", tune = "bic", 
#' varISIS = "aggr", seed = 31)
#'
#' predict(model3, testX, type = "response")
#' predict(model3, testX, type = "link")
#' }
#'
predict.SIS <- function(object, newx, lambda = object$lambda, which = NULL, type = c("response", "link", "class"), ...) {
  if (!is.null(which)) {
    lambda <- object$fit$lambda[which]
  }
  if (class(object$fit)[1] == "ncvreg") {
    pred <- predict(object$fit, newx[, object$ix0], lambda = lambda, type = type)
  } else {
    pred <- predict(object$fit, newx[, object$ix0], s = lambda, type = type)
  }
  return(pred)
}
# {
#     if (class(object$fit)[1] == "cv.ncvreg") {
#         family = object$fit$fit$family
#         coefs = object$fit$fit$beta[-1, s]
#         intercept = object$fit$fit$beta[1, s]
#     } else if (class(object$fit)[1] == "cv.glmnet") {
#         family = switch(class(object$fit$glmnet.fit)[1], elnet = "gaussian", lognet = "binomial", fishnet = "poisson",
#             coxnet = "cox")
#         coefs = object$fit$glmnet.fit$beta[, s]
#         intercept = object$fit$glmnet.fit$a0[s]
#     } else if (class(object$fit)[1] == "ncvreg") {
#         family = object$fit$family
#         coefs = object$fit$beta[-1, s]
#         intercept = object$fit$beta[1, s]
#     } else {
#         family = switch(class(object$fit)[1], elnet = "gaussian", lognet = "binomial", fishnet = "poisson",
#             coxnet = "cox")
#         coefs = object$fit$beta[, s]
#         intercept = object$fit$a0[s]
#     }
#
#     newx = newx[, object$ix]
#
#     if (length(coefs) == 1 || length(newx) == 1)
#         eta = newx * coefs else if (length(coefs) == length(s))
#         eta = newx %*% t(coefs) else eta = newx %*% coefs
#
#     if (family != "cox")
#         eta = t(intercept + t(eta))
#
#     if (family == "gaussian")
#         pihat = eta
#     if (family == "binomial")
#         pihat = exp(eta)/(1 + exp(eta))
#     if (family %in% c("poisson", "cox"))
#         pihat = exp(eta)
#
#     if (match.arg(type) == "response")
#         return(pihat)
#     if (match.arg(type) == "link")
#         return(eta)
#     if (match.arg(type) == "class" && family == "binomial")
#         return(matrix(as.integer(eta > 0), ncol = length(s)))
#     if (match.arg(type) == "class" && family == "poisson")
#         stop("Choose response for Poisson, or choose binomial for response")
#
# }