4.5 Ordered outcomes

Ratings such as “disagree, neutral, agree” are ordered, but the steps between them need not be equal. The ordered logistic model treats them as cut points on a hidden continuous scale, and lets predictors shift people along that scale.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from statsmodels.miscmodels.ordinal_model import OrderedModel

from statly import style

style.use()
rng = np.random.default_rng(1)

anes = sm.datasets.anes96.load_pandas().data  # 944 American voters, 1996
PARTY = [
    "strong Dem.",
    "weak Dem.",
    "lean Dem.",
    "independent",
    "lean Rep.",
    "weak Rep.",
    "strong Rep.",
]
party = anes.PID.astype(int).to_numpy()  # party identification, 0 to 6
ideology = (
    anes.selfLR.to_numpy()
)  # self-placement from 1 (extremely liberal) to 7 (extremely conservative)
CATEGORY_COLORS = plt.cm.RdBu_r(np.linspace(0.05, 0.95, 7))

4.5.1 Ordered categories

Many outcomes are categories with a natural order: pain from none to severe, a tumour grade, a survey answer from “strongly disagree” to “strongly agree”, a party identification from strong Democrat to strong Republican. Treating the codes 0 to 6 as numbers in a linear regression assumes that every step is the same size, which nobody claims. Treating them as unordered categories (multinomial regression, Chapter 4.4) throws away the order and needs six sets of coefficients.

counts = np.bincount(party, minlength=7)

fig, axes = style.subplots(2, height=3.3, width=10.5)
axes[0].barh(range(7), counts, color=CATEGORY_COLORS)
axes[0].set_yticks(range(7), PARTY)
axes[0].set(xlabel="voters")
table = pd.crosstab(ideology, party, normalize="index").to_numpy()
bottom = np.zeros(table.shape[0])
levels = np.sort(np.unique(ideology))
for j in range(7):
    axes[1].bar(
        levels, table[:, j], bottom=bottom, color=CATEGORY_COLORS[j], width=0.8, label=PARTY[j]
    )
    bottom += table[:, j]
axes[1].yaxis.set_major_formatter(style.PERCENT)
axes[1].legend(fontsize=7, bbox_to_anchor=(1.02, 1), loc="upper left")
axes[1].set(xlabel="self-placement (1 liberal … 7 conservative)", ylabel="share of voters")
fig.tight_layout()
plt.show()

Figure 4.5.1: Party identification of 944 American voters in 1996 (left) and its make-up by the voters’ own placement on a liberal (1) to conservative (7) scale (right, each bar the voters with that placement). The more conservative the self-placement, the more Republican the party identification, but the categories shift gradually, not in equal steps.

4.5.2 A hidden scale with cut points

The ordered logistic (or proportional odds) model imagines a continuous hidden score $Y^*$, here a voter’s leaning, equal to $\beta x$ plus logistic noise. The observed category is the interval of $Y^*$ between two cut points $\theta_{j-1} < Y^* \le \theta_j$. In terms of probabilities:

$$P(Y \le j \mid x) = \frac{1}{1 + e^{-(\theta_j - \beta x)}}, \qquad j = 0, 1, \dots, 5$$

One slope β serves all the cut points: a predictor slides the whole distribution of the hidden score up or down, and with it the chance of every category.

model = OrderedModel(party, ideology[:, None], distr="logit").fit(method="bfgs", disp=False)
slope = model.params[0]
cuts = model.model.transform_threshold_params(model.params)[1:-1]
latent = np.linspace(-3, 13, 500)
SHORT = [
    "SD",
    "WD",
    "LD",
    "I",
    "LR",
    "WR",
    "SR",
]  # strong, weak, leaning Democrat; independent; …

fig, ax = plt.subplots(figsize=(9, 3.4))
for x_value, color in ((1, style.BLUE), (4, style.GREY), (7, style.CORAL)):
    ax.plot(
        latent,
        stats.logistic.pdf(latent, loc=slope * x_value),
        color=color,
        lw=2,
        label=f"self-placement {x_value}",
    )
for c in cuts:
    ax.axvline(c, color=style.INK, lw=1, ls="--")
edges = np.r_[latent[0], cuts, latent[-1]]
for j in range(7):
    height = 0.27 if j % 2 == 0 else 0.22  # alternate heights so narrow bands stay readable
    ax.text(
        (edges[j] + edges[j + 1]) / 2,
        height,
        SHORT[j],
        ha="center",
        fontsize=8,
        color=style.INK,
    )
ax.legend(fontsize=8, loc="center right")
ax.set(xlabel="hidden score", ylabel="how common", yticks=[], ylim=(0, 0.31))
plt.show()

