R/predict.R

In HDPenReg: High-Dimensional Penalized Regression

#' Predict response of a new sample Xnew at a given level of penalty
#'
#' @title Prediction of response
#' @author Quentin Grimonprez
#' @param object a LarsParth object
#' @param Xnew a matrix (of size n*object@@p) of covariates.
#' @param lambda If mode ="norm", lambda represents the l1-norm of the coefficients with which we want to predict.
#' If mode="fraction", lambda represents the ratio
#' (l1-norm of the coefficients with which we want to predict)/(l1-norm maximal of the LarsPath object).
#' @param mode "fraction", "lambda" or "norm".
#' @param ... other arguments. Not used.
#' @return The predicted response
#' @aliases predict.LarsPath
#' @method predict LarsPath
#' @examples
#' dataset <- simul(50, 10000, 0.4, 10, 50, matrix(c(0.1, 0.8, 0.02, 0.02), nrow = 2))
#' result <- HDlars(dataset$data[1:40, ], dataset$response[1:40])
#' y <- predict(result, dataset$data[41:50, ], 0.3, "fraction")
#' @export
predict.LarsPath <- function(object, Xnew, lambda, mode = c("fraction", "lambda", "norm"), ...) {
  mode <- match.arg(mode)

  if (missing(object)) {
    stop("object is missing.")
  }
  if (!inherits(object, "LarsPath")) {
    stop("object must be a LarsPath object.")
  }

  if (!is.numeric(lambda)) {
    stop("lambda must be a positive real.")
  }
  if (length(lambda) > 1) {
    stop("lambda must be a positive real.")
  }
  if (lambda < 0) {
    stop("lambda must be a positive real.")
  }

  fraction <- lambda
  if (mode == "norm") {
    fraction <- lambda / object@l1norm[object@nbStep + 1]
  }

  yPred <- rep(object@mu, nrow(Xnew))

  if (mode == "lambda") {
    if (lambda == 0) {
      yPred <- yPred + Xnew[, object@variable[[object@nbStep + 1]]] %*% object@coefficient[[object@nbStep + 1]] - sum(object@meanX[object@variable[[object@nbStep + 1]]] * object@coefficient[[object@nbStep + 1]])
      return(yPred)
    }

    if (lambda >= object@lambda[1]) {
      return(yPred)
    }

    ## fraction >0 and <1
    coeff <- computeCoefficients(object, fraction, mode = "lambda")

    yPred <- yPred + Xnew[, coeff$variable, drop = FALSE] %*% coeff$coefficient - sum(object@meanX[coeff$variable] * coeff$coefficient)

    return(yPred)
  }


  ## fraction = 0 : all coefficients are equal to 0
  if (fraction == 0) {
    return(yPred)
  }

  ## fraction = 1 : coefficients of the last step
  if (fraction >= 1) {
    yPred <- yPred + Xnew[, object@variable[[object@nbStep + 1]]] %*% object@coefficient[[object@nbStep + 1]] - sum(object@meanX[object@variable[[object@nbStep + 1]]] * object@coefficient[[object@nbStep + 1]])
    return(yPred)
  }

  ## fraction >0 and <1
  coeff <- computeCoefficients(object, fraction)

  yPred <- yPred + Xnew[, coeff$variable, drop = FALSE] %*% coeff$coefficient - sum(object@meanX[coeff$variable] * coeff$coefficient)

  return(yPred)
}


