Tunes the regularisation parameter lambda of an L0glm fit

Tunes the lambda parameter for the L0glm.fit function. It uses a sequence of penalty parameter values (lambda) to fit L0-penalized GLMs, and selects the best lambda by an information criterion or held-out deviance from cross-validation. It can handle multiple output variables (channels), and can produce a plots of the regularization path, true vs. estimated coefficients (if true coefficients are passed, e.g. from a simulation), estimated support size, deviance, log likelihood, mean absolute deviation, accuracy, precision, recall and the IC of interest, all in function of the level of regularization (L0 penalty lambda). Some plots are only produced if the user passes the true coefficients coefs.true (e.g. from a simulation). Each lambda fit is initialized independently. The preset fit is used only to construct the grid; its selected support is not imposed on other fits. A common supplied start or pilot-derived seed can be reused independently at every lambda, avoiding dependence on the order of the grid. The preset calibrates adaptive ridge to the selected information criterion. Path tuning can improve its fit for correlated predictors or unrelaxed estimates that retain shrinkage.

Usage

L0glm.tune(
  X,
  y,
  family,
  weights = matrix(1, nrow = nrow(as.matrix(y)), ncol = ncol(as.matrix(y))),
  offset = Matrix::Matrix(0, nrow = nrow(X), ncol = ncol(as.matrix(y)), sparse = TRUE),
  dispersion = NULL,
  tune = c("gic", "aic", "bic", "ebic", "mbic", "cv"),
  nlambda = 10,
  prop_range = 1,
  coefs.true = NULL,
  printbest = TRUE,
  plot = TRUE,
  mu_floor = 1e-08,
  nfolds = 5,
  foldid = NULL,
  lambda_min = 2,
  start = NULL,
  ...
)

Arguments

X

A matrix of predictor variables.

y

Response data with n observations, aligned with the rows of X. For ordinary GLMs, a vector of length n or an n by m matrix of response variables (tasks/channels). Special-family response formats are:

mgcv::nb(): a vector or n by m matrix of nonnegative integer counts, with one column per response channel.
mgcv::multinom(K): a factor/class vector of length n with K+1 categories (one trial per observation), or an n by K+1 nonnegative integer count matrix with one column per category. Row sums give the numbers of trials.
mgcv::cox.ph(): a right-censored survival::Surv(time, status) response, or a numeric two-column matrix with time in the first column and status in the second. Use status=1 for an observed event and status=0 for right censoring. Event indicators belong in y, not in weights.

family

Model family and link function. All standard GLM families and their supported link functions are accepted; see family. Special families are mgcv::nb() for negative binomial regression with fixed or estimated shape, mgcv::multinom(K) for multinomial regression with K+1 categories, and mgcv::cox.ph() for right-censored Cox proportional hazards models. See Special families under Details. Fixed-noise scalar models nb1(), Gamma_scale() and Gamma_censored() use the ordinary single- or multichannel response formats with exact likelihood-based IRLS; see signal_families.

weights

Optional finite, nonnegative prior/case observation weights. For ordinary GLMs, supply a vector of length n or an n by m matrix; for binomial proportions, weights give trial totals. Special families accept a scalar, a vector of length n, or an n by 1 matrix; NULL uses unit weights.

mgcv::nb(): these formats share weights across response channels; an n by m matrix supplies channel-specific weights.
mgcv::multinom(K): one weight per observation, shared across categories. For a count-matrix response, trial totals come from row sums of y; weights multiply each observation's log-likelihood contribution.
mgcv::cox.ph(): case weights enter both weighted event counts and risk sets. Event indicators belong in y, not in weights.

At least one positive weight is required.

offset

Optional finite, known additive offsets on the linear-predictor (link) scale. For ordinary GLMs, supply a vector of length n or an n by m matrix. For special families, NULL uses zero offsets; a scalar is replicated over observations and response columns, while a length-n vector or n by 1 matrix supplies one offset per observation, shared across response columns.

