How to compute marginal effects for multinomial and ordered outcome models

Prerequisites

Problem statement

You have a multinomial logit (MNLogit) or ordered logit/probit (OrderedModel) and you want marginal effects for each outcome class probability. Every statistic returns K values (one per class) with full joint covariance across both variables and classes.

Minimal working solution

Margins auto-detects multi-outcome models. Use outcome= to subset to specific classes.

import numpy as np
import pandas as pd
import statsmodels.api as sm
from smmargins import Margins

rng = np.random.default_rng(7)
n = 500
df = pd.DataFrame({
    "x1": rng.standard_normal(n),
    "x2": rng.standard_normal(n),
})
# 3-class outcome
z0 = 0.5 + 1.0 * df["x1"] - 0.5 * df["x2"] + 0.3 * rng.standard_normal(n)
z1 = -0.3 + 0.5 * df["x1"] + 0.5 * df["x2"] + 0.3 * rng.standard_normal(n)
df["y"] = np.argmax(np.column_stack([np.zeros(n), z0, z1]), axis=1)

# Fit multinomial logit
X = sm.add_constant(df[["x1", "x2"]])
fit_mn = sm.MNLogit(df["y"], X).fit(disp=False)

M = Margins(fit_mn)
print(f"Number of outcomes: {M.n_outcomes}")
print(f"Outcome labels: {M.outcome_labels}")

# AAP for all classes (probabilities sum to 1)
aap = M.predict()
print(f"\nAAP (sum = {aap.estimate.sum():.6f}):")
print(aap.summary())

# AME of x1 for all outcome classes
ame = M.dydx("x1")
print("\nAME of x1 per outcome class:")
print(ame.summary())

Variations

Subset to specific outcomes

# By index
print("Outcome 1 only:")
print(M.predict(outcome=1))

# By label
print("\nOutcomes 0 and 2:")
print(M.predict(outcome=[0, 2]))

# Post-hoc slicing of an existing result
full = M.dydx("x1")
sub = full.outcome(1)
print("\nPost-hoc slice for outcome 1:")
print(sub.summary())

Multiple variables jointly

# 2 variables * 3 outcomes = 6 rows; full joint vcov
joint = M.dydx(["x1", "x2"])
print(f"Shape: {joint.estimate.shape}")
print(f"Vcov shape: {joint.vcov.shape}")

Elasticities on multi-outcome models

# Elasticities work on multi-outcome models too
eyex = M.dydx("x1", method="eyex")
print(eyex.summary())

⚠️ Trade-off: MNLogit uses an analytic softmax-derivative (no extra forward predict calls), making it fast. OrderedModel uses central finite differences for the cumulative-link parameterization, which is slower. Elasticities on multi-outcome models emit a RuntimeWarning when any predicted class probability falls below 1e-12 (numerically unstable).

When to use this

Use multi-outcome margins when your dependent variable has more than two unordered categories (MNLogit) or ordered categories (OrderedModel). The default behavior returns all K outcome classes — subset with outcome= when you only care about specific classes. The joint covariance across classes lets you form cross-outcome contrasts (e.g., “effect on P(class=1) minus effect on P(class=0)”).

When NOT to use this

⚠️ Trade-off: Do not treat the K AMEs as independent — they share the same parameter covariance and are constrained to sum to zero across outcomes for any given variable (in MNLogit). Single-outcome models silently ignore outcome=, so the same code paths work uniformly, but be careful not to assume outcome= does something on a binary logit.

See also