Ordinal and Multinomial Regression

The missing regression family

Every dichotomous and continuous outcome in this catalog has a regression home — logistic for binary, Cox for time-to-event, gamma/Poisson for counts and costs — but a very common RWE outcome shape has none: a fixed-window endpoint with three or more discrete levels. Functional-outcome scales (modified Rankin Scale, mRS, 0-6), performance-status grades (ECOG 0-4), adverse-event severity grades (CTCAE 1-5), patient-satisfaction Likert items, and discharge disposition (home / skilled nursing facility / hospice-death) are all routinely dichotomized in practice ("mRS 0-2 vs 3-6", "any grade 3+ AE") purely because logistic regression is the tool analysts reach for by default. Ordinal and multinomial regression are the correct, information-preserving alternatives, and both are direct extensions of the same maximum-likelihood machinery as binary logistic regression.

The proportional-odds (cumulative logit) model

For an ordered outcome Y with K categories (1 = best ... K = worst), the model specifies K-1 cumulative logits, one at each cutpoint between adjacent categories:

logit[P(Y <= k | X)] = alpha_k - beta'X, for k = 1, ..., K-1

This is the parameterization used by the standard software (R's MASS::polr, Python's statsmodels OrderedModel — confirmed by fitting both directly). Each cutpoint gets its own intercept (threshold) alpha_k, but critically, the same covariate vector beta is shared across every cutpoint — this is the proportional-odds assumption. Because of the minus sign, a positive beta_j lowers P(Y<=k) at every cutpoint, i.e. shifts probability mass toward the higher, worse categories: exp(beta_j) is therefore the common odds ratio for a worse outcome per one-unit increase in X_j, and its reciprocal, exp(-beta_j), is the common odds ratio for a better outcome. These are two directions of the same single fitted number, so which one is being reported must always be stated explicitly — this entry reports exp(-beta_j), the better-outcome direction, throughout (see the worked example). Software sign conventions are not universal: SAS PROC LOGISTIC is commonly documented with the opposite-sign form logit[P(Y<=k)] = alpha_k + beta'X, so a coefficient carried across R/Python and SAS without re-checking its sign can silently flip which direction "better" points — always confirm the direction against a manual cross-tab or the software's own printed odds-ratio labels before reporting. Whichever direction is chosen, the same beta_j applies at every cutpoint simultaneously — the effect on the odds is the same no matter where the scale is split into better-versus-worse: the odds ratio for "good vs moderate-or-poor" equals the odds ratio for "good-or-moderate vs poor." This is a constraint on the cumulative cutpoints, not on adjacent-category (poor->moderate vs moderate->good) steps, which need not be equal even when proportional odds holds — in this entry's worked example the cumulative ORs are both 2.0 while the adjacent-category ORs are 1.6 (good vs moderate) and 1.5 (moderate vs poor). This is what lets an ordinal endpoint be summarized by one number instead of K-1 separate binary-logistic odds ratios, one per arbitrary dichotomization, each with its own multiplicity problem.

Testing the proportional-odds assumption — the Brant test

The proportional-odds assumption is testable, not free. Brant's (1990) test fits the K-1 separate binary logistic models implied by each cutpoint, then checks whether their coefficient vectors are statistically distinguishable from the single shared beta the proportional-odds model estimates, using a Wald chi-square test (overall, and per-covariate). A large, significant Brant statistic means the constant-effect assumption is violated for at least one covariate — commonly the treatment arm itself, if a drug shifts patients out of the worst category more than it shifts them between the two best categories. Two caveats analysts routinely mishandle: the Brant test's power scales with sample size, so with tens of thousands of claims records even a clinically trivial deviation from proportionality can be "significant," while in a modest registry sample a real violation may go undetected; and with many categories or covariates, running the test on every covariate individually inflates the false positive rate for at least one violation, so a global test should be checked first. The Brant test is one diagnostic, not a certification: a non-significant result does not prove proportional odds holds — it may simply reflect limited power — and it should be paired with cutpoint-specific estimates and predicted-probability plots rather than trusted alone. The exposure/treatment variable should be assessed specifically and pre-specified as such regardless of what the global test shows, because a clinically important departure confined to that one covariate can be masked within an aggregate global statistic.