mgcv::nb(): an n by m matrix supplies channel-specific offsets on the chosen link scale.
mgcv::multinom(K): an n by K+1 matrix supplies category-specific offsets, including every category. Offsets are centered across categories within each row. A scalar or shared observation offset therefore cancels from the multinomial probabilities.
mgcv::cox.ph(): use a scalar, length-n vector or n by 1 matrix of offsets on the log-relative-hazard scale.

dispersion

Passed to L0glm.fit: NULL automatically initializes and estimates dispersion, a positive value fixes it, and a negative value supplies its absolute value as an initial estimate. A vector can specify separate values for the response channels. Fixed-scale families retain their family dispersion.

tune

A case-insensitive method for model selection: "aic", "bic", "gic", "ebic", "mbic", or "cv". Defaults to "gic". Cross-validation minimizes mean held-out deviance. See information_criteria for formulas, guidance and references.

nlambda

Number of lambda values to use for tuning (at least 2; default 10). With at least three points and an interior reference, approximately half lie on either side of lambda_reference, which is included exactly. A two-point grid contains only the displayed endpoints. If the reference equals the lower bound, all points lie at or above it.

prop_range

Numeric in $(0,\,1]$. Fraction of each log10 distance from lambda_reference to the full endpoints lambda_min and lambda_max retained in the grid. prop_range = 1 keeps the full range; smaller values zoom toward the reference, with each endpoint moving inward on the log scale.

coefs.true

A vector or matrix of true coefficient values, for comparison with the estimated coefficients.

printbest

A logical value indicating whether to print the best lambda value according to the selected criterion or CV (default is TRUE).

plot

A logical value indicating whether to produce the plots mentioned in the description. Both panels use $\log_{10}(\lambda)$ on a reversed horizontal axis, with weaker regularization to the right. Coefficient paths have a separate color per predictor and are labeled $\hat\beta$. Multichannel paths show coefficient row means. The top axes show support sizes at 16 equally spaced log-scale positions, using the nearest fitted lambda at each position. For multichannel fits these count rows with any nonzero coefficient, independently of the plotted row means. Paths and support counts exclude mandatory terms; no additional models are fitted.

mu_floor

Nonnegative prediction floor used only when scoring positive-mean families, such as Poisson (default 1e-8). It retains all observations and makes scores finite when a fitted mean is zero. Set to zero for exact likelihood scoring, which can give infinite criteria.

nfolds

Number of cross-validation folds (default 5), used when tune = "cv". Must be at least 2 and no larger than the number of positive-weight observation rows.

foldid

Optional integer fold assignment for each row of X. Labels must run from 1 to the number of folds. When supplied without nfolds, the number of folds is inferred. All response channels of a row are held out together. Otherwise balanced random folds are generated; use set.seed() for reproducibility.

lambda_min

Finite positive lower endpoint of the full-data lambda grid (default 2, the AIC preset). Set a smaller positive value to explore weaker regularization or increase it to focus on sparser models. CV training penalties retain their training-weight scaling, so they can be smaller than this full-data bound. A zoomed grid can start above it.

start

NULL (default) automatically initializes each lambda fit. Numeric starts can be a scalar, a length-p vector repeated across response columns, or a p by m matrix, including any added intercept row. Negative binomial models use one column per response channel; Cox uses one column and no intercept; mgcv::multinom(K) uses all K+1 category columns in response-category order, centered rowwise to CLR. Numeric scalar 1 supplies ones in each fit's normalized working coordinates; other numeric starts use original coefficient units. For ordinary GLMs, negative binomial and Cox models, exact starting zeros exclude individual coefficients. For multinomial models, rows supplied as all zero are excluded before CLR centering; zero category entries in other rows remain eligible. Scalar 1 centers to zero without excluding rows, whereas scalar 0 excludes all selectable rows. Intercepts and always.include terms remain eligible. These restrictions also apply during relaxed refits, and output coefficients retain their original sparse dimensions. See L0glm.fit for further starting-value details. Use "pilot" to build a common seed from ten independent unrelaxed fits at log-spaced penalties between lambda_min and the GIC preset (or one fit when these bounds coincide). The magnitude is the mean absolute coefficient across pilots, with the sign of its mean; exact cancellation uses the sign of the first largest-magnitude pilot coefficient. For ordinary GLMs, negative binomial and Cox models this retains the union of pilot supports separately in each channel. Multinomial seeds are centered to CLR, then rows that are all zero in the pooled seed are excluded. CV constructs each pilot seed using training observations only. No preceding path fit is used to initialize the next fit.

