Ordered logit, ordered probit and multinomial logit in Python

Some outcomes are categories, not quantities. Survey answers from "very unhappy" to "very happy", credit ratings and the highest level of schooling have a natural order, but the distance between neighbouring categories means nothing. Others, such as the way people travel to work or the occupation they choose, have no order at all. Coding the categories 1, 2, 3 and running OLS treats the codes as measurements. Ordered and multinomial models treat them as categories and estimate the probability of each one.

What ordered and multinomial models estimate

An ordered model assumes a continuous latent variable y* = x′β + ε behind the observed category. The outcome is category j when y* falls between two cutpoints, κⱼ₋₁ < y* ≤ κⱼ:

Pr(y = j | x) = F(κⱼ − x′β) − F(κⱼ₋₁ − x′β),   κ₀ = −∞, κ_J = +∞

With a logistic F this is the ordered logit, also called the proportional-odds model (McCullagh 1980). With a standard normal F it is the ordered probit (McKelvey and Zavoina 1975). The J − 1 cutpoints are estimated together with β and replace the intercept. A positive coefficient moves probability towards the higher categories: the top category becomes more likely and the bottom one less likely, while a middle category can gain or lose. Only marginal effects on each probability show which.

The model uses one β for every cutpoint. In the ordered logit, exp(β) is therefore the odds ratio of being above any cutpoint, whichever cutpoint is chosen. This is the proportional-odds, or parallel-regression, assumption, and Brant (1990) proposed a test of it. When it fails, the multinomial logit is an option that ignores the order instead of imposing it.

A multinomial logit handles unordered categories. Each category except a base category has its own coefficient vector, and the base is normalized to zero:

Pr(y = j | x) = exp(x′βⱼ) / Σₗ exp(x′βₗ),   β_base = 0

It is the model implied by random utility with independent type-I extreme value (Gumbel) errors. exp(βⱼ) is a relative-risk ratio: the change in the odds of category j against the base. A coefficient compares a category only with the base, so its sign need not match the effect on that category's probability; marginal effects settle the question.

The multinomial logit assumes the independence of irrelevant alternatives (IIA): the odds between two categories do not depend on which other categories exist. Adding a second bus line that is a near copy of the first should mostly draw passengers from the first bus, but under IIA it draws from every mode in proportion. Hausman and McFadden (1984) give a test. When IIA is implausible, a nested logit or a multinomial probit relaxes it; OpenEconometrics has oe.nlogit and oe.mprobit (two or three alternatives) for data with one row per case and alternative.

In OpenEconometrics

oe.ologit and oe.oprobit take the same arguments. The ordered model in the example corresponds to the Stata command ologit satisfaction lwage i.night, and to oprobit with the same variables:

• oe.ologit(data=df, y="satisfaction", x=["lwage", "night"], categorical=["night"]). A numeric outcome is ordered by value, and the values need not be consecutive. An ordered pandas Categorical is ordered by its categories. An unordered text column is refused with invalid_ordered_outcome.
• There is no constant. The cutpoints are reported after the slopes as /cut1 … /cut{J-1}, with standard errors. offset= adds a column to the index with coefficient 1.

oe.mlogit(data=df, y="mode", x=["income", "distance"], base="car") corresponds to mlogit mode income distance, baseoutcome(1) in Stata, where Stata needs mode coded as numbers with 1 for car. The outcome may hold numbers, strings, booleans or a Categorical. Without base=, the most frequent category is the base. Terms are named <category>:<term>, for example bus:income.

All three share these conventions:

• covariance="nonrobust" (default), "opg", "robust" or "cluster". Giving cluster= with one or two columns selects clustering; two-way clustering is an OpenEconometrics extension.
• weights= with weight_type= ("fweight", "aweight", "pweight", "iweight"). Sampling weights default to robust. Analytic weights are an OpenEconometrics extension: Stata's commands do not accept them.
• tests["model"] is the likelihood-ratio test against the cutpoints-only or constant-only model under nonrobust, and a Wald test otherwise. The pseudo R² is McFadden's.
• A regressor that predicts some outcomes perfectly raises separation_detected. In an ordered model this includes a regressor that separates the two sides of one cutpoint.