When proportional odds fails: partial and generalized ordered models

If the Brant test flags a specific covariate (very commonly the treatment/exposure indicator, which is exactly the covariate an RWE study cares most about getting right), the fix is not to abandon the ordinal structure and dichotomize — it is a partial proportional-odds model (Peterson & Harrell, 1990), which lets the flagged covariate(s) have cutpoint-specific slopes while holding the rest of the covariates to the shared beta. The fully general version, sometimes called the generalized ordered logit model, frees every covariate's slope at every cutpoint (fit via `VGAM::vglm` in R with `parallel = FALSE`, or Stata's user-written `gologit2`); this sacrifices the one-number summary but is honest about a treatment effect that concentrates at one end of the scale (e.g., a therapy that mostly prevents death but does little to convert "moderate disability" into "no disability").

Multinomial (baseline-category) logistic regression for unordered outcomes

When the categories have no meaningful order — which of three or more index treatments a clinician selected, the cause-specific reason for a hospital discharge (home / SNF-rehab / hospice-or-death), the drug class chosen after a first-line failure — the proportional-odds structure does not apply at all. Multinomial logistic regression instead picks one category as the reference and fits a separate logit for every other category against it:

log[P(Y = j | X) / P(Y = ref | X)] = alpha_j + beta_j'X, for each non-reference category j

exp(beta_j) is a relative risk ratio (RRR) — technically still an odds ratio, comparing category j to the reference category — and each non-reference category gets its own, entirely independent coefficient vector. Interpretation is pairwise: "a one-unit increase in age multiplies the odds of choosing Drug B over Drug A by exp(beta_B), holding the other covariates fixed." That same one-unit increase in age simultaneously changes the odds of choosing Drug C over Drug A too, through Drug C's own, independently estimated coefficient beta_C — the model holds the other covariates fixed, not the other category contrasts, which can move at the same time, each according to its own coefficient. The model rests on the independence of irrelevant alternatives (IIA) assumption — the relative odds between any two categories do not depend on what other categories are available or how attractive they are. IIA is frequently implausible when one "category" is really a close substitute for another (e.g., two brand formulations of the same molecule). A drop-one-category re-fit with a comparison of the remaining coefficients (a Hausman-McFadden-style check) is a limited, sometimes unstable diagnostic — passing it does not make a structurally implausible choice set acceptable. When categories are known close substitutes, prefer a model that does not assume IIA (nested logit, multinomial probit, or a mixed/random-parameters logit) rather than defending multinomial logistic on a drop-one check alone; pre-specifying multinomial logistic as an accepted simplification is standard RWE channeling practice only when the treatment options are not close substitutes for one another.

Interpreting the common OR — connection to the win ratio

The proportional-odds common OR and the win ratio (win-ratio-generalized-pairwise-comparisons-rwe) solve a similar problem — summarizing an ordered comparison with a single number — from different directions. The common OR is a model-based, covariate-adjusted parameter: exp(beta) from a regression that assumes proportional odds and lets you control for confounders directly. The win ratio is a design-based, pairwise statistic: for every possible pair of one treated and one comparator patient, count who "wins" on a pre-specified hierarchy of outcome components (which may include an ordinal scale as one tier), with ties handled explicitly, and report the ratio of wins to losses. The win ratio needs no proportional-odds assumption and extends naturally to a hierarchical composite of different outcome types (death, then hospitalization, then a symptom score), whereas the common OR is confined to a single already-graded ordinal scale but integrates far more easily into a standard covariate-adjusted regression, g-computation, or propensity-weighted framework. In practice, some ordinal-endpoint trials and RWE analyses (stroke, COVID-19) report both a common OR from a proportional-odds model and a win ratio or win odds as complementary summaries of the same shift across the scale, though neither is a required companion to the other.

RWE uses

