R/extract.R

In robregcc: Robust Regression with Compositional Covariates

#' Extract coefficients estimate from the sparse version of the robregcc fitted object.
#'
#' S3 methods extracting estimated coefficients for objects generated by
#' \code{robregcc}. Robust coeffcient estimate.
#'
#' @param object Object generated by \code{robregcc}.
#' @param type 0/1 residual estimate before/after sanity check
#' @param s 0/1 no/yes 1se rule
#' @return coefficient estimate
#' @importFrom stats coef
#' @export
#' @examples  
#' 
#' library(magrittr)
#' library(robregcc)
#' 
#' data(simulate_robregcc)
#' X <- simulate_robregcc$X;
#' y <- simulate_robregcc$y
#' C <- simulate_robregcc$C
#' n <- nrow(X); p <- ncol(X); k <-  nrow(C)
#' 
#' Xt <- cbind(1,X)                         # accounting for intercept in predictor
#' C <- cbind(0,C)                           # accounting for intercept in constraint
#' bw <- c(0,rep(1,p))                       # weight matrix to not penalize intercept 
#' 
#' example_seed <- 2*p+1               
#' set.seed(example_seed) 
#' \donttest{
#' # Breakdown point for tukey Bisquare loss function 
#' b1 = 0.5                    # 50% breakdown point
#' cc1 =  1.567                # corresponding model parameter
#' b1 = 0.25; cc1 =  2.937   
#' 
#' # Initialization [PSC analysis for compositional data]
#' control <- robregcc_option(maxiter=1000,tol = 1e-4,lminfac = 1e-7)
#' fit.init <- cpsc_sp(Xt, y,alp = 0.4, cfac = 2, b1 = b1, 
#' cc1 = cc1,C,bw,1,control)  
#' 
#' ## Robust model fitting
#' 
#' # control parameters
#' control <- robregcc_option()
#' beta.wt <- fit.init$betaR    # Set weight for model parameter beta
#' beta.wt[1] <- 0
#' control$gamma = 1            # gamma for constructing  weighted penalty
#' control$spb = 40/p           # fraction of maximum non-zero model parameter beta
#' control$outMiter = 1000      # Outer loop iteration
#' control$inMiter = 3000       # Inner loop iteration
#' control$nlam = 50            # Number of tuning parameter lambda to be explored
#' control$lmaxfac = 1          # Parameter for constructing sequence of lambda
#' control$lminfac = 1e-8       # Parameter for constructing sequence of lambda 
#' control$tol = 1e-20;         # tolrence parameter for converging [inner  loop]
#' control$out.tol = 1e-16      # tolerence parameter for convergence [outer loop]
#' control$kfold = 10           # number of fold of crossvalidation
#' control$sigmafac = 2#1.345
#' # Robust regression using adaptive lasso penalty
#' fit.ada <- robregcc_sp(Xt,y,C,
#'                        beta.init = beta.wt,  cindex = 1, 
#'                        gamma.init = fit.init$residuals,
#'                        control = control, 
#'                        penalty.index = 1, alpha = 0.95)
#' 
#' # Robust regression using lasso penalty [Huber equivalent]   
#' fit.soft <- robregcc_sp(Xt,y,C, cindex = 1, 
#'                         control = control, penalty.index = 2, 
#'                         alpha = 0.95)
#' 
#' 
#' # Robust regression using hard thresholding penalty
#' control$lmaxfac = 1e2               # Parameter for constructing sequence of lambda
#' control$lminfac = 1e-3              # Parameter for constructing sequence of lambda
#' control$sigmafac = 2#1.345
#' fit.hard <- robregcc_sp(Xt,y,C, beta.init = fit.init$betaf, 
#'                         gamma.init = fit.init$residuals,
#'                         cindex = 1, 
#'                         control = control, penalty.index = 3, 
#'                         alpha = 0.95)
#' 
#' coef_cc(fit.ada)
#' coef_cc(fit.soft)
#' coef_cc(fit.hard)
#' 
#' }
coef_cc = function(object, type = 0,  s = 0) {
  if (type == 0) {
    return(object$beta0[, s + 1])
  } else {
    return(object$betaE[, s + 1])
  }
}