...

Other arguments passed to L0glm.fit, including:

screen, nsis, screen.maxit, screen.nthreads: the screening strategy, candidate batch size, outer iteration limit and OpenMP screening thread count from L0glm.fit(). The default screen="auto" enables ISIS exactly when $p_{\mathrm{eligible}}>10n$. Here $n$ is the number of training rows in the current fit, and the eligible count excludes mandatory terms and predictors excluded by explicit zero starts in every response channel. No separate minimum sample size or predictor count is imposed. Screening is recomputed inside each full-data or CV training fit; validation responses never determine its candidates. screen="none" disables it.
precluster, cluster.refit, control.precluster: optional Poisson channel aggregation followed by original-channel refinement; see preclustering. Tuning requires cluster.refit=TRUE. A partition is computed once for the full-data path and separately for each CV training fold. Only memberships are shared across lambda fits; candidate supports and coefficients are fitted independently at every lambda.
always.include: unpenalized columns, including optional smooth background bases. Nonnegative fits refine these in a constrained conditional block independently within each training fit; see compulsory_terms.
grouped: TRUE (default) selects predictor rows across response channels or multinomial categories; FALSE selects individual channel coefficients or CLR category entries. Multinomial fitting always uses the joint likelihood, and Cox has one coefficient per predictor so the choices coincide there.
relaxed: refit the selected support before scoring each candidate full-data or CV training model. Special families use lambda = 0 for this refit. Grouped refits admit eligible entries within retained rows; ungrouped refits retain only selected entries. The default is TRUE for grouped multichannel ordinary GLMs and NB, and grouped multinomial models with either factor or count-matrix responses; it is FALSE for Cox and ungrouped or other single-channel fits. Set it explicitly to refit an ungrouped model. With screening enabled the default is TRUE, as in L0glm.fit().
thresh: relative hard cutoff during selection (default 1e-5), applied to working coefficient row norms for grouped multichannel NB/multinomial fits and individual magnitudes for ungrouped or single-channel fits. The absolute 1e-20 cutoff remains when thresh = 0. Multinomial rows are restored to zero-sum CLR coordinates after thresholding. Mandatory terms are exempt, and relaxed refits do not apply these cutoffs. Cutoffs use each fit's normalized working coefficients, not returned original-unit coefficients, and do not smooth adaptive weights.
prior.mean: ridge-penalty centers in original input coefficient units.

See Special families below and L0glm.fit for support restrictions and further details of these arguments.

Value

L0glm.tune returns a list containing the following elements:

tune: chosen selection method.
p: original number of selectable predictors, retained despite structural starting zeros and pilot support restriction.
start: the common numeric start, including any intercept row, or NULL for independent automatic starts.
pilot: pilot penalties, elapsed seconds and row support sizes, or NULL when no pilot was requested.
fit_args: fitting arguments, including the resolved common start, excluding X, y and lambda. For example, do.call(L0glm.fit, c(list(X = X, y = y, lambda = result$lambda_best), result$fit_args)) reproduces fit_best with the same package version.
preprocessing: training preprocessing of the full-data fit.
metrics: a data frame containing metrics evaluating the models fit with different regularization parameters; CV adds cv_deviance and cv_se.
lambda_best: best lambda by the chosen selection method.
lambda_1se: largest evaluated lambda within one SE of the minimum CV mean deviance, using the SE at that minimum; NULL for information-criterion tuning.
lambda_preset: preset lambda value that a priori would be expected to approximately optimise the chosen criterion, or GIC when tuning by CV.
lambda_reference: positive grid reference; normally identical to lambda_preset, but at least lambda_min.
lambdas: a vector containing all the lambdas that were used.
lambda_min: lower endpoint of the full numerical path.
lambda_max: upper endpoint of the full numerical path.
lambda_nonzero: nonzero lower bracket for the refined endpoint, or NA when the preset is already zero or refinement is unavailable.
endpoint_zero: whether the fitted upper endpoint has zero penalized coefficients.
mu_floor: prediction floor used for scoring.
coefficients: list of coefficient matrices in lambda order.
fit_best: the fitted model at lambda_best.
fit_1se: the full-data fitted model at lambda_1se, or NULL for information-criterion tuning.
cv: CV results, or NULL for criterion tuning. Contains lambda, mean, se, the lambda-by-fold deviance matrix, foldid, nfolds, fold_weights, and the training penalties lambda_train. fit_best, fit_1se and coefficients are full-data fits, not individual fold models. Also contains each fold's starts, pilots and training preprocessing.