(1) Ordinal clinical endpoints. The 2007 wave of stroke-trial methodology work (Saver 2007; Bath et al. 2007) established shift analysis of the mRS via proportional-odds regression as more statistically efficient than dichotomizing at any single cutpoint. This approach was subsequently adopted by COVID-19 therapeutic trials that used a 7- to 11-point ordinal clinical-status scale (ambulatory, hospitalized-no-oxygen, hospitalized-oxygen, ICU, death) as either the primary endpoint or, as in the NIAID ACTT-1 remdesivir trial, a pre-specified key secondary analysis of day-15 clinical status run alongside a time-to-recovery primary endpoint (Beigel et al. 2020). EHR- and registry-based RWE emulations and post-authorization safety/effectiveness studies of COVID-19 therapeutics have reused the same ordinal scale and proportional-odds machinery where the underlying clinical-status detail is directly observed; claims-only studies generally cannot construct this scale directly, because it has no standard claims code (see Data-source operational depth below). (2) Channeling and treatment-choice models. When patients are routed to one of three or more treatment options based on measured (and unmeasured) characteristics — the defining feature of confounding-by-indication-channeling — a multinomial logistic regression of the treatment-choice variable on baseline covariates is the direct regression analogue of a binary propensity-score model, generalized to a generalized propensity score (GPS), which then feeds inverse-probability weighting or matching across 3+ treatment arms (Feng et al. 2012; McCaffrey et al. 2013 for a gradient-boosting alternative to the multinomial-logit GPS model). Multinomial logistic is also the standard model for a cause-specific unordered categorical outcome such as discharge disposition.

Pros, cons, and trade-offs

• vs binary logistic regression (logistic-regression-for-binary-outcomes): Ordinal regression retains all the information in a 3+ level scale instead of forcing an often-arbitrary dichotomization, and yields one common OR instead of several binary-cut odds ratios each carrying its own multiplicity problem. It adds the proportional-odds assumption as a genuine structural requirement that must be tested; binary logistic has no analogous assumption to fail. Prefer ordinal when the scale genuinely has 3+ ordered levels and no single cutpoint is clinically privileged; prefer binary logistic when one specific threshold (e.g., "responder at week 12") is the pre-specified clinical question, or when the Brant test badly fails and no partial-PO tooling is available.
• vs the win ratio (win-ratio-generalized-pairwise-comparisons-rwe): The common OR is model-based and integrates covariate adjustment natively; the win ratio is design-based, needs no proportionality assumption, and handles a hierarchical composite of different outcome types directly. Prefer proportional-odds regression for a single, already-graded ordinal scale where covariate-adjusted inference is the goal; prefer the win ratio for a composite hierarchy of distinct clinical events or when avoiding the proportional-odds assumption entirely is a priority.
• vs treating the ordinal score as continuous and using OLS (ols-linear-regression): OLS on the raw numeric scale score assumes the categories are equally spaced (that mRS 0->1 and mRS 4->5 are the same size of clinical change), which is rarely defensible for a clinical grading scale with a small number of categories. Ordinal regression makes no such assumption. OLS is most useful as a quick, easily communicated sensitivity check; it is a defensible primary model only when the scale has many categories, equal spacing is supported by prior evidence, and the estimand and SAP justify a mean-difference interpretation.

When NOT to use

Do not force a proportional-odds model onto an outcome with no defensible order (treatment choice, cause of death category) — use multinomial logistic instead. Do not report a single common OR from a proportional-odds model without first assessing proportionality (the Brant test plus cutpoint-specific estimates and predicted-probability plots) for the exposure of interest specifically; if that assessment shows a material departure — not merely a statistically significant global p-value — report a partial-proportional-odds or generalized ordered-logit result instead of a misleading single number. Do not use multinomial logistic when categories are close substitutes for one another and IIA is implausible; a drop-one-category re-fit is a limited diagnostic, not sufficient justification on its own — prefer an IIA-free model (nested logit, multinomial probit, mixed logit) when the choice set has structurally close substitutes. As with binary logistic, treat differential or incomplete follow-up across arms as a threat to validity, not a default reason to restrict to completers: define the estimand and ascertainment window prospectively, check whether the loss is informative (differs by arm or by outcome severity), and prefer inverse-probability-of- censoring weighting, multiple imputation, or a linkage-based outcome-recovery strategy over a bare complete-case restriction, which should be a labeled sensitivity analysis rather than the default approach. When the loss itself is the clinical event of interest (e.g., death before assessment), move to a time-to-event or competing-risks framework (cox-ph-regression; competing-risks-cause-specific-fine-gray-rwe) instead of collapsing it into a fixed-window outcome.