#' Extract residuals estimate from the sparse version of the robregcc fitted object. 
#'
#' Robust residuals estimate
#' 
#' @name residuals
#'
#' @param object robregcc fitted onject
#' @param ... Other argumnts for future usage.
#' @return residuals estimate
#' @importFrom stats residuals
#' @rdname residuals
#' @export
#' @examples  
#' 
#' library(magrittr)
#' library(robregcc)
#' 
#' data(simulate_robregcc)
#' X <- simulate_robregcc$X;
#' y <- simulate_robregcc$y
#' C <- simulate_robregcc$C
#' n <- nrow(X); p <- ncol(X); k <-  nrow(C)
#' 
#' Xt <- cbind(1,X)                         # accounting for intercept in predictor
#' C <- cbind(0,C)                           # accounting for intercept in constraint
#' bw <- c(0,rep(1,p))                       # weight matrix to not penalize intercept 
#' 
#' \donttest{
#' example_seed <- 2*p+1               
#' set.seed(example_seed) 
#' 
#' # Breakdown point for tukey Bisquare loss function 
#' b1 = 0.5                    # 50% breakdown point
#' cc1 =  1.567                # corresponding model parameter
#' b1 = 0.25; cc1 =  2.937   
#' 
#' # Initialization [PSC analysis for compositional data]
#' control <- robregcc_option(maxiter=1000,tol = 1e-4,lminfac = 1e-7)
#' fit.init <- cpsc_sp(Xt, y,alp = 0.4, cfac = 2, b1 = b1, 
#' cc1 = cc1,C,bw,1,control)  
#' 
#' ## Robust model fitting
#' 
#' # control parameters
#' control <- robregcc_option()
#' beta.wt <- fit.init$betaR    # Set weight for model parameter beta
#' beta.wt[1] <- 0
#' control$gamma = 1            # gamma for constructing  weighted penalty
#' control$spb = 40/p           # fraction of maximum non-zero model parameter beta
#' control$outMiter = 1000      # Outer loop iteration
#' control$inMiter = 3000       # Inner loop iteration
#' control$nlam = 50            # Number of tuning parameter lambda to be explored
#' control$lmaxfac = 1          # Parameter for constructing sequence of lambda
#' control$lminfac = 1e-8       # Parameter for constructing sequence of lambda 
#' control$tol = 1e-20;         # tolrence parameter for converging [inner  loop]
#' control$out.tol = 1e-16      # tolerence parameter for convergence [outer loop]
#' control$kfold = 10           # number of fold of crossvalidation
#' control$sigmafac = 2#1.345
#' # Robust regression using adaptive lasso penalty
#' fit.ada <- robregcc_sp(Xt,y,C,
#'                        beta.init = beta.wt,  cindex = 1, 
#'                        gamma.init = fit.init$residuals,
#'                        control = control, 
#'                        penalty.index = 1, alpha = 0.95)
#' 
#' # Robust regression using lasso penalty [Huber equivalent]   
#' fit.soft <- robregcc_sp(Xt,y,C, cindex = 1, 
#'                         control = control, penalty.index = 2, 
#'                         alpha = 0.95)
#' 
#' 
#' # Robust regression using hard thresholding penalty
#' control$lmaxfac = 1e2               # Parameter for constructing sequence of lambda
#' control$lminfac = 1e-3              # Parameter for constructing sequence of lambda
#' control$sigmafac = 2#1.345
#' fit.hard <- robregcc_sp(Xt,y,C, beta.init = fit.init$betaf, 
#'                         gamma.init = fit.init$residuals,
#'                         cindex = 1, 
#'                         control = control, penalty.index = 3, 
#'                         alpha = 0.95)
#'                         
#'                         
#' residuals(fit.ada)
#' residuals(fit.soft)
#' residuals(fit.hard)
#' 
#' 
#' }
residuals.robregcc <- function(object, ... ) {
  type = 0; s = 1
  if (type == 0) {
    return(object$residuals0[, s + 1])
  } else {
    return(object$residualsE[, s + 1])
  }
}