---
title: "Introduction to sglssnal"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to sglssnal}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(sglssnal)
```

## 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

```{r}
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)
```

`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`:

```{r}
fit_path <- sglssnal(A, b, group, nlambda = 10, alpha = 0.5, verbose = 0)
dim(coef(fit_path)) # one column per lambda
fit_path$lambda
```

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:

```{r}
cvfit <- cv.sglssnal(A, b, group, nlambda = 10, alpha = 0.5, nfolds = 5, verbose = 0)
cvfit$cv_info$cv_lambda_id
coef(cvfit)[, cvfit$cv_info$cv_lambda_id]
```

`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.

## predict()

```{r}
Anew <- matrix(rnorm(5 * p), 5, p)
predict(fit, Anew)
```

## 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):

```{r}
dim(riboflavin$A)
length(unique(riboflavin$group))
```

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.

```{r}
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]
```

The cross-validated error along the lambda path is shown below, with the selected optimum marked in red.

```{r, fig.width=6, fig.height=4}
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:

```{r}
beta <- coef(cv_ribo)[-1, best] # drop the intercept row
active <- unique(riboflavin$group[beta != 0])
length(active) # of 36
sum(beta != 0) # of 1199
table(riboflavin$group[beta != 0]) # nonzero genes per active term
```

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`).