Sparse negative binomial regression

For overdispersed counts the NB2 model has

E(Y_i)=\mu_i,\qquad \operatorname{Var}(Y_i)=\mu_i+\mu_i^2/\theta.

The Poisson limit is approached as theta increases. mgcv::nb() estimates theta alongside the coefficients; mgcv::nb(theta=2) keeps it fixed. The default link is log.

Fit counts with a sparse signal

set.seed(21)
X <- matrix(rnorm(350*20), 350, 20)
mu <- exp(.4 + X[, 2] - .8*X[, 7])
y <- rnbinom(350, mu=mu, size=2)
fit <- L0glm.fit(X, y, family=mgcv::nb(), intercept=TRUE,
  lambda="gic", relaxed=TRUE, maxit=100, tol=1e-8)
fit
#> L0glm fit
#> Family: negative binomial   Link: log 
#> Observations: 350   Selectable predictors: 20   Channels: 1
#> Selected rows: 2   Nonzero selectable coefficients: 2   Mandatory rows: 1
#> Lambda: 5.3 
#> Convergence (last pass per stage/channel):
#>      stage channels converged not_converged unavailable max_iterations
#>  selection        1         1             0           0             29
#>    relaxed        1         1             0           0              7
fit$theta
#> [1] 2.9
predict(fit, X[1:5, ], type="response")
#>       [,1]
#> [1,] 3.581
#> [2,] 1.875
#> [3,] 0.320
#> [4,] 4.392
#> [5,] 0.557

For multichannel counts, supply an observation-by-channel matrix. grouped=TRUE selects predictor rows jointly; relaxed=TRUE refits eligible coefficients in the selected rows. Shapes may be estimated separately by channel.

Tune by GIC or CV

gic <- L0glm.tune(X, y, family=mgcv::nb(), intercept=TRUE,
  tune="gic", nlambda=10, relaxed=TRUE, maxit=100)
cv <- L0glm.tune(X, y, family=mgcv::nb(), intercept=TRUE,
  tune="cv", nfolds=5, nlambda=10, relaxed=TRUE, maxit=100)
predict(cv$fit_1se, X[1:5, ])

Each training fold estimates its own theta. Observation weights and offsets are supported. The negative-binomial reference details the likelihood and shape updates. The family convention follows mgcv::nb().