Ordered categorical (proportional odds / cumulative logit) family.
OrderedCategorical(n_categories)
OrderedCategorical models an ordinal response — one with a small number of categories that have a natural order but no meaningful numeric spacing, such as a Likert scale (“disagree” / “neutral” / “agree”) or a severity grade — using the proportional-odds cumulative logit model. A single set of K - 1 ordered cutpoints (thresholds) alpha_1 < alpha_2 < ... < alpha_{K-1} is estimated jointly with the smooth/linear predictor eta, and each cutpoint defines a binary split between “category k or below” versus “above category k”. Because the same eta (and hence the same covariate effects) is shared across all thresholds, covariate effects are assumed to shift the log-odds of being in a higher category by the same amount at every threshold — the proportional-odds assumption. Use this family for ordinal outcomes with three or more ordered levels; for a two-level (binary) response, use Binomial instead, and for unordered categorical responses, use Multinomial.
Parameters
n_categories: int-
Number of ordered response categories
K (must be >= 2). Responses passed to GAM.fit() should be integer-coded 1, 2, ..., K.
Notes
The response, without an intercept in the design matrix (since the cutpoints absorb it), is modeled through cumulative probabilities:
P(Y \le k \mid \eta) = \operatorname{expit}(\alpha_k - \eta), \qquad k = 1, \dots, K-1,
which is equivalent to a logit link on each cumulative probability, g(P(Y \le k)) = \alpha_k - \eta. Category probabilities follow by differencing:
P(Y = 1) = \operatorname{expit}(\alpha_1 - \eta), \qquad
P(Y = K) = 1 - \operatorname{expit}(\alpha_{K-1} - \eta),
and for interior categories 1 < k < K,
P(Y = k) = \operatorname{expit}(\alpha_k - \eta) - \operatorname{expit}(\alpha_{k-1} - \eta).
Because this loss does not fit the standard GLM deviance framework, link and link_inverse are the identity on eta, and fitting instead uses a custom irls_update together with an inner maximum-likelihood step (_update_cutpoints) that re-estimates the cutpoints alpha at each P-IRLS iteration. The deviance reported is -2 times the multinomial log-likelihood of the observed categories under the fitted probabilities.
Examples
Fit a GAM to a four-level ordinal response with a smooth covariate effect:
import numpy as np
import whittaker as wk
from scipy.special import expit
rng = np.random.default_rng(0)
n = 300
x = np.linspace(-3, 3, n)
eta = np.sin(x)
cutpoints = np.array([-1.5, 0.0, 1.5])
y = np.empty(n)
for i in range(n):
probs = np.diff(
np.concatenate([[0.0], expit(cutpoints - eta[i]), [1.0]])
)
y[i] = rng.choice([1, 2, 3, 4], p=probs)
data = {"x": x, "y": y}
model = wk.GAM("y ~ s(x)", family=wk.OrderedCategorical(n_categories=4))
model.fit(data, method="REML")
print(model.summary())
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: OrderedCategorical(K=4)
Inference: REML
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.0143 0.0887 0.161 0.8722
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 5.02 6 51.396 2.466e-09
Total EDF: 6.02
Scale est: 1.000000
Deviance: 785.2424
Null dev: 1365.4479
Dev. expl: 42.5%
GCV score: 2.725803
AIC: 797.29
BIC: 819.59
Attributes
|
Name
|
Description
|
|
cutpoints
|
Fitted cutpoints (thresholds) alpha_1 < alpha_2 < ... < alpha_{K-1}.
|
|
n_categories
|
Number of ordered response categories K.
|
|
scale_known
|
Whether the dispersion (scale) parameter is fixed rather than estimated.
|
cutpoints
Fitted cutpoints (thresholds) alpha_1 < alpha_2 < ... < alpha_{K-1}.
cutpoints: NDArray | None
Each cutpoint alpha_k defines the boundary of the cumulative logit P(Y \le k \mid \eta) = \operatorname{expit}(\alpha_k - \eta) separating category k or below from categories above k. The cutpoints are estimated jointly with eta and are updated at every P-IRLS iteration by _update_cutpoints.
n_categories
Number of ordered response categories K.
This is the value supplied to the constructor and equals 1 + the number of cutpoints (thresholds) fitted by the model.
scale_known
Whether the dispersion (scale) parameter is fixed rather than estimated.
Overrides the base Family.scale_known (which returns False) because the multinomial log-likelihood underlying the proportional-odds model has no free dispersion parameter — the scale is always 1 and is never estimated during fitting.
Methods
|
Name
|
Description
|
|
deviance()
|
Total deviance, -2 times the multinomial log-likelihood of the observed categories.
|
|
initialize()
|
Initialize cutpoints and the starting linear predictor before P-IRLS.
|
|
irls_update()
|
Working response and weights for the cumulative-logit log-likelihood.
|
|
link()
|
Identity link, g(mu) = mu.
|
|
link_derivative()
|
Derivative of the identity link, g'(mu) = 1.
|
|
link_inverse()
|
Identity inverse link, g^{-1}(eta) = eta.
|
|
log_likelihood()
|
Log-likelihood of the observed categories, -0.5 * deviance(y, mu).
|
|
simulate()
|
Simulate ordinal category labels from the fitted cumulative-logit model.
|
|
unit_deviance()
|
Per-observation deviance contributions, overriding the base (y - mu)^2 default.
|
|
variance()
|
Placeholder variance function, V(mu) = 1.
|
deviance()
Total deviance, -2 times the multinomial log-likelihood of the observed categories.
deviance(
y,
mu,
*,
weights=None,
)
Implements the family-specific version of the abstract Family.deviance. Because the proportional-odds model has a saturated log-likelihood of 0 (each observation can be assigned probability 1 to its own category), the deviance reduces to
D(y, \eta) = -2 \sum_i \log P(Y_i = y_i \mid \eta_i),
with the cumulative category probabilities P(Y_i = k \mid \eta_i) computed by _category_probs from the current cutpoints. The weights argument is accepted for interface compatibility but is not applied.
Parameters
y: NDArray-
Observed category labels, coded 1, ..., K, shape (n,).
mu: NDArray-
Linear predictor values eta (see link), shape (n,).
weights: NDArray | None = None-
Unused.
Returns
float-
The total deviance, or
len(y) as a placeholder if the model has not yet been fitted (cutpoints not initialized).
initialize()
Initialize cutpoints and the starting linear predictor before P-IRLS.
Overrides the base Family.initialize (which returns y unchanged). Sets the initial cutpoints via _init_cutpoints, which derives each threshold alpha_k from the empirical logit of P(Y <= k + 1), and starts the linear predictor eta at zero for every observation (i.e. an intercept-only model at the first iteration).
Parameters
y: NDArray-
Observed category labels, coded
1, ..., K, shape (n,).
Returns
NDArray-
Starting linear predictor values, an array of zeros with shape
(n,).
irls_update()
Working response and weights for the cumulative-logit log-likelihood.
irls_update(
y,
mu,
eta,
)
Implements the family-specific override of the base Family.irls_update (whose default returns None and falls back to the standard GLM formula). On the first call the cutpoints are initialized via _init_cutpoints; on every call the cutpoints are then re-estimated by maximum likelihood via _update_cutpoints given the current eta, implementing the two-block (cutpoints, then smooth) coordinate ascent that fits the proportional-odds model. Category probabilities are computed by _category_probs, and for each observation i in its observed category k = y_int[i] - 1 the first and second derivatives of the log-likelihood
\ell_i(\eta_i) = \log\bigl[P(Y_i = k+1 \mid \eta_i)\bigr]
with respect to eta_i are computed in closed form (differing at the boundary categories k = 0 and k = K - 1 versus interior categories). The working weight is the (clipped, non-negative) negative second derivative, and the working response is the Newton step z_i = \eta_i + (d\ell_i/d\eta_i) / W_i.
Parameters
y: NDArray-
Observed category labels, coded 1, ..., K, shape (n,).
mu: NDArray-
Unused; OrderedCategorical uses eta directly (see link).
eta: NDArray-
Current linear predictor values, shape
(n,).
Returns
tuple of NDArray, NDArray-
(z, W), the working response and working weights, each shape (n,).
link()
Identity link, g(mu) = mu.
OrderedCategorical does not use a conventional mean-to-linear-predictor link: the object passed around as mu throughout P-IRLS is already the linear predictor eta, and the actual cumulative-logit transformation is applied internally by _category_probs and irls_update. This method and link_inverse are therefore the identity, present only to satisfy the Family interface.
Parameters
mu: NDArray-
Values to pass through unchanged, shape
(n,).
link_derivative()
Derivative of the identity link, g'(mu) = 1.
Since link is the identity, its derivative is constant. This method is unused in practice because OrderedCategorical overrides irls_update to compute the working response and weights directly from the cumulative-logit log-likelihood, bypassing the standard GLM link/variance formula.
Parameters
mu: NDArray-
Values used only to determine the output shape, shape
(n,).
Returns
NDArray-
Array of ones, shape
(n,).
link_inverse()
Identity inverse link, g^{-1}(eta) = eta.
Mirrors link: since OrderedCategorical treats mu and eta as the same quantity and derives category probabilities directly from eta via _category_probs, no transformation is needed here.
Parameters
eta: NDArray-
Values to pass through unchanged, shape
(n,).
log_likelihood()
Log-likelihood of the observed categories, -0.5 * deviance(y, mu).
log_likelihood(
y,
mu,
scale,
*,
weights=None,
)
Implements the family-specific version of the abstract Family.log_likelihood. Since deviance is already defined as -2 times the multinomial log-likelihood, this simply inverts that relationship; the scale parameter is ignored because OrderedCategorical has no free dispersion parameter (see scale_known).
Parameters
y: NDArray-
Observed category labels, coded 1, ..., K, shape (n,).
mu: NDArray-
Linear predictor values eta (see link), shape (n,).
scale: float-
Unused; the scale is fixed at 1 for this family.
weights: NDArray | None = None-
Optional prior weights, forwarded to
deviance.
Returns
float-
The total multinomial log-likelihood.
simulate()
Simulate ordinal category labels from the fitted cumulative-logit model.
simulate(
mu,
scale,
rng,
)
Implements the family-specific version of the abstract Family.simulate. For each observation, category probabilities are computed from the current cutpoints and eta via _category_probs, and a category in 1, ..., K is drawn from the resulting categorical distribution using rng.choice.
Parameters
mu: NDArray-
Linear predictor values eta (see link), shape (n,).
scale: float-
Unused; this family has no free dispersion parameter (see scale_known).
rng: np.random.Generator-
A
numpy.random.Generator instance used to draw the simulated categories.
Returns
NDArray-
Simulated category labels, coded
1, ..., K, shape (n,).
Notes
Raises RuntimeError if called before the model has been fitted (i.e. before initialize has set cutpoints).
unit_deviance()
Per-observation deviance contributions, overriding the base (y - mu)^2 default.
Implements the family-specific per-observation deviance so that the total returned by deviance decomposes as sum(unit_deviance(y, mu)). For observation i in category y_i, the contribution is
d_i = -2 \log P(Y_i = y_i \mid \eta_i),
using the cumulative category probabilities from _category_probs.
Parameters
y: NDArray-
Observed category labels, coded 1, ..., K, shape (n,).
mu: NDArray-
Linear predictor values
eta (see link), shape (n,).
Returns
NDArray-
Per-observation deviance contributions d_i, shape
(n,), or an array of ones if the model has not yet been fitted (cutpoints not initialized).
variance()
Placeholder variance function, V(mu) = 1.
OrderedCategorical does not fit the standard GLM variance-function framework: its working weights are instead derived from the second derivative of the multinomial log-likelihood with respect to eta inside irls_update. This method returns a constant and is not used by the fitting loop.
Parameters
mu: NDArray-
Values used only to determine the output shape, shape
(n,).
Returns
NDArray-
Array of ones, shape
(n,).