Interpreting the output

Consider a proportional-odds model for 90-day post-stroke functional outcome (three ordered categories: good, moderate, poor) comparing Drug A (arm=1) vs Drug B (arm=0), fit as `outcome_ord ~ arm` under the alpha_k - beta'X parameterization above with no other covariates. The fitted treatment coefficient is beta_arm = -0.693 (negative, because Drug A shifts probability mass toward the better categories under this parameterization) — confirmed by fitting both MASS::polr in R and statsmodels OrderedModel in Python directly on this dataset, both returning beta_arm = -0.693. exp(-beta_arm) = exp(0.693) = 2.0 is therefore the common odds ratio for a better outcome, identical at both cutpoints of the scale (the worked example below reproduces this exactly from raw counts).

Formal interpretation: patients on Drug A had twice the unadjusted (crude, treatment-only) odds, at every cutpoint of the 90-day outcome scale, of being in a better outcome category than patients on Drug B. This example fits `arm` alone with no covariates, so it is a crude comparison, not an adjusted effect — an adjusted analysis would add baseline covariates to X, and the arm coefficient would then need its own Brant assessment. This is the common odds ratio because it is constrained to be identical whether the cutpoint is "good vs moderate-or-worse" or "good-or-moderate vs poor" — a claim that is only credible because, in this constructed example, both cutpoint-specific odds ratios independently equal 2.0. In real data the Brant test is one diagnostic for how far apart the cutpoint-specific odds ratios have drifted; it does not by itself certify that a single-number summary is adequate.

Practical interpretation: for clinical or HTA communication, the common OR is usually supplemented with the model's predicted probability of each outcome category by arm (a form of standardization analogous to the marginal risk difference reported for binary logistic models) — e.g., "Drug A patients have a predicted 50% chance of a good outcome vs 33% for Drug B" — because a single OR, however well-supported, does not by itself convey the absolute size of the shift in any one category that a clinician or payer cares about. Regulatory relevance here is concrete rather than aspirational: the ACTT-1 remdesivir trial (Beigel et al. 2020), whose secondary endpoints include the ordinal day-15 clinical-status analysis described above, was among the trials FDA reviewed for the drug's emergency use authorization and later full approval; HTA bodies are the primary audience for the predicted-probability-by-category output described above, and are the setting where the common OR vs win-ratio choice discussed earlier in this entry is most often actually made.

Data-source operational depth

• Claims (FFS vs MA vs commercial): Ordinal clinical scales (mRS, ECOG, CTCAE grade) have no standard claims code and are essentially unobservable in pure claims data; claims-only ordinal analyses are largely restricted to administratively derived categories such as discharge disposition (a genuinely unordered multinomial outcome built from the discharge-status code on the index inpatient claim: home, home-health, SNF/IRF, hospice, in-hospital death) or a severity-grade proxy built from a validated diagnosis/procedure algorithm. Multinomial treatment-choice models (channeling) are natural in claims: the exposure variable itself (which of 3+ drugs was dispensed at the index fill) is multinomial, built from the dispensed NDC, and used as the outcome of a generalized-propensity-score model fit on pre-index covariates.
• EHR: Structured severity fields (ECOG in oncology flowsheets, NYHA class in cardiology notes) or free-text/NLP-extracted grades are richer than claims but frequently only recorded for a subset of encounters. Because sicker, more closely monitored patients tend to have more complete grading, the observation process itself may be informative — never assume missing-at-random by default, but do not assume MNAR either without investigating. Report capture rates by arm and by baseline severity, and model predictors of being graded at all, before choosing an analysis: multiple imputation or inverse-probability-of-observation weighting when the missingness is plausibly explained by observed covariates (MAR-conditional), complete-case analysis only as a documented sensitivity check (it can be biased under informative observation, not a default remedy), and explicit missing-not-at-random sensitivity analyses (delta-adjusted pattern-mixture or selection models) when the observation process is suspected to depend on the unobserved severity itself.
• Registry: The cleanest and most complete source of a true clinical ordinal scale — mRS in stroke registries, tumor grade/stage in oncology registries, adjudicated CTCAE grading in post-marketing safety registries — and the benchmark for validating a claims/EHR-derived ordinal or multinomial outcome algorithm.
• Linked claims-EHR-registry: The strongest substrate for an ordinal or multinomial RWE analysis: registry-adjudicated grading for outcome fidelity, claims for complete longitudinal exposure and covariate capture, and EHR for granular clinical covariates; reconcile date discrepancies before fixing the outcome-assessment window.