oe.margins(fit, "lwage", data=df) returns the average marginal effect on every category probability, one row per outcome, with delta-method standard errors from the full covariance of slopes and cutpoints. outcome=4 restricts it to one category, selected by its actual label rather than its position. oe.predict(fit, new_rows, interval="mean") returns all category probabilities. The discrete-choice guide gives the likelihoods and conventions, and the prediction guide covers the probability and effect calculations.

Supported options

Options of oe.mlogit

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 3,000 workers. Job satisfaction on a four-point scale follows a true ordered logit with latent index 1.0·lwage − 0.6·night plus a logistic error, and cutpoints 1.5, 2.7 and 3.9; log wages are normal with mean 3 and standard deviation 0.6, and half the workers do night shifts. The same workers choose a commuting mode from random utilities with Gumbel errors. Against the car, the bus has utility 1.0 − 0.30·income + 0.05·distance and the bike 2.0 − 0.20·income − 0.40·distance, with income in tens of thousands and distance in kilometres.

import numpy as np
import pandas as pd
import openecon as oe

CUTS = [1.5, 2.7, 3.9]                     # true cutpoints of the ordered model
MODE = {"bus": (1.0, -0.30, 0.05),         # true (intercept, income, distance) against car
        "bike": (2.0, -0.20, -0.40)}


def simulate(n=3000, seed=11):
    rng = np.random.default_rng(seed)
    lwage = rng.normal(3.0, 0.6, size=n)
    night = rng.integers(0, 2, size=n)
    latent = 1.0 * lwage - 0.6 * night + rng.logistic(size=n)
    satisfaction = 1 + np.digitize(latent, CUTS)           # 1 = very unhappy ... 4 = very happy
    income = rng.normal(5.0, 1.5, size=n)                   # tens of thousands
    distance = rng.exponential(5.0, size=n)                 # km to work
    utility = np.column_stack([np.zeros(n)] + [a + b * income + c * distance
                                               for a, b, c in MODE.values()])
    choice = (utility + rng.gumbel(size=(n, 3))).argmax(axis=1)
    mode = np.array(["car", "bus", "bike"])[choice]
    return pd.DataFrame({"satisfaction": satisfaction, "lwage": lwage, "night": night,
                         "mode": mode, "income": income, "distance": distance})


if __name__ == "__main__":
    df = simulate()
    ordered = oe.ologit(data=df, y="satisfaction", x=["lwage", "night"], categorical=["night"])
    print(ordered.summary())
    probit = oe.oprobit(data=df, y="satisfaction", x=["lwage", "night"], categorical=["night"])
    print("oprobit slopes:", {c.term: round(c.estimate, 4) for c in probit.coefficients[:2]})

    ame = oe.margins(ordered, "lwage", data=df)
    print("\nAME of lwage on Pr(satisfaction = j):")
    print(ame[["outcome", "estimate", "std_error", "ci_low", "ci_high"]].round(4).to_string(index=False))

    choice = oe.mlogit(data=df, y="mode", x=["income", "distance"], base="car")
    truth = {f"{mode}:{term}": value for mode, values in MODE.items()
             for term, value in zip(["Intercept", "income", "distance"], values)}
    rows = [(c.term, truth[c.term], c.estimate, c.std_error) for c in choice.coefficients]
    print(f"\nmlogit, base = {choice.extra['base']}, shares:",
          df["mode"].value_counts(normalize=True).round(3).to_dict())
    print(pd.DataFrame(rows, columns=["term", "true", "estimate", "std_error"]).round(4).to_string(index=False))
    test = choice.tests["model"]
    print(f"{test['label']}: chi2({test['df']}) = {test['statistic']:.1f}")

With openecon 0.3.20a1 from PyPI the script prints:

Ordered logistic regression — satisfaction
Observations: 3000  |  Covariance: nonrobust  |  Confidence: 95%

Term             Estimate  Std. error         z     P>|stat|   CI lower   CI upper
--------------  ---------  ----------  --------  -----------  ---------  ---------
[satisfaction]                                                                    
lwage             1.00532   0.0583157   17.2392  1.34832e-66   0.891022    1.11962
night[1]        -0.719069   0.0672868  -10.6866  1.17567e-26  -0.850948  -0.587189
/cut1             1.49805    0.178505   8.39218  4.77191e-17    1.14818    1.84791
/cut2             2.68793     0.18277   14.7066  5.84748e-49    2.32971    3.04615
/cut3             3.92128    0.189518   20.6908  4.18867e-95    3.54983    4.29273

