What this package solves
sglssnal fits the sparse-group lasso
\[\min_x \tfrac12\|Ax-b\|_2^2 +
\alpha\lambda\|x\|_1 + (1-\alpha)\lambda\sum_i
w_i\|x_{G_i}\|_2\]
where predictor columns are partitioned into groups \(G_i\). The \(\ell_1\) term drives ordinary
coordinate-level sparsity; the group \(\ell_2\) term drives whole-group
sparsity, dropping entire groups to zero at once. alpha
controls the mix between the two.
Unlike most lasso solvers (Simon et al. 2013, Liang et al. 2024)
which descend directly on the primal via coordinate descent or proximal
gradient, sglssnal applies a second-order semismooth Newton
method to the problem’s dual – the algorithm of Zhang, Zhang,
Sun & Toh (2020). This tends to matter most exactly when a
lasso-type solver is slowest: at small lambda,
with correlated predictors, or where first-order methods take many
iterations to fully sparsify.
A minimal fit
set.seed(1)
n <- 50
p <- 20
A <- matrix(rnorm(n * p), n, p)
bstar <- c(2, -3, rep(0, p - 2)) # only the first group is truly active
b <- as.numeric(A %*% bstar + rnorm(n, sd = 0.1))
group <- rep(1:4, each = 5) # 4 groups of 5 columns each
fit <- sglssnal(A, b, group, lambda = 0.3, alpha = 0.5, verbose = 0)
coef(fit)
#> 21 x 1 sparse Matrix of class "dgCMatrix"
#> s0
#> (Intercept) -0.12645053
#> V1 1.96560174
#> V2 -2.92455897
#> V3 -0.01738276
#> V4 .
#> V5 .
#> V6 .
#> V7 .
#> V8 .
#> V9 .
#> V10 .
#> V11 .
#> V12 .
#> V13 .
#> V14 .
#> V15 .
#> V16 .
#> V17 .
#> V18 .
#> V19 .
#> V20 .
coef() returns a sparse matrix with an intercept row
prepended. Here group 1 (columns 1-5, which contain the only two nonzero
true coefficients) survives; groups 2-4 are zeroed out entirely.
Fitting a lambda path
A path of lambda may be provided. By default,
sglssnal() fits a full descending path with
nlambda values, auto-generated down to
lambda_min_ratio of the largest, reusing each fit’s
solution as the warm start for the next, smaller
lambda:
fit_path <- sglssnal(A, b, group, nlambda = 10, alpha = 0.5, verbose = 0)
dim(coef(fit_path)) # one column per lambda
#> [1] 21 10
fit_path$lambda
#> [1] 41.389136946 14.874484589 5.345612596 1.921113559 0.690412416
#> [6] 0.248121357 0.089170192 0.032046106 0.011516773 0.004138914
Both sglssnal() and cv.sglssnal() default
alpha to \(0.05\) (mostly
group penalty), matching Simon et al. (2013)’s own real-data examples
and the sparsegl package’s asparse default
(Liang et al. 2024).
Cross-validation
cv.sglssnal() fits the same kind of path but chooses
lambda by \(k\)-fold
cross-validated prediction error:
cvfit <- cv.sglssnal(A, b, group, nlambda = 10, alpha = 0.5, nfolds = 5, verbose = 0)
cvfit$cv_info$cv_lambda_id
#> [1] 6
coef(cvfit)[, cvfit$cv_info$cv_lambda_id]
#> (Intercept) V1 V2 V3 V4
#> -1.534080e-01 1.974676e+00 -2.935600e+00 -2.168284e-02 1.254773e-05
#> V5 V6 V7 V8 V9
#> 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> V10 V11 V12 V13 V14
#> 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> V15 V16 V17 V18 V19
#> 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#> V20
#> 0.000000e+00
cv.sglssnal() fits the full dataset first to establish
the lambda path – that fit (not a fold-average) is what’s returned, with
cv_info (the lambda sequence, cross-validated error
cvm, and the selected index) attached.
A real worked example: predicting riboflavin yield from gene
expression
The bundled riboflavin dataset comes from DSM’s
industrial strain-improvement program for riboflavin (vitamin B2)
production via fermentation with engineered Bacillus subtilis:
71 strain variants, each profiled for the expression of 4088 genes, with
the production yield (log-transformed) as the outcome. See
?riboflavin for the full source.
1199 of those genes have a Biological Process annotation mapping to a
term in the “generic GO Slim” (a curated, coarse-grained subset of the
Gene Ontology; CC BY 4.0) and are included here, grouped into 36 terms
(sizes ranging from 284 genes down to singletons):
dim(riboflavin$A)
#> [1] 71 1199
length(unique(riboflavin$group))
#> [1] 36
13% of these 1199 genes have annotations spanning more than one Slim
term. sglssnal’s group argument requires a
strict partition, so each such gene here was assigned to its
most-frequently annotated term.
cv_ribo <- cv.sglssnal(riboflavin$A, riboflavin$b, riboflavin$group,
nlambda = 20, lambda_min_ratio = 1e-3, alpha = 0.75, nfolds = 5, verbose = 0
)
best <- cv_ribo$cv_info$cv_lambda_id
cv_ribo$cv_info$lambda[best]
#> [1] 0.02092215
The cross-validated error along the lambda path is shown below, with
the selected optimum marked in red.
plot(cv_ribo$cv_info$lambda, cv_ribo$cv_info$cvm,
log = "x", type = "b", pch = 16,
xlab = "lambda", ylab = "cross-validated error",
main = "riboflavin: CV error along the lambda path"
)
abline(v = cv_ribo$cv_info$lambda[best], lty = 2, col = "red")
At that optimum, a minority of the 36 terms survive, at varying
density:
beta <- coef(cv_ribo)[-1, best] # drop the intercept row
active <- unique(riboflavin$group[beta != 0])
length(active) # of 36
#> [1] 10
sum(beta != 0) # of 1199
#> [1] 136
table(riboflavin$group[beta != 0]) # nonzero genes per active term
#>
#> GO:0005975 GO:0006281 GO:0006520 GO:0006766 GO:0007059 GO:0042254 GO:0055085
#> 11 8 24 9 1 3 43
#> GO:0055086 GO:0071941 GO:1901135
#> 23 3 11
The selected terms include the vitamin metabolic process itself
(GO:0006766 – riboflavin’s own biosynthesis pathway), its
direct precursor supply (GO:0055086 nucleobase-containing
small molecule metabolism; GTP is riboflavin’s biosynthetic precursor),
energy and carbon metabolism (GO:0005975,
GO:1901135), amino acid metabolism
(GO:0006520, needed to build the biosynthetic enzymes),
transmembrane transport (GO:0055085), and machinery for
growth and division (GO:0006281 DNA repair,
GO:0007059 chromosome segregation, GO:0042254
ribosome biogenesis) and nitrogen assimilation
(GO:0071941).