Worked claims/registry example

Question: does 90-day post-stroke functional outcome (mRS collapsed to three ordered categories: good, moderate, poor) differ between new initiators of antithrombotic Drug A vs Drug B, in a registry-linked-claims cohort. (1) Eligibility/time zero: first qualifying dispensed fill with 365 days of continuous enrollment beforehand and no prior fill of either drug (incident users); assign `arm` from the dispensed NDC. (2) Outcome: the 90-day mRS grade from the linked stroke registry, collapsed to good/moderate/poor. Define the assessment window and estimand prospectively; check whether 90-day assessment completion differs by arm or by early severity (a signal of informative loss) and prefer inverse-probability-of-censoring weighting or multiple imputation for any differential missingness over silently restricting to completers — reserve a complete-case restriction as a labeled sensitivity analysis, not the primary approach. (3) Covariates measured strictly in the half-open pre-index window [index_date - 365 days, index_date): age, sex, NIHSS at admission, prior stroke, comorbidity burden — excluding the index date itself so no post-initiation or same-day, order-ambiguous measurement can leak into the covariate set. (4) Fit the proportional-odds model `outcome_ordinal ~ arm + covariates`; assess proportionality (Brant test plus cutpoint-specific estimates) before reporting a single common OR. (5) If the arm coefficient shows a material Brant departure (common in practice when a drug prevents death much more than it converts moderate to good outcomes), refit as a partial-proportional-odds model that frees only the arm coefficient across cutpoints, and report cutpoint-specific ORs alongside predicted-probability-by-category plots instead of a single potentially misleading common OR.

Decision diagram

flowchart TD
  Q[Categorical outcome with 3+ levels] --> O{Are the categories<br/>naturally ordered?}
  O -->|"Yes - mRS, ECOG,<br/>CTCAE grade, satisfaction scale"| PO[Proportional-odds<br/>cumulative logit model]
  O -->|"No - treatment choice,<br/>discharge disposition"| MN[Multinomial baseline-category<br/>logit model]
  PO --> BR{Brant test:<br/>are cutpoint slopes equal?<br/>(one diagnostic, not a certification)}
  BR -->|"Not significant - PO plausible"| COR[Report single common OR<br/>across all cutpoints]
  BR -->|"Material departure - PO violated"| PPO["Partial / generalized ordered logit -<br/>free the violating covariate's slopes"]
  MN --> IIA{Independence of<br/>irrelevant alternatives plausible?}
  IIA -->|Yes| RRR[Report relative-risk ratios<br/>vs a reference category]
  IIA -->|"No - categories are<br/>close substitutes"| NEST["Prefer an IIA-free model -<br/>nested logit, multinomial probit,<br/>or mixed logit"]

flowchart TD
  Y["Ordered outcome Y with K categories<br/>(e.g., mRS collapsed to good / moderate / poor)"] --> C1["Cutpoint 1:<br/>logit P(Y<=good) = alpha1 - beta'X"]
  Y --> C2["Cutpoint 2:<br/>logit P(Y<=good-or-moderate) = alpha2 - beta'X"]
  C1 --> PROP["Same beta at every cutpoint<br/>(proportional-odds assumption)"]
  C2 --> PROP
  PROP --> OR["exp(-beta) = common OR:<br/>odds of a BETTER outcome for treated vs<br/>comparator, constant across all cutpoints<br/>(exp(+beta) is the reciprocal, WORSE-outcome, OR)"]
  OR --> WR["Conceptually related to the win ratio:<br/>both summarize an ordered comparison with<br/>one number, but the win ratio needs no<br/>proportional-odds assumption and handles<br/>ties/hierarchical outcomes directly"]