Figure 4.5.2: The ordered logistic model fitted to party identification with self-placement as the predictor. Each curve is the distribution of the hidden score for voters at one self-placement (1, 4 and 7); the vertical lines are the six fitted cut points, which split the hidden scale into the seven categories (SD: strong Democrat, WD: weak, LD: leaning, I: independent, LR, WR, SR: the same for Republican). Moving to the right on the ideology scale slides the whole curve right, shifting probability from the Democratic categories to the Republican ones. The cut points are not evenly spaced.

4.5.3 Reading the model

The slope is a log odds ratio, as in logistic regression, with a twist: it applies to every split of the categories into “this or lower” against “higher”. One step to the right on the ideology scale multiplies the odds of being more Republican, at any cut, by the same factor:

X = anes[["selfLR", "age", "income"]]
full = OrderedModel(party, X, distr="logit").fit(method="bfgs", disp=False)
rows = []
for name, label in (
    ("selfLR", "self-placement (per step)"),
    ("age", "age (per year)"),
    ("income", "income (per bracket)"),
):
    b, se, p = full.params[name], full.bse[name], full.pvalues[name]
    rows.append(f"| {label} | {b:.3f} | {se:.3f} | {np.exp(b):.2f} | {p:.3g} |")
print("| predictor | β | standard error | odds ratio | p-value |")
print("|:--|--:|--:|--:|--:|")
print("\n".join(rows))

Table 4.5.1: The ordered logistic model of party identification on self-placement, age and income (income in the survey’s 24 brackets).

predictor	β	standard error	odds ratio	p-value
self-placement (per step)	0.994	0.052	2.70	2.96e-80
age (per year)	-0.006	0.004	0.99	0.0789
income (per bracket)	0.066	0.010	1.07	6.17e-11

Each step to the conservative side multiplies the odds of a more Republican identification by about 2.7, holding age and income fixed. Higher income leans Republican too; age, here, hardly at all. Probabilities are easier to read than odds:

grid = np.linspace(1, 7, 61)
newdata = pd.DataFrame({"selfLR": grid, "age": anes.age.mean(), "income": anes.income.mean()})
probabilities = full.predict(newdata)

fig, ax = plt.subplots(figsize=(8.5, 3.4))
ax.stackplot(grid, probabilities.T, colors=CATEGORY_COLORS, labels=PARTY, alpha=0.9)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=7, bbox_to_anchor=(1.02, 1), loc="upper left")
ax.set(
    xlabel="self-placement (1 liberal … 7 conservative)",
    ylabel="predicted probability",
    xlim=(1, 7),
    ylim=(0, 1),
)
plt.show()

Figure 4.5.3: Predicted probability of each party identification across the self-placement scale, at average age and income, stacked so that each vertical slice adds to 100%. The bands for the strong identifications grow at the ends; the middle categories peak in between.

4.5.4 Checking proportional odds

The model’s key assumption is that one slope fits all the cuts. A direct check: fit a separate binary logistic regression at each cut (“this category or lower” against “higher”) and compare the six slopes. If they are similar, the assumption is reasonable; if they differ a lot, a model with separate slopes (a partial proportional odds model) or a multinomial model is safer.

binary_slopes = []
for cut in range(6):
    higher = (party > cut).astype(float)
    fit = sm.Logit(higher, sm.add_constant(X)).fit(disp=0)
    binary_slopes.append((fit.params["selfLR"], fit.bse["selfLR"]))
binary_slopes = np.array(binary_slopes)

fig, ax = plt.subplots(figsize=(8, 3.2))
ax.errorbar(
    range(6),
    binary_slopes[:, 0],
    yerr=1.96 * binary_slopes[:, 1],
    fmt="o",
    color=style.BLUE,
    capsize=4,
    lw=2,
)
ax.axhline(full.params["selfLR"], color=style.CORAL, ls="--", lw=1.5)
ax.set_xticks(range(6), [f"{PARTY[j]}\nor lower" for j in range(6)], fontsize=7.5)
ax.set(ylabel="slope of self-placement (log odds)")
plt.show()

Figure 4.5.4: The slope of self-placement in six separate logistic regressions, one per cut of party identification (with 95% intervals), against the single slope of the ordered model (dashed). The slopes rise from about 0.75 at the first cut to about 1.2 at the middle and Republican cuts, and the first cut’s interval excludes the single slope. One slope is a compromise here; a model letting the slope vary by cut would be worth comparing.

4.5.5 In one picture

• Ordered categories have an order but not equal steps.
• The ordered logistic model: $P(Y \le j) = \text{logistic}(\theta_j - \beta x)$, cut points $\theta_j$ on a hidden scale, one slope β.
• $e^\beta$ is the odds ratio for a higher category, the same at every cut.
• Show predicted probabilities of each category.
• Check proportional odds with separate logistic regressions per cut.

Next: counts, and the family of generalised linear models.