#' Compute coefficients at a given level of penalty
#'
#' @title Compute coefficients
#' @author Quentin Grimonprez
#' @param x a LarsParth object
#' @param lambda If mode ="norm", lambda represents the l1-norm of the coefficients with which we want to predict.
#' If mode="fraction", lambda represents the ratio l1-norm/(l1-norm maximal) of the LarsPath object of the coefficients
#' with which we want to predict).
#' @param mode "fraction" or "norm" or "lambda".
#' @return A list containing
#' \describe{
#'   \item{variable}{Index of non-zeros coefficients.}
#'   \item{coefficient}{non-zeros coefficients.}
#' }
#' @examples
#' dataset <- simul(50, 10000, 0.4, 10, 50, matrix(c(0.1, 0.8, 0.02, 0.02), nrow = 2))
#' result <- HDlars(dataset$data[1:40, ], dataset$response[1:40])
#' coeff <- computeCoefficients(result, 0.3, "fraction")
#' @export
computeCoefficients <- function(x, lambda, mode = "fraction") {
  if (missing(x)) {
    stop("x is missing.")
  }
  if (!inherits(x, "LarsPath")) {
    stop("x must be a LarsPath object.")
  }
  if (!(mode %in% c("fraction", "norm", "lambda"))) {
    stop("mode must be \"fraction\" or \"norm\" or \"lambda\".")
  }
  # lambda
  if (!is.numeric(lambda)) {
    stop("lambda must be a positive real.")
  }
  if (length(lambda) > 1) {
    stop("lambda must be a positive real.")
  }
  if (lambda < 0) {
    stop("lambda must be a positive real.")
  }

  if (mode == "fraction") {
    lambda <- lambda * x@l1norm[x@nbStep + 1]
  }

  abscissa <- c()

  if (mode == "lambda") {
    abscissa <- c(x@lambda, 0)
    if (lambda == 0) {
      return(list(variable = x@variable[[x@nbStep + 1]], coefficient = x@coefficient[[x@nbStep + 1]]))
    }
    if (lambda >= x@lambda[1]) {
      return(list(variable = c(), coefficient = c()))
    }
  } else {
    abscissa <- x@l1norm
    if (lambda >= x@l1norm[x@nbStep + 1]) {
      return(list(variable = x@variable[[x@nbStep + 1]], coefficient = x@coefficient[[x@nbStep + 1]]))
    }
    if (lambda == 0) {
      return(list(variable = c(), coefficient = c()))
    }
  }

  index <- 1
  if (mode == "lambda") {
    while (abscissa[index] > lambda) {
      index <- index + 1
    }
  } else {
    while (abscissa[index] < lambda) {
      index <- index + 1
    }
  }

  index <- index - 1

  addId <- c()
  for (i in seq_along(x@addIndex[[index]])) {
    addId <- c(addId, which(x@variable[[index + 1]] == x@addIndex[[index]][i]))
  }

  dropId <- c()
  for (i in seq_along(x@dropIndex[[index]])) {
    dropId <- c(dropId, which(x@variable[[index]] == x@dropIndex[[index]][i]))
  }

  normalId <- c()
  for (i in seq_along(x@variable[[index]])) {
    if (!(x@variable[[index]][i] %in% x@dropIndex[[index]])) {
      normalId <- c(normalId, i)
    }
  }


  coeff <- c()
  if (length(addId) != 0) {
    coeff <- c(
      .computeOrdinate(abscissa[index], abscissa[index + 1], lambda,
                       x@coefficient[[index]][normalId], x@coefficient[[index + 1]][-addId]),
      .computeOrdinate(abscissa[index], abscissa[index + 1], lambda,
                       x@coefficient[[index]][dropId], rep(0, length(dropId))),
      .computeOrdinate(abscissa[index], abscissa[index + 1], lambda,
                       rep(0, length(addId)), x@coefficient[[index + 1]][addId])
    )
  } else {
    coeff <- c(
      .computeOrdinate(abscissa[index], abscissa[index + 1], lambda,
                       x@coefficient[[index]][normalId], x@coefficient[[index + 1]]),
      .computeOrdinate(abscissa[index], abscissa[index + 1], lambda, x@coefficient[[index]][dropId], rep(0, length(dropId)))
    )
  }

  variable <- c(x@variable[[index]][normalId], x@variable[[index]][dropId], x@variable[[index + 1]][addId])

  return(list(variable = variable, coefficient = coeff))
}


.computeOrdinate <- function(x1, x2, x3, y1, y2) {
  return(y1 + (y2 - y1) * ((x3 - x1) / (x2 - x1)))
}

# beta=rep(0,x@p)
#  if(x@fusion)
#  {
#    a=0
#    index=sort(variable,index.return=TRUE)$ix
#    for(i in 1:(length(variable)-1))
#    {
#      a=a+coefficient[index[i]]
#      beta[variable[index[i]]:(variable[index[i+1]]-1)]=rep(a, variable[index[i+1]]-variable[index[i]])
#
#    }
#    a=a+coefficient[index[length(index)]]
#    beta[variable[index[length(index)]]:x@p]=rep(a,x@p-variable[index[length(index)]]+1)
#  }