Worked example

Scenario

A team compares 90-day functional outcome after stroke for patients starting Drug A vs Drug B. The outcome has three ordered categories — good, moderate, and poor — and the team wants a single number summarizing whether Drug A shifts patients toward better outcomes across the whole scale, using a proportional-odds model. Here are the raw category counts by arm before any model is fit.

Dataset

90-day functional outcome counts by treatment arm (120 patients per arm).

Steps

Result

Common OR = 2.0 (crude, unadjusted, treatment-only): patients on Drug A had twice the odds, at every cutpoint of the 90-day outcome scale, of being in a better outcome category than patients on Drug B; the fitted model's arm coefficient is beta_arm = -0.693, and exp(-beta_arm) = 2.0 is the direction reported here (better outcome). In real data the two cutpoint-specific odds ratios almost never line up this exactly, and an adjusted model would add covariates before this number is reported as a primary result. The Brant test is one diagnostic for how far the cutpoint-specific odds ratios have drifted apart, not a pass/fail certification that proportional odds holds; a material — not merely statistically significant — departure, especially for the exposure variable, signals that a partial-proportional-odds model (which lets the treatment effect differ by cutpoint) should be reported alongside or instead of the single common OR.

Trade-offs

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel

# ── 1. Proportional-odds (cumulative logit) model ────────────────────────────────────────
# cohort: person_id, arm (0/1), outcome_ord ('good' < 'moderate' < 'poor'), covariates...
def fit_proportional_odds(cohort: pd.DataFrame, covariates: list[str],
                          treat: str = "arm", outcome: str = "outcome_ord") -> dict:
    d = cohort.copy()
    d[outcome] = pd.Categorical(d[outcome], categories=["good", "moderate", "poor"],
                                 ordered=True)
    X = d[[treat] + covariates]
    mod = OrderedModel(d[outcome], X, distr="logit")   # no intercept: thresholds serve as alphas
    res = mod.fit(method="bfgs", disp=False)

    # OrderedModel fits logit[P(Y<=k)] = alpha_k - beta'X, so exp(beta) is the OR for a WORSE
    # (higher-category) outcome per unit increase in `treat`; exp(-beta) is the OR for a BETTER
    # outcome -- the direction this catalog entry reports throughout (verified on the
    # worked-example dataset: beta_arm = -0.693, better_outcome_or = exp(0.693) = 2.0).
    beta = res.params[treat]
    ci = res.conf_int().loc[treat]
    better_outcome_or = (float(np.exp(-beta)), float(np.exp(-ci[1])), float(np.exp(-ci[0])))
    return {"beta": float(beta), "better_outcome_or": better_outcome_or, "n": len(d),
            "llf": res.llf, "fit": res}

# ── 2. Descriptive per-cutpoint coefficient diagnostic (NOT a formal Brant test) ─────────
# statsmodels has no packaged Brant test. This only compares the treatment coefficient from
# separate binary logistic fits at each cutpoint; it has no correlation-adjusted covariance
# matrix, degrees of freedom, or p-value, so treat it as descriptive only -- use R's
# brant::brant() below for the real Wald joint test before making a proportionality decision.
def cutpoint_coefficient_diagnostic(cohort: pd.DataFrame, covariates: list[str],
                          treat: str = "arm", outcome: str = "outcome_ord") -> pd.DataFrame:
    d = cohort.copy()
    levels = ["good", "moderate", "poor"]
    rows = []
    for k in range(len(levels) - 1):
        y_bin = (d[outcome].map({lv: i for i, lv in enumerate(levels)}) <= k).astype(int)
        rhs = sm.add_constant(d[[treat] + covariates])
        fit = sm.Logit(y_bin, rhs).fit(disp=0)
        rows.append({"cutpoint": f"<= {levels[k]}", "coef": fit.params[treat],
                     "or": float(np.exp(fit.params[treat]))})
    return pd.DataFrame(rows)   # widely divergent 'or' values across rows is descriptive
                                 # evidence of possible PO violation -- confirm with a formal test

