## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(sglssnal)

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

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

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

## -----------------------------------------------------------------------------
Anew <- matrix(rnorm(5 * p), 5, p)
predict(fit, Anew)

## -----------------------------------------------------------------------------
dim(riboflavin$A)
length(unique(riboflavin$group))

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

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

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