When channel preclustering is requested, preclustering retains the full-data partition and cut diagnostics, and cv$preclustering retains the corresponding training-only diagnostics for every fold.

Details

For prior.mean = NULL or an explicitly all-zero prior mean, the upper endpoint is expanded until an actual fit has no nonzero penalized coefficients, excluding an intercept and always.include terms. When the preset fit is nonzero, logarithmic bisection refines the nonzero/zero bracket to within 5%, so a conservative working-model estimate does not leave a long flat tail in the plot. This endpoint is determined by the fitted path. For example, an all-zero Poisson identity-link model with zero offset has infinite deviance when any observed count is positive. The default scoring floor provides numerical comparisons along the path; it does not change the fitted model. With a nonzero prior mean or an unsuccessful endpoint search, a warning is issued and endpoint_zero is false. A zoomed grid need not reach the full endpoint. Supplied observation weights, offsets and known dispersions are included in scoring; unknown Gaussian scales are profiled separately for each response channel. Quasi-Poisson information criteria use a common dispersion across the path. When unspecified, it is estimated by a preset pilot fit and then held fixed; estimated scales retain their nuisance-parameter count. With CV, dispersion estimation uses training observations independently in each fold.

With tune = "cv", the full-data GIC preset defines the reference grid point. If that preset is nonpositive, the grid uses AIC's positive reference value 2 while lambda_preset still reports GIC. The reference is raised to lambda_min if necessary; the reported preset is unchanged. The common grid is constructed from full-data fits; all coefficient initialization, dispersion estimation and grouped selection in each fold use training rows only. Folds are fitted sequentially; parallelize still controls response-channel parallelism within each fit. With start = "pilot", each fold constructs its own seed from its training rows, using training-weight-scaled pilot penalties. Full-data coefficients never initialize a CV training fit. A user-supplied numeric start is reused as supplied; the user is responsible for its provenance. Automatic channel preclustering likewise uses each fold's training responses and training column totals as size factors. Externally supplied control.precluster$clusterids or size_fac are reused as supplied; the user controls their provenance. Pilot support restriction fits the union selected by the pilots. Features outside that union remain zero throughout the path.

The fitter's penalty is on a total-likelihood scale. In each fold lambda is multiplied by the fraction of total prior weight in the training rows, so the returned lambda refers to the full-data penalty. Validation deviance is divided by validation prior weight, then averaged across folds with those weights. Its standard error is the square root of the weighted mean squared deviation of fold scores divided by nfolds - 1; for equal weights this is sd(fold_scores) / sqrt(nfolds). Error bars are plus/minus one SE. Quasi-families can be tuned by CV even when likelihood-based ICs are unavailable. The one-standard-error choice lambda_1se is the largest evaluated lambda with finite mean deviance no greater than the minimum mean deviance plus the SE at that minimum. Both panels mark this value in blue alongside the minimum-deviance choice and the GIC preset. fit_best continues to correspond to lambda_best. fit_1se contains the full-data fit at lambda_1se, retained from the same path without any additional fits.

Special families

