Negative Binomial family with log link (NB2 parameterization).
NegativeBinomial(theta=1.0)
The Negative Binomial family models count data that is overdispersed relative to the Poisson distribution, i.e. where the observed variance exceeds the mean. This commonly arises when counts are driven by unobserved heterogeneity across observations (e.g. some individuals or locations are systematically more prone to events than others). Whittaker uses the NB2 parameterization, where the variance function is V(mu) = mu + mu^2/theta and theta controls the degree of overdispersion: as theta -> infinity the distribution converges to Poisson, while smaller theta implies heavier overdispersion. The canonical log link is used, giving the same multiplicative interpretation of coefficients as Poisson.
Parameters
theta: float = 1.0-
Overdispersion (size) parameter, must be positive. Smaller values of
theta imply greater overdispersion; larger values make the distribution approach Poisson. The value supplied at construction is used as the starting point and is refined during fitting via an outer iteration around P-IRLS unless explicitly held fixed by the caller.
Notes
The canonical link is the natural logarithm:
g(\mu) = \log(\mu)
The variance function is
V(\mu) = \mu + \frac{\mu^2}{\theta},
so the variance always exceeds the mean by the extra term \mu^2/\theta. The deviance is
D(y, \hat\mu) = 2 \sum_i \left[ y_i \log\!\left(\frac{y_i}{\hat\mu_i}\right)
- (y_i + \theta) \log\!\left(\frac{y_i + \theta}{\hat\mu_i + \theta}\right) \right] .
Examples
Fit a GAM to overdispersed count data:
import numpy as np
import whittaker as wk
rng = np.random.default_rng(0)
n = 200
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(0.5 * np.sin(x))
theta = 3.0
y = rng.negative_binomial(theta, theta / (theta + mu))
data = {"x": x, "y": y}
model = wk.GAM("y ~ s(x)", family=wk.NegativeBinomial(theta=theta))
model.fit(data, method="REML")
print(model.summary())
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: NegativeBinomial(theta=4.721, link='log')
Inference: REML
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) -0.0054 0.0792 -0.069 0.9451
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 3.54 4 15.493 0.003781
Total EDF: 4.54
Scale est: 1.000000
Deviance: 210.5567
Null dev: 229.3245
Dev. expl: 8.2%
GCV score: 1.102291
AIC: 543.06
BIC: 558.05
Attributes
|
Name
|
Description
|
|
scale_known
|
Whether the dispersion parameter is fixed. Always True for NegativeBinomial.
|
|
theta
|
Overdispersion (size) parameter theta of the NB2 distribution.
|
scale_known
Whether the dispersion parameter is fixed. Always True for NegativeBinomial.
The NB2 dispersion is governed entirely by theta rather than a separate scale parameter estimated during P-IRLS, so the scale is treated as fixed at 1. theta itself is refined by an outer iteration around P-IRLS rather than the usual scale estimation.
theta
Overdispersion (size) parameter theta of the NB2 distribution.
Smaller values imply greater overdispersion relative to Poisson; larger values make the distribution approach Poisson as theta -> infinity. This value is used by variance, unit_deviance, log_likelihood, and simulate.
Methods
|
Name
|
Description
|
|
deviance()
|
Total NB2 deviance, the (weighted) sum of unit_deviance.
|
|
initialize()
|
Starting values for mu: y nudged away from zero.
|
|
link()
|
Apply the log link: \eta = \log(\mu).
|
|
link_derivative()
|
Derivative of the log link: g'(\mu) = 1/\mu.
|
|
link_inverse()
|
Apply the inverse log link: \mu = e^{\eta}.
|
|
log_likelihood()
|
NB2 log-likelihood evaluated using the current theta.
|
|
simulate()
|
Simulate NB2-distributed response values with mean mu and the current theta.
|
|
unit_deviance()
|
Per-observation NB2 deviance contributions.
|
|
variance()
|
NB2 variance function: V(\mu) = \mu + \mu^2/\theta.
|
deviance()
Total NB2 deviance, the (weighted) sum of unit_deviance.
deviance(
y,
mu,
*,
weights=None,
)
Parameters
y: NDArray-
Observed response values (counts), shape (n,).
mu: NDArray-
Fitted conditional mean values, shape (n,).
weights: NDArray | None = None-
Optional prior weights, shape
(n,).
Returns
float-
The total (weighted) deviance.
initialize()
Starting values for mu: y nudged away from zero.
Since the log link requires strictly positive mu, small or zero counts are pushed away from zero to avoid log(0) on the first P-IRLS iteration.
Parameters
y: NDArray-
Observed response values (counts), shape
(n,).
Returns
NDArray-
Starting values for
mu, shape (n,).
link()
Apply the log link: \eta = \log(\mu).
Parameters
mu: NDArray-
Conditional mean values, shape
(n,). Must be positive.
Returns
NDArray-
Linear predictor values \eta = \log(\mu), shape
(n,).
link_derivative()
Derivative of the log link: g'(\mu) = 1/\mu.
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Derivative values 1/\mu, shape
(n,).
link_inverse()
Apply the inverse log link: \mu = e^{\eta}.
The linear predictor is clipped to [-30, 30] before exponentiating to guard against overflow while fitting.
Parameters
eta: NDArray-
Linear predictor values, shape
(n,).
Returns
NDArray-
Conditional mean values \mu = e^{\eta}, shape
(n,).
log_likelihood()
NB2 log-likelihood evaluated using the current theta.
log_likelihood(
y,
mu,
scale,
*,
weights=None,
)
The scale argument is accepted for interface compatibility but ignored; the overdispersion is instead governed by theta (see scale_known).
Parameters
y: NDArray-
Observed response values (counts), shape (n,).
mu: NDArray-
Fitted conditional mean values, shape (n,).
scale: float-
Ignored.
weights: NDArray | None = None-
Optional prior weights, shape
(n,).
Returns
float-
The total log-likelihood.
simulate()
Simulate NB2-distributed response values with mean mu and the current theta.
simulate(
mu,
scale,
rng,
)
Parameters
mu: NDArray-
Mean (fitted values), shape (n,).
scale: float-
Ignored; overdispersion is governed by theta.
rng: np.random.Generator-
A
numpy.random.Generator instance.
Returns
NDArray-
Simulated response values, shape
(n,).
unit_deviance()
Per-observation NB2 deviance contributions.
Computes d_i = 2 \left[ y_i \log(y_i/\hat\mu_i)
- (y_i + \theta) \log\!\left(\frac{y_i+\theta}{\hat\mu_i+\theta}\right) \right], with the usual convention that the y_i \log(\cdot) term vanishes when y_i = 0.
Parameters
y: NDArray-
Observed response values (counts), shape (n,).
mu: NDArray-
Fitted conditional mean values, shape
(n,).
Returns
NDArray-
Per-observation deviance contributions, shape
(n,).
variance()
NB2 variance function: V(\mu) = \mu + \mu^2/\theta.
The variance always exceeds the mean by the extra term \mu^2/\theta, which vanishes as theta -> infinity (recovering the Poisson variance V(mu) = mu).
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Variance values \mu + \mu^2/\theta, shape
(n,).