Gamma family with log link.
The Gamma family models strictly positive, continuous, right-skewed responses — for example, insurance claim sizes, waiting times, rainfall amounts, or other quantities that are bounded below by zero and become more variable as their mean grows. Although the canonical link for the Gamma distribution is the inverse, Whittaker uses the log link by default, since it guarantees positive fitted values and is generally easier to interpret (coefficients act multiplicatively on the response, as with Poisson). Use Gamma when the response is positive and continuous and its coefficient of variation is roughly constant across the range of fitted values; if instead the variance grows linearly with the mean, Poisson or Tweedie with 1 < p < 2 may fit better.
Notes
The (non-canonical, but default) link is the natural logarithm:
g(\mu) = \log(\mu)
The variance function grows with the square of the mean:
V(\mu) = \mu^2, \qquad \operatorname{Var}(Y) = \phi \, \mu^2,
so the coefficient of variation \sqrt{\operatorname{Var}(Y)} / \mu = \sqrt{\phi} is constant. The deviance is
D(y, \hat\mu) = 2 \sum_i \left[ -\log\!\left(\frac{y_i}{\hat\mu_i}\right)
+ \frac{y_i - \hat\mu_i}{\hat\mu_i} \right] .
Examples
Fit a GAM to positive, right-skewed data with a smooth trend:
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 + 0.4 * np.sin(x))
shape = 4.0
y = rng.gamma(shape, mu / shape)
data = {"x": x, "y": y}
model = wk.GAM("y ~ s(x)", family=wk.Gamma())
model.fit(data, method="REML")
print(model.summary())
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Gamma(link='log')
Inference: REML
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.4618 0.0360 12.819 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 4.79 5 59.894 1.278e-11
Total EDF: 5.79
Scale est: 0.259540
Deviance: 50.4055
Null dev: 67.0458
Dev. expl: 24.8%
GCV score: 0.267277
AIC: 446.54
BIC: 465.64
Attributes
|
Name
|
Description
|
|
scale_known
|
Whether the dispersion parameter is fixed. Always False for Gamma.
|
scale_known
Whether the dispersion parameter is fixed. Always False for Gamma.
The Gamma shape parameter (and hence the dispersion phi = 1/shape) is estimated from the data during fitting rather than fixed, unlike Poisson or Binomial where the scale is known in advance.
Methods
|
Name
|
Description
|
|
deviance()
|
Total Gamma deviance, the (weighted) sum of unit_deviance.
|
|
initialize()
|
Starting values for mu: y nudged away from zero.
|
|
link()
|
Apply the (default, non-canonical) 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()
|
Gamma log-likelihood parameterized by mean mu and shape \alpha = 1/\phi.
|
|
simulate()
|
Simulate Gamma-distributed response values with mean mu and dispersion scale.
|
|
unit_deviance()
|
Per-observation Gamma deviance contributions.
|
|
variance()
|
Gamma variance function: V(\mu) = \mu^2.
|
deviance()
Total Gamma deviance, the (weighted) sum of unit_deviance.
deviance(
y,
mu,
*,
weights=None,
)
Parameters
y: NDArray-
Observed response values (positive), 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, values are pushed away from zero to avoid log(0) on the first P-IRLS iteration.
Parameters
y: NDArray-
Observed response values (positive), shape
(n,).
Returns
NDArray-
Starting values for
mu, shape (n,).
link()
Apply the (default, non-canonical) 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()
Gamma log-likelihood parameterized by mean mu and shape \alpha = 1/\phi.
log_likelihood(
y,
mu,
scale,
*,
weights=None,
)
Evaluates the Gamma log-density at each observation using the shape parameter \alpha = 1/\text{scale} and rate \alpha/\mu_i, then sums (optionally weighted) over observations.
Parameters
y: NDArray-
Observed response values (positive), shape (n,).
mu: NDArray-
Fitted conditional mean values, shape (n,).
scale: float-
Dispersion parameter \phi = 1/\alpha.
weights: NDArray | None = None-
Optional prior weights, shape
(n,).
Returns
float-
The total log-likelihood.
simulate()
Simulate Gamma-distributed response values with mean mu and dispersion scale.
simulate(
mu,
scale,
rng,
)
Parameters
mu: NDArray-
Mean (fitted values), shape (n,).
scale: float-
Estimated dispersion parameter \phi = 1/\alpha.
rng: np.random.Generator-
A
numpy.random.Generator instance.
Returns
NDArray-
Simulated response values, shape
(n,).
unit_deviance()
Per-observation Gamma deviance contributions.
Computes d_i = 2 \left[ -\log(y_i/\hat\mu_i) + (y_i - \hat\mu_i)/\hat\mu_i \right].
Parameters
y: NDArray-
Observed response values (positive), shape (n,).
mu: NDArray-
Fitted conditional mean values, shape
(n,).
Returns
NDArray-
Per-observation deviance contributions, shape
(n,).
variance()
Gamma variance function: V(\mu) = \mu^2.
The variance grows with the square of the mean, so the coefficient of variation \sqrt{\operatorname{Var}(Y)}/\mu = \sqrt\phi is constant across the range of mu.
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Variance values \mu^2, shape
(n,).