NB, Cox and multinomial fits use independent initializations at each lambda. A supplied prior.mean is reused in the full-data, pilot and CV training fits, with normalization handled within each fit. For NB it centers regression coefficients, not theta; for Cox it centers log-relative-hazard coefficients, not the baseline hazard, and no intercept is included. Multinomial targets have K+1 category columns and are centered rowwise to CLR; a scalar or vector shared across categories therefore gives zero CLR targets. Nonzero targets cannot reactivate coefficients or rows excluded by start. Mandatory terms remain unpenalized, and special-family relaxed refits remove the ridge penalty on retained support. Adaptive weights and support counts continue to use coefficients relative to zero. The endpoint search tests the supplied prior for literal zero before CLR centering; use NULL or explicit zero targets for automatic refinement. Multinomial coefficient paths and seeds are CLR; all-zero starting rows exclude predictors, while individual category zeroes remain eligible. With $b$ denoting normalized working coefficients, grouped NB and multinomial adaptive ridge weights are $D_{jc}=\lambda/(4\|b_j\|_2^2)$; ungrouped weights are $D_{jc}=\lambda/(4|b_{jc}|^2)$. Cox uses the scalar rule. Grouped hard thresholding compares each row norm with thresh times the largest selectable row norm or the absolute floor 1e-20; it does not threshold individual entries within retained rows. Ungrouped multinomial thresholding compares individual CLR magnitudes within each category column, then recenters surviving entries within their allowed set. Excluded entries stay zero, and rows with fewer than two allowed entries are zero. Mandatory terms are excluded from cutoff maxima. Inverse-square weights and the same thresh are used in full-data candidates, pilot selection fits and CV training fits; each fit computes its own training normalization. Special-family cutoffs apply only during selection with positive lambda and are disabled during relaxed refits. thresh = 0 retains the absolute cutoff and does not smooth adaptive weights. With relaxed = TRUE, grouped multinomial refits jointly estimate all category coefficients in retained predictor rows; zero category entries within those rows may become nonzero. Ungrouped multinomial refits jointly estimate only the selected CLR entries and keep the others fixed at zero. Both preserve zero-sum rows: $s$ admitted category entries give $s-1$ free directions. Starting-row exclusions remain in force, mandatory terms remain eligible, and no other predictor rows can enter the refit. NB refits follow the same row-versus-entry support rule and continue updating estimated channel shapes. Cox refits selected log-relative-hazard effects. Refits apply inside each CV training fold, before validation scoring; full-data paths, fit_best and fit_1se contain the corresponding refitted models. Pilot initialization fits remain unrelaxed. Grouped multinomial IC tuning charges active rows and ungrouped tuning charges nonzero CLR entries. Parameter-count ICs account for the zero-sum constraint separately from these L0 objective charges. Grouped IC tuning charges rows; metrics$tuning_score is the selected row criterion plus mandatory/nuisance parameter constants, whereas the ordinary gic metric retains its parameter-count definition. CV trains theta, normalization and pilot starts using training observations only. Cox CV uses the difference between full-data and training partial log likelihood at training coefficients, divided by held-out event weight; this retains full risk sets. Multinomial CV uses deviance per weighted trial, and NB uses deviance per weighted response. The common lambda grid uses full data. An independent test set or nested CV can assess prediction error after tuning.

Examples

# Small versions of the L0glm.fit deconvolution and Gaussian ISIS examples.
# Five lambda values and three folds keep these examples short; increase
# nlambda and nfolds for a finer search. True coefficients are only used
# for plots and recovery metrics, never for fitting or choosing lambda.
old_par <- par(mgp=c(2, .7, 0))

# Single-channel Poisson deconvolution.
# Zoom around the GIC preset; use prop_range=1 for the full lambda range.
single <- simul(method=2, n=500, p=500, s=.1, seed=2)
single_gic <- L0glm.tune(single$X, single$y, family=poisson(identity),
  tune="gic", nlambda=5, prop_range=.25, coefs.true=single$beta_true,
  nonnegative=TRUE, relaxed=TRUE, parallelize=FALSE, plot=TRUE)

# Venetian-blind CV: train on two out of every three time points.
folds <- rep(1:3, length.out=nrow(single$X))
single_cv <- L0glm.tune(single$X, single$y, family=poisson(identity),
  tune="cv", nlambda=5, prop_range=.25, nfolds=3, foldid=folds,
  coefs.true=single$beta_true,
  nonnegative=TRUE, relaxed=TRUE, parallelize=FALSE, plot=TRUE)
