SFSI (version 1.2.0)

solveEN: Coordinate Descent algorithm to solve Elastic-Net-type problems

Description

Computes the entire Elastic-Net solution for the regression coefficients for all values of the penalization parameter, via the Coordinate Descent (CD) algorithm (Friedman, 2007). It uses as inputs a variance matrix among predictors and a covariance vector between response and predictors

Usage

solveEN(Sigma, Gamma, X = NULL, alpha = 1, lambda = NULL,
        nlambda = 100, common.lambda = TRUE, 
        lambda.min = .Machine$double.eps^0.5, dfmax = NULL, 
        scale = TRUE, tol = 1E-5, maxiter = 1000, mc.cores = 1L,
        return.beta = TRUE, save.beta = FALSE, verbose = FALSE)

Value

Returns a list object containing the elements:

• lambda: (vector) all the sequence of values of the penalty.

• beta: (matrix) regression coefficients for each predictor (in rows) associated to each value of the penalization parameter lambda (in columns).

• df: (vector) degrees of freedom, number of non-zero predictors associated to each value of lambda.

• yHat: (matrix) fitted values calculated using a matrix of predictors (when argument X is not NULL), associated to each value of lambda (in columns).

The returned object is of the class 'LASSO' for which methods coef and fitted exist. Function 'path.plot' can be also used

Arguments

Sigma

(numeric matrix) Variance-covariance matrix of predictors. It can be of the "float32" type as per the 'float' R-package

Gamma

(numeric matrix) Covariance between response variable and predictors. If it contains more than one column, the algorithm is applied to each column separately as different response variables

X

(numeric matrix) Optional matrix of predictors to obtain fitted values

lambda

(numeric vector) Penalization parameter sequence. Default is lambda=NULL, in this case a decreasing grid of 'nlambda' lambdas will be generated starting from a maximum equal to

max(abs(Gamma)/alpha)

to a minimum equal to zero. If alpha=0 the grid is generated starting from a maximum equal to 5

nlambda

(integer) Number of lambdas generated when lambda=NULL

lambda.min

(numeric) Minimum value of lambda that are generated when lambda=NULL

common.lambda

TRUE or FALSE to whether computing the coefficients for a grid of lambdas common to all columns of Gamma or for a grid of lambdas specific to each column of Gamma. Default is common.lambda=TRUE

alpha

(numeric) Value between 0 and 1 for the weights given to the L1 and L2-penalties

scale

TRUE or FALSE to scale matrix Sigma for variables with unit variance and scale Gamma by the standard deviation of the corresponding predictor taken from the diagonal of Sigma

tol

(numeric) Maximum error between two consecutive solutions of the CD algorithm to declare convergence

maxiter

(integer) Maximum number of iterations to run the CD algorithm at each lambda step before convergence is reached

dfmax

(integer) Maximum number of non-zero coefficients in the last solution. Default dfmax=NULL will calculate solutions for the entire lambda grid

mc.cores

(integer) Number of cores used. The analysis is run in parallel when mc.cores is greater than 1. Default is mc.cores=1

return.beta

TRUE or FALSE to whether return regression coefficients in the output object

save.beta

TRUE or FALSE to whether save regression coefficients (in a temporary folder). When TRUE coefficients are not returned in the output object and instead the path where coefficients were saved is returned. They can be retrieved using coef method if at least one return.beta or save.beta is TRUE

verbose

TRUE or FALSE to whether printing progress

Author

Marco Lopez-Cruz () and Gustavo de los Campos

Details

Finds solutions for the regression coefficients in a linear model

y_i = x'_i β + e_i

where y_i is the response for the i^th observation, x_i=(x_i1,...,x_ip)' is a vector of \(p\) predictors assumed to have unit variance, β=(β₁,...,β_p)' is a vector of regression coefficients, and e_i is a residual.

The regression coefficients β are estimated as function of the variance matrix among predictors (Σ) and the covariance vector between response and predictors (Γ) by minimizing the penalized mean squared error function

-Γ' β + 1/2 β' Σ β + λ J(β)

where λ is the penalization parameter and J(β) is a penalty function given by

1/2(1-α)||β||₂² + α||β||₁

where 0 ≤ α ≤ 1, and ||β||₁ = ∑_j=1|β_j| and ||β||₂² = ∑_j=1β_j² are the L1 and (squared) L2-norms, respectively.

The "partial residual" excluding the contribution of the predictor x_ij is

e_i^(j) = y_i - x'_i β + x_ijβ_j

then the ordinary least-squares (OLS) coefficient of x_ij on this residual is (up-to a constant)

β_j^(ols) = Γ_j - Σ'_j β + β_j

where Γ_j is the j^th element of Γ and Σ_j is the j^th column of the matrix Σ.

Coefficients are updated for each \(j=1,...,p\) from their current value β_j to a new value β_j(α,λ), given α and λ, by "soft-thresholding" their OLS estimate until convergence as fully described in Friedman (2007).

References

Friedman J, Hastie T, Höfling H, Tibshirani R (2007). Pathwise coordinate optimization. The Annals of Applied Statistics, 1(2), 302–332.

Hoerl AE, Kennard RW (1970). Ridge Regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67.

Tibshirani R (1996). Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society B, 58(1), 267–288.

Zou H, Hastie T (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society B, 67(2), 301–320.

Examples

  require(SFSI)
  data(wheatHTP)
  
  y = as.vector(Y[,"E1"])  # Response variable
  X = scale(X_E1)          # Predictors

  # Training and testing sets
  tst = which(Y$trial %in% 1:10)
  trn = seq_along(y)[-tst]

  # Calculate covariances in training set
  XtX = var(X[trn,])
  Xty = cov(X[trn,],y[trn])
  
  # Run the penalized regression
  fm1 = solveEN(XtX,Xty,alpha=0.5,nlambda=100) 
  
  # Predicted values
  yHat1 = fitted(fm1, X=X[trn,])  # training data
  yHat2 = fitted(fm1, X=X[tst,])  # testing data
  
  # Penalization vs correlation
  plot(-log(fm1$lambda[-1]),cor(y[trn],yHat1[,-1]), main="training")
  plot(-log(fm1$lambda[-1]),cor(y[tst],yHat2[,-1]), main="testing")
  # \donttest{
  if(requireNamespace("float")){
   # Using a 'float' type variable
   XtX2 = float::fl(XtX)
   fm2 = solveEN(XtX2,Xty,alpha=0.5)  
   max(abs(fm1$beta-fm2$beta))   # Check for discrepances in beta
  }
  # }