# OrderedCategorical


Ordered categorical (proportional odds / cumulative logit) family.


Usage

``` python
OrderedCategorical(n_categories)
```


[OrderedCategorical](OrderedCategorical.md#whittaker.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](Binomial.md#whittaker.Binomial) instead, and for unordered categorical responses, use [Multinomial](Multinomial.md#whittaker.Multinomial).


## Parameters


`n_categories: int`  
Number of ordered response categories `K` (must be `>= 2`). Responses passed to [GAM.fit()](GAM.md#whittaker.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:


``` python
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](#cutpoints) | Fitted cutpoints (thresholds) `alpha_1 < alpha_2 < ... < alpha_{K-1}`. |
| [n_categories](#n_categories) | Number of ordered response categories `K`. |
| [scale_known](#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`.


`n_categories: int`


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.


`scale_known: bool`


Overrides the base [Family.scale_known](Family.md#whittaker.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()](#deviance) | Total deviance, -2 times the multinomial log-likelihood of the observed categories. |
| [initialize()](#initialize) | Initialize cutpoints and the starting linear predictor before P-IRLS. |
| [irls_update()](#irls_update) | Working response and weights for the cumulative-logit log-likelihood. |
| [link()](#link) | Identity link, `g(mu) = mu`. |
| [link_derivative()](#link_derivative) | Derivative of the identity link, `g'(mu) = 1`. |
| [link_inverse()](#link_inverse) | Identity inverse link, `g^{-1}(eta) = eta`. |
| [log_likelihood()](#log_likelihood) | Log-likelihood of the observed categories, `-0.5 * deviance(y, mu)`. |
| [simulate()](#simulate) | Simulate ordinal category labels from the fitted cumulative-logit model. |
| [unit_deviance()](#unit_deviance) | Per-observation deviance contributions, overriding the base `(y - mu)^2` default. |
| [variance()](#variance) | Placeholder variance function, `V(mu) = 1`. |

------------------------------------------------------------------------


### deviance()


Total deviance, -2 times the multinomial log-likelihood of the observed categories.


Usage

``` python
deviance(
    y,
    mu,
    *,
    weights=None,
)
```


Implements the family-specific version of the abstract [Family.deviance](Family.md#whittaker.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.


Usage

``` python
initialize(y)
```


Overrides the base [Family.initialize](Family.md#whittaker.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.


Usage

``` python
irls_update(
    y,
    mu,
    eta,
)
```


Implements the family-specific override of the base [Family.irls_update](Family.md#whittaker.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](OrderedCategorical.md#whittaker.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`.


Usage

``` python
link(mu)
```


[OrderedCategorical](OrderedCategorical.md#whittaker.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](Family.md#whittaker.Family) interface.


#### Parameters


`mu: NDArray`  
Values to pass through unchanged, shape `(n,)`.


#### Returns


`NDArray`  
`mu`, unchanged.


------------------------------------------------------------------------


### link_derivative()


Derivative of the identity link, `g'(mu) = 1`.


Usage

``` python
link_derivative(mu)
```


Since `link` is the identity, its derivative is constant. This method is unused in practice because [OrderedCategorical](OrderedCategorical.md#whittaker.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`.


Usage

``` python
link_inverse(eta)
```


Mirrors `link`: since [OrderedCategorical](OrderedCategorical.md#whittaker.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,)`.


#### Returns


`NDArray`  
`eta`, unchanged.


------------------------------------------------------------------------


### log_likelihood()


Log-likelihood of the observed categories, `-0.5 * deviance(y, mu)`.


Usage

``` python
log_likelihood(
    y,
    mu,
    scale,
    *,
    weights=None,
)
```


Implements the family-specific version of the abstract [Family.log_likelihood](Family.md#whittaker.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](OrderedCategorical.md#whittaker.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.


Usage

``` python
simulate(
    mu,
    scale,
    rng,
)
```


Implements the family-specific version of the abstract [Family.simulate](Family.md#whittaker.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](OrderedCategorical.md#whittaker.OrderedCategorical.cutpoints)).

------------------------------------------------------------------------


### unit_deviance()


Per-observation deviance contributions, overriding the base `(y - mu)^2` default.


Usage

``` python
unit_deviance(
    y,
    mu,
)
```


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`.


Usage

``` python
variance(mu)
```


[OrderedCategorical](OrderedCategorical.md#whittaker.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,)`.