R/Prediction_Functions.R

In SplitGLM: Split Generalized Linear Models

#' 
#' @title Predictions for SplitGLM Object
#' 
#' @description \code{predict.SplitGLM} returns the predictions for a SplitGLM object.
#' 
#' @method predict SplitGLM
#'
#' @param object An object of class SplitGLM.
#' @param newx New data for predictions.
#' @param group_index The group for which to return the coefficients. Default is the ensemble.
#' @param type The type of predictions for binary response. Options are "prob" (default) and "class".
#' @param ... Additional arguments for compatibility.
#' 
#' @return The predictions for the SplitGLM object.
#' 
#' @export
#' 
#' @author Anthony-Alexander Christidis, \email{anthony.christidis@stat.ubc.ca}
#' 
#' @seealso \code{\link{SplitGLM}}
#' 
#' @examples
#' \donttest{
#' # Data simulation
#' set.seed(1)
#' n <- 50
#' N <- 2000
#' p <- 1000
#' beta.active <- c(abs(runif(p, 0, 1/2))*(-1)^rbinom(p, 1, 0.3))
#' # Parameters
#' p.active <- 100
#' beta <- c(beta.active[1:p.active], rep(0, p-p.active))
#' Sigma <- matrix(0, p, p)
#' Sigma[1:p.active, 1:p.active] <- 0.5
#' diag(Sigma) <- 1
#' 
#' # Train data
#' x.train <- mvnfast::rmvn(n, mu = rep(0, p), sigma = Sigma) 
#' prob.train <- exp(x.train %*% beta)/
#'               (1+exp(x.train %*% beta))
#' y.train <- rbinom(n, 1, prob.train)
#' mean(y.train)
#' # Test data
#' x.test <- mvnfast::rmvn(N, mu = rep(0, p), sigma = Sigma)
#' prob.test <- exp(x.test %*% beta)/
#'              (1+exp(x.test %*% beta))
#' y.test <- rbinom(N, 1, prob.test)
#' mean(y.test)
#' 
#' # SplitGLM - CV (Multiple Groups)
#' split.out <- SplitGLM(x.train, y.train,
#'                       glm_type="Logistic",
#'                       G=10, include_intercept=TRUE,
#'                       alpha_s=3/4, alpha_d=1,
#'                       lambda_sparsity=1, lambda_diversity=1,
#'                       tolerance=1e-3, max_iter=1e3,
#'                       active_set=FALSE)
#' split.coef <- coef(split.out)
#' # Predictions
#' split.prob <- predict(split.out, newx=x.test, type="prob", group_index=NULL)
#' split.class <- predict(split.out, newx=x.test, type="class", group_index=NULL)
#' plot(prob.test, split.prob, pch=20)
#' abline(h=0.5,v=0.5)
#' mean((prob.test-split.prob)^2)
#' mean(abs(y.test-split.class))
#' 
#' }
#' 
predict.SplitGLM <- function(object, newx, group_index = NULL, type = c("prob", "class")[1], ...){
  
  # Check input data
  if(!any(class(object) %in% "SplitGLM"))
    stop("The object should be of class \"SplitGLM\"")
  
  split.coef <- coef(object, group_index=group_index)
  
  if(object$glm_type=="Linear"){
    
    return(split.coef[1] + newx %*% split.coef[-1])
    
  } else if(object$glm_type=="Logistic"){
    
    if(!(type %in% c("prob", "class")))
      stop("The variable \"type\" must be one of: \"prob\", or \"class\".")
    
    logistic.prob <- exp(split.coef[1] + newx %*% split.coef[-1])/(1+exp(split.coef[1] + newx %*% split.coef[-1]))
    
    if(type=="prob")
      return(logistic.prob) else if(type=="class")
        return(round(logistic.prob, 0))
    
  } else if(object$glm_type=="Gamma"){
    
    return(-1/(split.coef[1] + newx %*% split.coef[-1]))
    
  } else if(object$glm_type=="Poisson"){
    
    return(exp(split.coef[1] + newx %*% split.coef[-1]))
  
  }
  
}
#'
#' @title Predictions for cv.SplitGLM Object
#' 
#' @description \code{predict.cv.SplitGLM} returns the predictions for a SplitGLM object.
#' 
#' @method predict cv.SplitGLM
#'
#' @param object An object of class cv.SplitGLM.
#' @param newx New data for predictions.
#' @param group_index The group for which to return the coefficients. Default is the ensemble.
#' @param type The type of predictions for binary response. Options are "prob" (default) and "class".
#' @param ... Additional arguments for compatibility.
#' 
#' @return The predictions for the cv.SplitGLM object.
#' 
#' @export
#' 
#' @author Anthony-Alexander Christidis, \email{anthony.christidis@stat.ubc.ca}
#' 
#' @seealso \code{\link{cv.SplitGLM}}
#' 
#' @examples 
#' \donttest{
#' # Data simulation
#' set.seed(1)
#' n <- 50
#' N <- 2000
#' p <- 1000
#' beta.active <- c(abs(runif(p, 0, 1/2))*(-1)^rbinom(p, 1, 0.3))
#' # Parameters
#' p.active <- 100
#' beta <- c(beta.active[1:p.active], rep(0, p-p.active))
#' Sigma <- matrix(0, p, p)
#' Sigma[1:p.active, 1:p.active] <- 0.5
#' diag(Sigma) <- 1
#' 
#' # Train data
#' x.train <- mvnfast::rmvn(n, mu = rep(0, p), sigma = Sigma) 
#' prob.train <- exp(x.train %*% beta)/
#'               (1+exp(x.train %*% beta))
#' y.train <- rbinom(n, 1, prob.train)
#' mean(y.train)
#' # Test data
#' x.test <- mvnfast::rmvn(N, mu = rep(0, p), sigma = Sigma)
#' prob.test <- exp(x.test %*% beta)/
#'              (1+exp(x.test %*% beta))
#' y.test <- rbinom(N, 1, prob.test)
#' mean(y.test)
#' 
#' # SplitGLM - CV (Multiple Groups)
#' split.out <- cv.SplitGLM(x.train, y.train,
#'                          glm_type="Logistic",
#'                          G=10, include_intercept=TRUE,
#'                          alpha_s=3/4, alpha_d=1,
#'                          n_lambda_sparsity=50, n_lambda_diversity=50,
#'                          tolerance=1e-3, max_iter=1e3,
#'                          n_folds=5,
#'                          active_set=FALSE,
#'                          n_threads=1)
#' split.coef <- coef(split.out)
#' # Predictions
#' split.prob <- predict(split.out, newx=x.test, type="prob", group_index=NULL)
#' split.class <- predict(split.out, newx=x.test, type="class", group_index=NULL)
#' plot(prob.test, split.prob, pch=20)
#' abline(h=0.5,v=0.5)
#' mean((prob.test-split.prob)^2)
#' mean(abs(y.test-split.class))
#' 
#' }
#' 
predict.cv.SplitGLM <- function(object, newx, group_index = NULL, type = c("prob", "class")[1], ...){
  
  # Check input data
  if(!any(class(object) %in% "cv.SplitGLM"))
    stop("The object should be of class \"cv.SplitGLM\".")
  
  split.coef <- coef(object, group_index=group_index)
  
  if(object$glm_type=="Linear"){
    
    return(split.coef[1] + newx %*% split.coef[-1])
    
  } else if(object$glm_type=="Logistic"){
    
    if(!(type %in% c("prob", "class")))
      stop("The variable \"type\" must be one of: \"prob\", or \"class\".")
    
    logistic.prob <- exp(split.coef[1] + newx %*% split.coef[-1])/(1+exp(split.coef[1] + newx %*% split.coef[-1]))
    
    if(type=="prob")
      return(logistic.prob) else if(type=="class")
        return(round(logistic.prob, 0))
    
  } else if(object$glm_type=="Gamma"){
    
    return(-1/(split.coef[1] + newx %*% split.coef[-1]))
    
  } else if(object$glm_type=="Poisson"){
    
    return(exp(split.coef[1] + newx %*% split.coef[-1]))
    
  }
  
}