Fit a negative-binomial GLM with simultaneous (beta, theta) MLE

`fastglm_nb()` fits a negative-binomial regression model, jointly maximising the likelihood over the regression coefficients `beta` and the NB2 dispersion `theta`. It is the fastglm analogue of [MASS::glm.nb()], built on top of the native NB family kernel introduced in 0.0.6 so that all numerical loops – IRLS, the inner theta MLE Brent root-find, and the outer (beta, theta) alternation – run entirely in C++.

fastglm_nb(
  x,
  y,
  weights = NULL,
  offset = NULL,
  start = NULL,
  init.theta = NULL,
  link = c("log", "sqrt", "identity"),
  method = 2L,
  tol = 1e-08,
  maxit = 100L,
  outer.maxit = 25L,
  outer.tol = 1e-07,
  theta.tol = 1e-08,
  theta.maxit = 100L
)

Arguments

x: design matrix (numeric matrix, dgCMatrix not yet supported here).
y: non-negative integer response vector.
weights: optional prior weights vector of length `length(y)`.
offset: optional offset vector of length `length(y)`.
start: optional starting values for `beta`.
init.theta: optional starting value for `theta`. If `NULL`, uses the method-of-moments estimator from a Poisson pilot fit.
link: character, one of `"log"` (default), `"sqrt"`, `"identity"`.
method: integer; `0..5`, see [fastglm()].
tol: convergence tolerance for the IRLS inner loop.
maxit: maximum number of inner-loop IRLS iterations.
outer.maxit: maximum number of `(beta, theta)` outer iterations.
outer.tol: convergence tolerance for the outer loop on the sup-norm of the beta update plus the relative change in `theta`.
theta.tol: Brent tolerance for the inner theta MLE.
theta.maxit: max iterations for the inner theta MLE.

Value

A list of class `c("fastglm_nb", "fastglm")` with the usual fastglm components plus `theta`, `SE.theta`, `iter.theta`, and `twologlik` (twice the maximised NB log-likelihood).

Examples

set.seed(1)
n <- 500
x <- cbind(1, matrix(rnorm(n * 2), n, 2))
eta <- x %*% c(0.3, 0.5, -0.2)
mu  <- exp(eta)
if (requireNamespace("MASS", quietly = TRUE)) {
  y   <- MASS::rnegbin(n, mu = mu, theta = 2)
  fit <- fastglm_nb(x, y)
  c(theta = fit$theta, fit$coefficients)
}
#>      theta         x1         x2         x3 
#>  1.9625479  0.3059098  0.5879277 -0.1601587