c(GIC=single_gic$lambda_best, CV=single_cv$lambda_best,
  CV_1se=single_cv$lambda_1se, GIC_preset=single_cv$lambda_preset)
single_cv$metrics[c("lambda", "cv_deviance", "cv_se")]
L0glm.metrics(single_gic$fit_best, beta_true=single$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]
L0glm.metrics(single_cv$fit_best, beta_true=single$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]
# fit_1se is a full-data model, ready for prediction or coefficient extraction.
fitted(single_cv$fit_1se)[1:6, , drop=FALSE]

# Grouped multichannel Poisson deconvolution, keeping every channel.
multi <- simul(method=4, n=500, p=500, s=.05, seed=2)
multi_gic <- L0glm.tune(multi$X, multi$y, family=poisson(identity),
  tune="gic", nlambda=5, prop_range=.25, coefs.true=multi$beta_true,
  nonnegative=TRUE, grouped=TRUE, relaxed=TRUE, parallelize=FALSE, plot=TRUE)

# All channels of each time point share one fold assignment.
folds <- rep(1:3, length.out=nrow(multi$X))
multi_cv <- L0glm.tune(multi$X, multi$y, family=poisson(identity),
  tune="cv", nlambda=5, prop_range=.25, nfolds=3, foldid=folds,
  coefs.true=multi$beta_true,
  nonnegative=TRUE, grouped=TRUE, relaxed=TRUE, parallelize=FALSE, plot=TRUE)
c(GIC=multi_gic$lambda_best, CV=multi_cv$lambda_best,
  CV_1se=multi_cv$lambda_1se, GIC_preset=multi_cv$lambda_preset)
L0glm.metrics(multi_gic$fit_best, beta_true=multi$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC",
    "within_row_F1_score", "within_row_MCC")]
L0glm.metrics(multi_cv$fit_best, beta_true=multi$beta_true)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC",
    "within_row_F1_score", "within_row_MCC")]
dim(coef(multi_cv$fit_1se))            # all original predictor rows and channels

# High-dimensional Gaussian regression with 60 positive/negative effects.
# X alone needs about 4 GB; set p <- 20000L for a smaller example.
if (FALSE) { # \dontrun{
n <- 500L
p <- 1000000L
set.seed(123)
X_hd <- matrix(rnorm(n*p), n, p)
support <- sample.int(p, 60L)
beta_hd <- numeric(p)
beta_hd[support] <- runif(60L, -.8, .8)
mean_hd <- as.numeric(X_hd[, support, drop=FALSE] %*% beta_hd[support])
set.seed(1)
y_hd <- mean_hd + rnorm(n)

# screen="auto" applies p_eligible > 10*n independently at each lambda,
# with n equal to the training row count, including within each CV fold.
hd_gic <- L0glm.tune(X_hd, y_hd, family=gaussian(), tune="gic",
  nlambda=5, coefs.true=beta_hd, screen.nthreads=8, nthreads=1, plot=FALSE)
# Screening and preprocessing are recomputed within each training fold.
set.seed(42)
hd_cv <- L0glm.tune(X_hd, y_hd, family=gaussian(), tune="cv",
  nlambda=5, nfolds=3, coefs.true=beta_hd,
  screen.nthreads=8, nthreads=1, plot=FALSE)
c(GIC=hd_gic$lambda_best, CV=hd_cv$lambda_best,
  CV_1se=hd_cv$lambda_1se, GIC_preset=hd_cv$lambda_preset)
hd_cv$metrics[c("lambda", "cv_deviance", "cv_se")]
L0glm.metrics(hd_gic$fit_best, beta_true=beta_hd)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]
L0glm.metrics(hd_cv$fit_best, beta_true=beta_hd)[
  c("gic", "TP", "FP", "FN", "TN", "F1_score", "MCC")]
hd_cv$fit_best$screening$history
dim(coef(hd_cv$fit_1se))              # screened-out predictors remain zero rows
} # }