log_likelihood: -3939.21  |  pseudo_r_squared: 0.0526631  |  aic: 7888.42  |  bic: 7918.45
n_categories: 4
LR chi2 test against the cutpoints-only model: chi2(2) = 437.966, p = 7.885e-96
oprobit slopes: {'lwage': 0.6034, 'night[1]': -0.4323}

AME of lwage on Pr(satisfaction = j):
 outcome  estimate  std_error  ci_low  ci_high
       1   -0.1758     0.0097 -0.1948  -0.1568
       2   -0.0489     0.0038 -0.0563  -0.0415
       3    0.0548     0.0040  0.0469   0.0626
       4    0.1699     0.0095  0.1513   0.1885

mlogit, base = car, shares: {'car': 0.388, 'bus': 0.328, 'bike': 0.284}
          term  true  estimate  std_error
bike:Intercept  2.00    1.8765     0.1885
   bike:income -0.20   -0.1704     0.0334
 bike:distance -0.40   -0.4060     0.0225
 bus:Intercept  1.00    0.9295     0.1665
    bus:income -0.30   -0.2929     0.0307
  bus:distance  0.05    0.0579     0.0086
LR chi2 test against the constant-only model: chi2(4) = 899.4

Reading the output

• The ordered logit recovers the latent model. The wage coefficient is 1.0053 (true 1.0) and the cutpoints are 1.498, 2.688 and 3.921 (true 1.5, 2.7 and 3.9). The night-shift coefficient of −0.7191 is 1.8 standard errors from the true −0.6, and its interval from −0.851 to −0.587 just contains it.
• Ordered probit coefficients are on a different scale. The oprobit slopes are 0.6034 and −0.4323; the logit slopes are about 1.66 times as large, because the logistic error has a larger standard deviation than the normal one.
• Marginal effects show where the probability moves. The four AMEs of log wage sum to zero. A rise of 0.1 in log wage, roughly 10% higher pay, lowers the probability of the lowest category by about 1.8 percentage points (−0.1758 × 0.1) and raises that of the highest by about 1.7 points. The second category also loses (−0.0489): a middle category can move either way.
• The multinomial logit recovers all six coefficients. Each estimate is within one standard error of its true value; for example, each kilometre lowers the log odds of cycling against driving by 0.4060 (true 0.40). The bus coefficient on distance, 0.0579, is positive, so the odds of the bus against the car rise slightly with distance.
• Both models fit far better than their null models, which contain only the cutpoints or the constants: the likelihood-ratio statistics are 438.0 and 899.4.

Limitations

• Cluster-robust covariance for oe.ologit allows at most 2 clustering dimensions.
• Cluster-robust covariance for oe.oprobit allows at most 2 clustering dimensions.
• Cluster-robust covariance for oe.mlogit allows at most 2 clustering dimensions.
• No Brant test of proportional odds, no Hausman–McFadden test of IIA and no generalized ordered model with category-specific slopes are provided.
• oe.mlogit takes no offset, and none of the three accepts linear constraints on the coefficients.
• The prediction adapters return category probabilities and their marginal effects; classification, scores and tests of index differences between pairs of categories are outside them.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Brant, R. (1990). Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics, 46(4), 1171–1178. doi:10.2307/2532457
• Hausman, J. and McFadden, D. (1984). Specification tests for the multinomial logit model. Econometrica, 52(5), 1219–1240. doi:10.2307/1910997
• McCullagh, P. (1980). Regression models for ordinal data. Journal of the Royal Statistical Society Series B: Statistical Methodology, 42(2), 109–127. doi:10.1111/j.2517-6161.1980.tb01109.x
• McKelvey, R. D. and Zavoina, W. (1975). A statistical model for the analysis of ordinal level dependent variables. The Journal of Mathematical Sociology, 4(1), 103–120. doi:10.1080/0022250X.1975.9989847