Ordered logit and probit

Download the script: ordered.py

\[ P(y_i \le j \mid x_i) = F(\kappa_j - x_i'\beta), \qquad \kappa_1 < \dots < \kappa_{J-1}, \]

with \(F\) the logistic distribution function (ologit_stream) or the standard normal’s (oprobit_stream): a positive coefficient moves probability toward the higher levels.

import polars as pl
from logit_stream import ologit_stream, oprobit_stream
from logit_stream.simulate import EDUCATION, simulate_people

people = simulate_people(200_000, seed=1)
X = "age10 + female + log_income + owner + urban"
fit = ologit_stream(f"education ~ {X}", people, vcov="hetero")
fit.tidy()

┌─────────────┬───────────┬────────────┬────────────┬─────────────┐
│ Coefficient ┆ Estimate  ┆ Std. Error ┆ t value    ┆ Pr(>|t|)    │
╞═════════════╪═══════════╪════════════╪════════════╪═════════════╡
│ age10       ┆ 0.045762  ┆ 0.002861   ┆ 15.996862  ┆ 1.3438e-57  │
│ female      ┆ 0.242993  ┆ 0.008052   ┆ 30.178821  ┆ 4.4924e-200 │
│ log_income  ┆ 0.811634  ┆ 0.005046   ┆ 160.833563 ┆ 0.0         │
│ owner       ┆ 0.281224  ┆ 0.011119   ┆ 25.291536  ┆ 3.9581e-141 │
│ urban       ┆ 0.265801  ┆ 0.008153   ┆ 32.603152  ┆ 3.7004e-233 │
│ cut1        ┆ 6.756594  ┆ 0.052894   ┆ 127.73735  ┆ 0.0         │
│ cut2        ┆ 8.192264  ┆ 0.053519   ┆ 153.071066 ┆ 0.0         │
│ cut3        ┆ 9.237465  ┆ 0.054303   ┆ 170.108772 ┆ 0.0         │
│ cut4        ┆ 10.483533 ┆ 0.055288   ┆ 189.617347 ┆ 0.0         │
└─────────────┴───────────┴────────────┴────────────┴─────────────┘

fit.cutpoints is \(\kappa_1, \dots, \kappa_4\) as an array. The slopes are close to the ones simulate_people draws from (0.05, 0.25, 0.85, 0.3, 0.25), less the PUMA shifts this model omits.

Parameterization

There is no intercept: the cutpoints carry the location (a constant in the formula is collinear with them and dropped). The cutpoints are the last coefficients, cut1 … cut{J-1}, estimated directly rather than as increments, so their standard errors are reported as they are by Stata’s ologit and R’s polr. statsmodels’ OrderedModel estimates the first cutpoint and the logs of the increments instead; its slopes, and their standard errors, are the same as these, and its cutpoints transform into these.

A Newton step that would put the cutpoints out of order is halved until it does not; the log-likelihood is concave in the slopes and cutpoints together, so this is safe.

Levels

The outcome’s levels are ordered as levels=, or else sorted: integers by value, an Enum by its categories, strings alphabetically (which is rarely what an ordered outcome wants: give levels= or make it an Enum). education is an Enum; as strings it needs levels= for the same fit:

strings = people.with_columns(pl.col("education").cast(pl.String))
ologit_stream(f"education ~ {X}", strings, levels=EDUCATION)

Residuals

resid() gives, for each row of the estimation sample, the probability of each level (prob_<level>) and two residuals for the outcome as a whole:

• resid_prob, Li and Shepherd’s (2012) probability-scale residual, P(Y < y) - P(Y > y): in [-1, 1], 0 on average under the model, and the same as statsmodels’ OrderedModel reports (resid_prob), which the test suite checks;

• resid_gen, the generalized residual (Gourieroux et al. 1987), the derivative of the row’s log-likelihood in its linear index: its products with the covariates sum to zero at the estimates, as the score does.

fit.resid().collect()

┌────────┬──────────────┬──────────────────┬────────────┬───────────┐
│ row_id ┆ education    ┆ prob_High school ┆ resid_prob ┆ resid_gen │
╞════════╪══════════════╪══════════════════╪════════════╪═══════════╡
│ 0      ┆ Less than HS ┆ 0.276661         ┆ -0.849079  ┆ -0.849079 │
│ 1      ┆ Some college ┆ 0.038366         ┆ -0.816462  ┆ -0.816462 │
│ 2      ┆ Bachelor's   ┆ 0.048776         ┆ -0.427149  ┆ -0.427149 │
│ 3      ┆ Bachelor's   ┆ 0.228077         ┆ 0.422086   ┆ 0.422086  │
└────────┴──────────────┴──────────────────┴────────────┴───────────┘

For the ordered logit the two are the same. With \(a = F(\kappa_{j-1} - x'\beta)\) and \(b = F(\kappa_j - x'\beta)\) for a row at level \(j\), the probability-scale residual is \(a - (1 - b)\), and since the logistic density is \(F(1 - F)\), the generalized residual is \((a(1-a) - b(1-b)) / (b - a) = a + b - 1\) too. For the ordered probit they differ (the generalized residual is unbounded); the script checks both:

logit: largest |resid_prob - resid_gen| 7.88e-15
probit: largest |resid_prob - resid_gen| 2.15e+00

Predictions

predict(type="response") (the default) returns each level’s probability, type="class" the most probable level, type="link" the index \(x'\beta\) (with the fixed effect, if any):

fit.predict(new).collect()
fit.predict(new, type="class").collect()

┌───────────┬───────────────────┬────────────────┬───────────────────┬─────────────────┬───────────────┐
│ person_id ┆ prob_Less than HS ┆ prob_High      ┆ prob_Some         ┆ prob_Bachelor's ┆ prob_Graduate │
│           ┆                   ┆ school         ┆ college           ┆                 ┆               │
╞═══════════╪═══════════════════╪════════════════╪═══════════════════╪═════════════════╪═══════════════╡
│ 0         ┆ 0.109424          ┆ 0.231098       ┆ 0.254371          ┆ 0.241318        ┆ 0.16379       │
│ 1         ┆ 0.028216          ┆ 0.080533       ┆ 0.148869          ┆ 0.289165        ┆ 0.453218      │
│ 2         ┆ 0.026633          ┆ 0.076496       ┆ 0.143304          ┆ 0.285609        ┆ 0.467958      │
└───────────┴───────────────────┴────────────────┴───────────────────┴─────────────────┴───────────────┘

The ordered probit

oprobit_stream is the same model with normal errors, and takes the same options except random intercepts and firth=. Its coefficients are the logit’s divided by about 1.7, the ratio of the two distributions’ scales, and it fits these (logistic) data slightly worse:

probit = oprobit_stream(f"education ~ {X}", people)

┌─────────────┬───────────┬──────────┬──────────┐
│ coefficient ┆ logit     ┆ probit   ┆ ratio    │
╞═════════════╪═══════════╪══════════╪══════════╡
│ age10       ┆ 0.045762  ┆ 0.02662  ┆ 1.719099 │
│ female      ┆ 0.242993  ┆ 0.140909 ┆ 1.724473 │
│ log_income  ┆ 0.811634  ┆ 0.475319 ┆ 1.707556 │
│ ...
log-likelihood: logit -300550.7, probit -300686.3

Fixed effects

With an absorbed fixed effect ("education ~ age10 + female | puma"), each PUMA’s effect shifts its rows’ latent index, and cut1 is normalized to 0 (the effects carry the location; cutpoints reports the 0), so the table starts at cut2. See fixed effects, and for random PUMA shifts, random intercepts.

Everything else – standard errors, weights, l2, the options – is as for the multinomial logit.

The full script: examples/ordered.py

"""The ordered logit and probit: cutpoints, residuals, predictions and a
fixed effect.

    python examples/ordered.py
"""

import polars as pl

from logit_stream import ologit_stream, oprobit_stream
from logit_stream.simulate import EDUCATION, simulate_people, truth_people

people = simulate_people(200_000, seed=1)
X = "age10 + female + log_income + owner + urban"

# education is a Polars Enum, so its levels are in the Enum's order
print(people["education"].dtype)

# ---- the ordered logit: slopes, then the cutpoints cut1 ... cut4
fit = ologit_stream(f"education ~ {X}", people, vcov="hetero", verbose=False)
print(fit.tidy())
print("cutpoints:", fit.cutpoints.round(3))

# simulate_people centres log income at 10.4 and adds PUMA shifts, which
# this model omits; the slopes are close to what it draws from
truth = truth_people()["education"]
print({k: (round(fit.coef()[k], 3), v) for k, v in truth.items() if not k.startswith("cut")})

# ---- a string outcome needs levels= (sorted alphabetically otherwise)
strings = people.with_columns(pl.col("education").cast(pl.String))
again = ologit_stream(f"education ~ {X}", strings, levels=EDUCATION, verbose=False)
print("same fit from strings with levels=:", bool(abs(again.beta - fit.beta).max() < 1e-10))

# ---- residuals: each level's probability, the probability-scale residual
# and the generalized residual
r = fit.resid().collect()
print(r.select("row_id", "education", "prob_High school", "resid_prob", "resid_gen").head(4))
print(r.select(pl.col("resid_prob", "resid_gen").mean()))

# ---- predictions for new data
new = simulate_people(1_000, seed=2)
p = fit.predict(new).collect()
print(p.select("person_id", pl.col("^prob_.*$")).head(3))
print(fit.predict(new, type="class").select("person_id", "education", "predicted")
         .head(3).collect())
print(fit.predict(new, type="link").select("person_id", "xb").head(3).collect())

# ---- the ordered probit: the same model with normal errors; coefficients
# about 1/1.7 of the logit's
probit = oprobit_stream(f"education ~ {X}", people, verbose=False)
print(pl.DataFrame({"coefficient": fit.coefnames, "logit": fit.beta, "probit": probit.beta,
                    "ratio": fit.beta / probit.beta}))
print(f"log-likelihood: logit {fit.loglik:.1f}, probit {probit.loglik:.1f}")

# the two residuals are the same for the logit (both are F(cut_{j-1} - xb) +
# F(cut_j - xb) - 1), not for the probit
for name, res in (("logit", r), ("probit", probit.resid().collect())):
    d = (res["resid_prob"] - res["resid_gen"]).abs().max()
    print(f"{name}: largest |resid_prob - resid_gen| {d:.2e}")

# ---- an absorbed fixed effect: each PUMA shifts the latent index, and
# cut1 is normalized to 0
fe = ologit_stream(f"education ~ {X} | puma", people, vcov={"CRV1": "puma"},
                   verbose=False)
print(fe.tidy())
print(fe.fixef().head(3).collect())