Ordinal Regression¶
Proportional odds (cumulative link) models matching R’s MASS::polr(). Supports logistic, probit, and complementary log-log links.
Ordinal regression (proportional odds / cumulative link models).
Provides the proportional odds model matching R’s MASS::polr(), with logistic, probit, and complementary log-log link functions.
- Public API:
polr(y, X, …) - Fit a proportional odds model OrdinalSolution - Result wrapper with R-style summary()
- pystatistics.ordinal.polr(y, X, *, link='logistic', names=None, category_names=None, tol=1e-08, max_iter=200, l2=0.0, backend=None, conf_level=0.95)[source]¶
Fit a proportional odds (cumulative link) model.
Fits the model P(Y <= j | x) = g^{-1}(alpha_j - x’beta), matching R’s MASS::polr(). The proportional odds assumption means all categories share the same slope coefficients beta; only the threshold (intercept) parameters alpha_j differ.
- Parameters:
y (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Ordered categorical response as integer codes 0, 1, …, K-1. Must contain all levels from 0 to max(y) with no gaps.
X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Design matrix of shape (n, p). Must NOT include an intercept column (thresholds serve as category-specific intercepts).
link (str) – Link function. One of ‘logistic’ (default, proportional odds), ‘probit’, or ‘cloglog’.
names (list[str] | None) – Labels for the p predictor columns. If None, defaults to ‘x1’, ‘x2’, etc.
category_names (list[str] | None) – Labels for the K ordered categories. If None, defaults to ‘0’, ‘1’, etc.
tol (float) – Convergence tolerance (gradient norm) for L-BFGS-B.
max_iter (int) – Maximum number of optimizer iterations.
l2 (float) – L2 (ridge) penalty coefficient on the slope coefficients beta (the threshold/intercept parameters are never penalized). The penalty
0.5 * l2 * ||beta||^2is added to the negative log-likelihood. The default 0.0 fits the exact, R-validated maximum-likelihood model. A small positive value makes the objective strongly convex, which keeps the fit finite and well-conditioned under (quasi-)complete separation — where the unpenalized slopes would diverge, L-BFGS-B would run tomax_iter, and the fit would be rejected. Reported coefficients, standard errors, and vcov are then those of the penalized fit. CPU backend only; the GPU backend applies its own ridge stabilization.backend (str | None) – Compute backend = (device, precision).
'cpu'(default): the R-validated reference path.'gpu': float32 (raises if no GPU).'gpu_fp64': float64, CUDA only (raises on MPS).'auto': GPU-float32 if CUDA is present, else CPU.conf_level (float)
- Returns:
OrdinalSolution wrapping the fitted model parameters, with properties for coefficients, thresholds, standard errors, z-values, p-values, and a summary() method matching R’s print.summary.polr output.
- Raises:
ValidationError – If inputs fail validation (wrong shapes, non-integer y, missing levels, etc.).
ConvergenceError – If the optimizer fails to converge.
- Return type:
Examples
>>> import numpy as np >>> from pystatistics.ordinal import polr >>> rng = np.random.default_rng(42) >>> n = 500 >>> X = rng.standard_normal((n, 2)) >>> eta = 1.0 * X[:, 0] - 0.5 * X[:, 1] >>> # Generate ordinal y with 3 levels >>> cum_p1 = 1 / (1 + np.exp(-(-1 - eta))) >>> cum_p2 = 1 / (1 + np.exp(-(1 - eta))) >>> u = rng.uniform(size=n) >>> y = np.where(u < cum_p1, 0, np.where(u < cum_p2, 1, 2)) >>> sol = polr(y, X, link='logistic', names=['x1', 'x2']) >>> print(sol.summary())
- class pystatistics.ordinal.OrdinalSolution(_result, _names, _level_names, _conf_level=0.95)[source]¶
Bases:
SolutionReprMixinUser-facing ordinal regression results.
Wraps Result[OrdinalParams] and provides convenient access to fitted model quantities. The summary() method produces output matching R’s summary(polr(…)).
- Parameters:
- property coef: dict[str, float]¶
Slope coefficients as a name-to-value dictionary.
- Returns:
Dict mapping predictor names to their fitted beta values.
- property coefficients: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Slope coefficient vector beta.
- Returns:
Array of shape (p,) with fitted slope parameters.
- property thresholds: dict[str, float]¶
Threshold parameters with level-pair labels.
- Returns:
Dict mapping threshold labels (e.g. ‘low|medium’) to their fitted alpha values.
- property threshold_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Raw threshold parameter vector alpha.
- Returns:
Array of shape (K-1,) with ordered threshold values.
- property vcov: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Full variance-covariance matrix.
Ordered as [thresholds, coefficients], shape (K-1+p, K-1+p). Reported in natural threshold coordinates (matching MASS::polr).
- Returns:
Variance-covariance matrix.
- property standard_errors: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Standard errors for slope coefficients beta.
- Returns:
Array of shape (p,) with SE(beta).
- property threshold_standard_errors: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Standard errors for threshold parameters.
These are in natural threshold coordinates (matching MASS::polr).
- Returns:
Array of shape (K-1,) with SE(alpha).
- property z_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
beta / SE(beta).
- Returns:
Array of shape (p,).
- Type:
Wald z-statistics for slope coefficients
- property p_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Two-sided p-values for slope coefficients from z-statistics.
- Returns:
Array of shape (p,).
- property conf_int: ndarray[tuple[Any, ...], dtype[floating[Any]]]¶
Wald confidence intervals for the slope coefficients, shape (p, 2).
beta ± z * SE(beta)with the normal quantile forconf_level(proportional-odds inference is asymptotic-normal).exp(conf_int)gives odds-ratio intervals. Thresholds are not included (seethresholds/threshold_standard_errors).