# ── 3. Multinomial baseline-category logit for an UNORDERED treatment-choice outcome ─────
# channel_cohort: person_id, drug_choice ('Drug_A'/'Drug_B'/'Drug_C'), covariates...
def fit_treatment_choice(channel_cohort: pd.DataFrame, covariates: list[str],
                         outcome: str = "drug_choice", reference: str = "Drug_A") -> dict:
    d = channel_cohort.copy()
    cats = [reference] + [c for c in d[outcome].unique() if c != reference]
    y = d[outcome].astype("category").cat.reorder_categories(cats).cat.codes  # ref = 0
    X = sm.add_constant(d[covariates])
    mod = sm.MNLogit(y, X)
    res = mod.fit(disp=0)
    rrr = np.exp(res.params)   # relative risk ratios vs the reference category
    # Generalized propensity score = each person's fitted probability of EVERY category, plus
    # the probability of the category actually chosen (needed for IPTW across 3+ arms).
    gps = pd.DataFrame(res.predict(X).to_numpy(), columns=cats, index=X.index)
    observed_treatment_gps = gps.to_numpy()[np.arange(len(d)), y.to_numpy()]
    return {"rrr": rrr, "gps": gps, "observed_treatment_gps": observed_treatment_gps,
            "n": len(d), "reference": reference, "fit": res}

library(MASS)
library(brant)
library(VGAM)
library(nnet)

# ── 1. Proportional-odds model + formal Brant test ────────────────────────────────────────
# cohort$outcome_ord must be an ordered factor: good < moderate < poor
cohort$outcome_ord <- factor(cohort$outcome_ord, levels = c("good", "moderate", "poor"),
                              ordered = TRUE)

fit_po <- polr(outcome_ord ~ arm + age + comorbidity_score, data = cohort, Hess = TRUE)

# polr() does not print p-values by default; compute them from the t-values.
ctable <- coef(summary(fit_po))
p_values <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2

# polr() fits logit[P(Y<=k)] = alpha_k - beta'X (same convention as statsmodels OrderedModel),
# so exp(coef) is the OR for a WORSE (higher-category) outcome; exp(-coef) is the OR for a
# BETTER outcome, the direction this catalog entry reports. Verified on the worked-example
# dataset: coef(fit_po)["arm"] = -0.693, so exp(-coef) = 2.0.
ci <- confint(fit_po)
worse_outcome_or_ci  <- exp(cbind(OR = coef(fit_po), ci))
better_outcome_or_ci <- exp(cbind(OR = -coef(fit_po),
                                  "2.5 %" = -ci[, "97.5 %"], "97.5 %" = -ci[, "2.5 %"]))

brant_result <- brant(fit_po)   # global + per-covariate proportionality test (the real Wald test)
# A significant row (commonly 'arm') is one diagnostic of a proportionality departure for that
# covariate, not a pass/fail certification; a non-significant result does not prove PO holds,
# and the exposure ('arm') should be assessed on its own regardless of the global test result.

# ── 2. Partial proportional odds when the Brant test flags `arm` ─────────────────────────
fit_partial <- vglm(outcome_ord ~ arm + age + comorbidity_score,
                    family = cumulative(parallel = FALSE ~ arm, link = "logitlink"),
                    data = cohort)
summary(fit_partial)   # arm now has a separate coefficient at each cutpoint

# ── 3. Multinomial logistic for an UNORDERED treatment-choice (channeling) outcome ───────
# channel_cohort$drug_choice: unordered factor with 3+ levels, reference set via relevel()
channel_cohort$drug_choice <- relevel(factor(channel_cohort$drug_choice), ref = "Drug_A")

fit_mn <- multinom(drug_choice ~ age + sex + comorbidity_score, data = channel_cohort,
                   trace = FALSE)
