This package provides tools for fitting regularization paths for
sparse group-lasso penalized learning problems. The model is fit for a
sequence of the regularization parameter lambda.
The strengths and improvements that this package achieved compared to
other sparse group-lasso packages are as follows:
Installing
This package is not on CRAN yet, so it can be installed using the
[devtools] (https://cran.r-project.org/package=devtools)
package:
devtools::install_github("dajmcdon/sparsegl", ref = "main")
Building the vignettes, such as this getting started guide, takes a
significant amount of time. They are not included in the package by
default. If you want to include vignettes, then use this modified
command:
devtools::install_github("dajmcdon/sparsegl", ref = "main",
build_vignettes = TRUE, dependencies = TRUE)
For this getting-started vignette, firstly, we will randomly generate
X, an input matrix of predictors of dimension
n-obs by p-feature. To initiate
y, a real-valued vector (make it a vector rather a
matrix),
- Linear Regression model: \(y = X\beta^* +
\epsilon\).
- Logistic regression model: \(y = (y_1,
y_2, \cdots, y_n)\), where \(y_i \sim
\text{Bernoulli}\left(\frac{1}{1 + \exp(-X_i \beta^*)}\right)\),
\(i = 1, 2, \cdots, n.\)
where the coefficient vector \(\beta^*\) is specified as below, and the
white noise \(\epsilon\) following
standard normal distribution serves as data variation. Then the sparse
group-lasso problem is formulated as the sum of mean squared error (
linear regression) or logistic loss (logistic regression) and a convex
combination of lasso penalty and group lasso penalty:
- Linear regression: \[
\min_{\beta\in\mathbb{R}^p}\left(\frac{1}{2n} \rVert y - \sum_g
X^{(g)}\beta^{(g)}\rVert_2^2 + (1-\alpha)\lambda\sum_g
\sqrt{\text{bs}}\rVert\beta^{(g)}\rVert_2 +
\alpha\lambda\rVert\beta\rVert_1 \right) \qquad (*).
\]
- Logistic regression: \[
\min_{\beta\in\mathbb{R}^p}\left(\frac{1}{2n}\sum_{i=1}^n \log\left(1 +
\exp\left(-y_i\beta^\top X_i\right)\right) + (1-\alpha)\lambda\sum_g
\sqrt{\text{bs}}\rVert\beta^{(g)}\rVert_2 +
\alpha\lambda\rVert\beta\rVert_1 \right) \qquad (**).
\]
where
\(X^{(g)}\) is the submatrix of
\(X\) with columns corresponding to the
features in group \(g\).
\(\beta^{(g)}\) is the
corresponding coefficients of the features in group \(g\).
bs is the length of \(\beta^{(g)}\).
\(\alpha\) adjusts the weight
between lasso penalty and group-lasso penalty.
\(\lambda\) fine-tunes the size
of penalty imposed on the model to control the number of nonzero
coefficients of features, and furthermore, to avoid
overfitting.
library(sparsegl)
set.seed(1010)
n <- 100
p <- 200
X <- matrix(data = rnorm(n*p, mean = 0, sd = 1), nrow = n, ncol = p)
beta_star <- c(rep(5, 5), c(5, -5, 2, 0, 0), rep(-5, 5), c(2, -3, 8, 0, 0), rep(0, (p - 20)))
groups <- rep(1:(p / 5), each = 5)
# Linear regression model
eps <- rnorm(n, mean = 0, sd = 1)
y <- X %*% beta_star + eps
# Logistic regression model
pr <- 1 / (1 + exp(-X %*% beta_star))
y_binary <- rbinom(n, 1, pr)
sparsegl()
Given an input matrix X, and a response vector
y (or a matrix with 1 column), a sparse group-lasso
regularized linear model is fitted for a sequence of penalty parameter
values in terms of penalized maximum likelihood. The penalty is composed
of lasso penalty and group lasso penalty. The other main arguments the
users might give are:
group: a vector with consecutive integers of length
p indicating the grouping of the features. By default, each
group only contains one feature if without initialization.
family: A character string specifying the likelihood
to use, could be either linear regression “gaussian” or logistic
regression loss “binomial”. Default is “gaussian”.
pf: a penalty vector of the length of the number of
groups in total. Default value for each entry is the square-root of the
corresponding size of each group, which is \(\sqrt{\text{bs}}\) in \((*)\) and \((**)\) above.
asparse: changes the weight of lasso penalty,
referring to \(\alpha\) in \((*)\) and \((**)\) above: asparse = \(1\) gives the lasso penalty only.
asparse = \(0\) gives the
group lasso penalty only. The default value of asparse is
\(0.05\).
lower_bnd: lower bound for coefficient values, a
vector in length of 1 or the number of groups including non-positive
numbers only. Default value for each entry is -\(\infty\).
upper_bnd: upper bound for coefficient values, a
vector in length of 1 or the number of groups including non-negative
numbers only. Default value for each entry is \(\infty\).
It returns a sparsegl object, and the main attribute of
this object is:
df: the number of nonzero coefficients for each lambda
value. This only represents an approximation to the exact
degree-of-freedom. Detailed explanation can be checked under function
estimate_risk() below.
fit1 <- sparsegl(X, y, group = groups)
Plotting function plot() for sparsegl
object
This function produces nonzero-only coefficient curves for each
penalty parameter lambda values in the regularization path
for a fitted sparsegl object. The arguments of this
function are:
To elaborate on y_axis and x_axis:
The plot with y_axis = “group” is the group norms
against the log-lambda or the scaled group norm vector.
Each group norm is defined by: \[
\text{asparse}\times\rVert\beta^{(g)}\rVert_1 + (1 -
\text{asparse})\times\sum_g\rVert\beta^{(g)}\rVert_2
\] Curves are plotted in the same color if the corresponding
features are in the same group. Notice that the number of curves shown
on the plots could be less than the actual number of groups since only
the groups containing nonzero features at least at one lambda in the
sequence lambda are included.
The plot with y_axis = “coef” is the coefficients of
features against the log-lambda or the scaled group norm
vector. Each curve with a distinct color represents a feature. Again,
only the features having nonzero coefficients at least at one lambda
value in the sequence lambda are considered.
The plot with x_axis = “lambda” indicates the
x_axis presenting with log-lambda.
The plot with x_axis = “penalty” indicates the
x_axis presenting with the scaled group norm vector. Each
element in this vector is defined by: \[
\frac{\text{asparse}\times\rVert \beta\rVert_1 +
(1-\text{asparse})\times\sum_g\rVert
\beta^{(g)}\rVert_2}{\max_\beta\left(\text{asparse}\times\ \rVert
\beta\rVert_1 + (1-\text{asparse})\times\sum_g\rVert
\beta^{(g)}\rVert_2\right)}
\]
plot(fit1, y_axis = "group", x_axis = "lambda")
plot(fit1, y_axis = "coef", x_axis = "penalty", add_legend = FALSE)
coef(), predict() and print()
sparsegl object
All three functions consume a fitted sparsegl object as
an arguments and
coef() and predict() return a matrix of
coefficients and predictions \(\hat{y}\) given a matrix X at
each lambda respectively, unless an optional argument s is
assigned with a specified vector or numeric value from (but not
necessarily) the original lambda sequence.
print() returns the number of nonzero features with
nonzero coefficients at each lambda.
coef <- coef(fit1, s = c(0.02, 0.03))
pred <- predict(fit1, newx = X[100,], s = fit1$lambda[2:3])
print(fit1)
cv.sparsegl()
This function does a k-fold cross-validation (cv) on
sparsegl. It takes the same arguments X,
y, group, which are specified above, with
additional argument pred.loss. It can be set with either
“L2” or “L1” for linear regression model, “loss” or “misclass” for
logistic regression model indicating the loss to use for cv error. This
will return a cv.sparseg object.
fit_l1 <- cv.sparsegl(X, y, group = groups, pred.loss = "L1")
plot(), coef() and predict()
for cv.sparsegl object
plot(): produces cross-validation curve with upper and
lower confidence bounds plots for each lambda in the
regularization path for a fitted cv.sparsegl object. In the
generated plot, the red points represent the cvm values at
each lambda, and the attached gray vertical error bars
indicate the uncertainty of cvm.
coef() and predict() for
cv.sparsegl object work similarly to those introduced
above. The only differences are that the argument object
passes a fitted cv.sparsegl object, and s can
also be set with lambda.1se, lambda.min in
addition.
coef <- coef(fit_l1, s = "lambda.1se")
pred <- predict(fit_l1, newx = X[50:80, ], s = "lambda.min")
estimate_risk()
This function returns the information criterion, which is the sum of
the maximum log-likelihood and a penalty term determined by the chosen
penalty type for a sparsegl model at each
lambda. It provides a means for model selection by
representing the trade-off between the goodness of fit of the model and
the complexity of the model. It takes the same arguments X
and y from the function sparsegl(). The
additional arguments it needs are:
where df is the degree-of-freedom, and n is the sample size.
approx_df: indicate if a vector of the approximation to
the degree-of-freedom at each penalty parameter lambda is
used. Default is FALSE and the program will compute the unbiased
estimate of the exact degree-of-freedom.
FYI: Degree-of-freedom is a tool to assess the complexity of a
statistical modeling procedure. object$df, the
approximation to the degree-of-freedom is the number of nonzero
coefficients of the model. Notice that it would take some time to
calculate the unbiased estimate of the exact degrees-of-freedom if
X is complicated. For more details about how to realize
this calculation, our method is implemented based on the paper https://arxiv.org/pdf/1212.6478.pdf.
risk <- estimate_risk(fit1, X, type = "AIC", approx_df = FALSE)