# CV plots show mean held-out deviance +/- one standard error. Vertical
# lines mark lambda_best, lambda_1se and the full-data GIC preset.
par(old_par)
# Sparse special families: GIC preset fits, GIC tuning and three-fold CV.
if (requireNamespace("mgcv", quietly=TRUE)) {
  set.seed(204)
  n <- 300L
  X_special <- matrix(rnorm(n*20), n, 20)
  folds_special <- sample(rep(1:3, length.out=n))
  family_par <- par(mgp=c(2, .7, 0))

  # Negative binomial: two nonzero slopes among 20 predictors.
  mu_nb <- exp(.4 + .8*X_special[, 1] - .6*X_special[, 2])
  y_nb <- rnbinom(n, mu=mu_nb, size=2)
  # nb() estimates theta; nb(theta=2) would keep it fixed at two.
  nb_fit <- L0glm.fit(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, lambda="gic", relaxed=TRUE, maxit=100, tol=1e-8)
  nb_fit$theta
  coef(nb_fit)
  nb_gic <- L0glm.tune(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, tune="gic", nlambda=5,
    relaxed=TRUE, maxit=100, tol=1e-8)
  nb_cv <- L0glm.tune(X_special, y_nb, family=mgcv::nb(),
    intercept=TRUE, tune="cv", nlambda=5, nfolds=3, foldid=folds_special,
    relaxed=TRUE, maxit=100, tol=1e-8)
  c(GIC=nb_gic$lambda_best, CV=nb_cv$lambda_best, CV_1se=nb_cv$lambda_1se)
  predict(nb_cv$fit_1se, X_special[1:5, , drop=FALSE], type="response")

  # Three-category multinomial: only two predictor rows affect probabilities.
  eta <- cbind(0, .9*X_special[, 1] - .7*X_special[, 2],
    -.8*X_special[, 1] + .6*X_special[, 2])
  probabilities <- exp(eta - apply(eta, 1, max))
  probabilities <- probabilities/rowSums(probabilities)
  y_multinomial <- factor(vapply(seq_len(n), function(i)
    sample.int(3, 1, prob=probabilities[i, ]), integer(1)), levels=1:3)
  multinomial_fit <- L0glm.fit(X_special, y_multinomial,
    family=mgcv::multinom(K=2), lambda="gic", grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  coef(multinomial_fit)               # original-scale, zero-sum CLR rows
  multinomial_gic <- L0glm.tune(X_special, y_multinomial,
    family=mgcv::multinom(K=2), tune="gic", nlambda=5, grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  multinomial_cv <- L0glm.tune(X_special, y_multinomial,
    family=mgcv::multinom(K=2), tune="cv", nlambda=5,
    nfolds=3, foldid=folds_special, grouped=TRUE,
    relaxed=TRUE, maxit=200, tol=1e-8)
  c(GIC=multinomial_gic$lambda_best, CV=multinomial_cv$lambda_best,
    CV_1se=multinomial_cv$lambda_1se)
  predict(multinomial_cv$fit_1se, X_special[1:5, , drop=FALSE], type="response")

  # Right-censored Cox model: two nonzero log-hazard effects, no intercept.
  if (requireNamespace("survival", quietly=TRUE)) {
    event_time <- rexp(n, rate=exp(.7*X_special[, 1] - .5*X_special[, 2]))
    censor_time <- rexp(n, rate=.4)
    y_cox <- survival::Surv(pmin(event_time, censor_time),
      event_time <= censor_time)
    cox_fit <- L0glm.fit(X_special, y_cox, family=mgcv::cox.ph(),
      lambda="gic", relaxed=TRUE, maxit=100, tol=1e-8)
    coef(cox_fit)
    cox_gic <- L0glm.tune(X_special, y_cox, family=mgcv::cox.ph(),
      tune="gic", nlambda=5, relaxed=TRUE, maxit=100, tol=1e-8)
    cox_cv <- L0glm.tune(X_special, y_cox, family=mgcv::cox.ph(),
      tune="cv", nlambda=5, nfolds=3, foldid=folds_special,
      relaxed=TRUE, maxit=100, tol=1e-8)
    c(GIC=cox_gic$lambda_best, CV=cox_cv$lambda_best, CV_1se=cox_cv$lambda_1se)
    predict(cox_cv$fit_1se, X_special[1:5, , drop=FALSE], type="risk")
  }
  par(family_par)
}