rrr <- exp(coef(fit_mn))                                     # relative risk ratios vs Drug_A
z   <- summary(fit_mn)$coefficients / summary(fit_mn)$standard.errors
p_mn <- (1 - pnorm(abs(z))) * 2

# The fitted probabilities from fit_mn are the generalized propensity scores for IPTW/matching
# across all 3+ treatment arms (Feng et al. 2012); extract with predict(fit_mn, type = "probs").

/* Proportional-odds (cumulative logit) model for the ORDERED outcome (good < moderate < poor).
   order=internal keeps the formatted order good/moderate/poor as coded; SAS auto-prints the
   'Score Test for the Proportional Odds Assumption' -- large chi-square / small p indicates a
   proportionality departure, not a pass/fail proof that the assumption holds either way.
   IMPORTANT -- cross-software sign check: SAS PROC LOGISTIC's cumulative-logit sign convention
   for beta is not guaranteed to match R's MASS::polr or Python's statsmodels OrderedModel (see
   the sign-convention note in the model equation above). Before reporting an odds ratio from
   this PROC LOGISTIC output, confirm its direction against a manual cross-tab of the raw outcome
   counts by arm (as in the worked example) rather than assuming the same beta sign as the R/
   Python code in this entry. */
proc logistic data=work.cohort;
  class arm (ref='Drug_B') / param=ref;
  model outcome_ord (order=internal) = arm age comorbidity_score / link=clogit clodds=pl;
run;

/* If the global score test (or a per-covariate check via separate binary models) flags `arm`,
   partial/generalized ordered models are not directly available in PROC LOGISTIC; fall back to
   fitting separate cumulative binary logistic models at each cutpoint, or use PROC NLMIXED /
   PROC GENMOD with a custom partial-proportional-odds likelihood, or move the analysis to R
   (VGAM::vglm) / Python (statsmodels + manual partial-PO specification). */

/* Generalized (baseline-category / multinomial) logit for an UNORDERED treatment-choice outcome,
   e.g. which of 3 index drugs was dispensed -- the channeling / generalized-propensity-score
   model feeding downstream IPTW across all arms. */
proc logistic data=work.channel_cohort;
  class drug_choice (ref='Drug_A') / param=ref;
  model drug_choice (ref='Drug_A') = age sex comorbidity_score / link=glogit;
  output out=work.gps_scores predprobs=individual;   /* generalized propensity score per arm */
run;

Citations

1. [1]McCullagh P. Regression models for ordinal data. Journal of the Royal Statistical Society: Series B (Methodological). 1980;42(2):109-142.
2. [2]Theil H. A multinomial extension of the linear logit model. International Economic Review. 1969;10(3):251-259.
3. [3]Brant R. Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics. 1990;46(4):1171-1178.
4. [4]Peterson B, Harrell FE Jr. Partial proportional odds models for ordinal response variables. Journal of the Royal Statistical Society: Series C (Applied Statistics). 1990;39(2):205-217.

1. [5]Saver JL. Novel end point analytic techniques and interpreting shifts across the entire range of outcome scales in acute stroke trials. Stroke. 2007;38(11):3055-3062.
2. [6]Optimising Analysis of Stroke Trials (OAST) Collaboration (Bath PM, Gray LJ, Collier T, Pocock S, Carpenter J). Can we improve the statistical analysis of stroke trials? Statistical reanalysis of functional outcomes in stroke trials. Stroke. 2007;38(6):1911-1915.

1. [7]Beigel JD, Tomashek KM, Dodd LE, et al. Remdesivir for the treatment of Covid-19 - final report. New England Journal of Medicine. 2020;383(19):1813-1826.
2. [8]Feng P, Zhou XH, Zou QM, Fan MY, Li XS. Generalized propensity score for estimating the average treatment effect of multiple treatments. Statistics in Medicine. 2012;31(7):681-697.
3. [9]McCaffrey DF, Griffin BA, Almirall D, Slaughter ME, Ramchand R, Burgette LF. A tutorial on propensity score estimation for multiple treatments using generalized boosted models. Statistics in Medicine. 2013;32(19):3388-3414.