順序型邏輯迴歸與 shift plot

分級的結果，被二分之後剩下什麼

臨床研究有一整類結果既不是連續的、也不是有無，而是有序的分級：

• mRS（modified Rankin Scale，改良 Rankin 量表）0–6，中風試驗的標準結果
• NYHA 心衰竭功能分級 I–IV
• CTCAE 不良事件嚴重度 1–5 級
• 疼痛評分 0–10、咳嗽嚴重度 0–2、內視鏡分級、病理分期

這些分級的共同性質是有順序但沒有間距：mRS 從 2 進到 3 與從 4 進到 5，都叫「差一級」， 但那兩件事對病人的意義完全不同。所以既不能當成連續變項去算平均（間距沒有意義）， 也不該直接二分成「好 vs 壞」（會把大量訊息丟掉）。

最常見的處理是二分：「mRS 0–2 算 good outcome」。這麼做的代價是 從 mRS 5 進步到 mRS 3 的病人，在分析裡與完全沒有進步的病人被算成同一件事。 中風試驗因此改用 shift analysis（位移分析），問的不是「多少人達標」， 而是「整個分布有沒有往好的方向移動」。它的統計工具就是這一頁的順序型邏輯迴歸 （ordinal logistic regression），而它報出來的那個數字，就是讀者最常撞到卻最少被解釋的 common odds ratio（共同勝算比）。

這一頁的例子

沿用群集與重複測量資料那一頁的甘草含漱液試驗 （medicaldata::licorice_gargle），但換一個結果：拔管當下的咳嗽嚴重度 extubation_cough，分成 0 = 沒有咳嗽、1 = 輕度、2 = 中重度 三級。 233 位病人有紀錄（2 位缺失）。

這張圖就是 shift analysis 的視覺形式。 它把「分布整體往哪邊移動」直接畫出來， 而不是先選一個切點再報一個比例。注意兩件事：無咳嗽的比例上升了 14 個百分點， 而中重度的比例從 13% 掉到 3%—— 後面這個變化的相對幅度大得多，但如果你只二分成「有沒有咳嗽」，它會完全消失。

模型長什麼樣

順序型邏輯迴歸最常用的形式是比例勝算模型（proportional odds model）。 它不是對每一個等級各配一個模型，而是對每一個切點同時建模。三個等級有兩個切點： 「≥ 1 vs 0」與「≥ 2 vs ≤ 1」。MASS::polr 是從每個切點的下側寫這個模型的， 下面印出來的結果就是這個形式：

關鍵在於：每個切點有自己的截距 $\zeta_k$，但共用同一個 $\beta$。 截距吸收「這個等級本來有多常見」，斜率吸收「治療把整個分布推了多遠」。

本頁模型的兩個截距是 0.417（切在 0|1） 與 2.104（切在 1|2）， 治療的係數是 -0.723（標準誤 0.283，t = -2.55）。

取指數並翻成「對照組相對於甘草組」的方向，得到 common OR 2.060（95% CI 1.18–3.59）。

動手跑一次

library(medicaldata)
library(MASS)

lg <- licorice_gargle
d  <- subset(data.frame(y = lg$extubation_cough, treat = lg$treat), !is.na(y))

table(d$treat, d$y)                     # 先看格子有多空

fit <- polr(factor(y, ordered = TRUE) ~ factor(treat), data = d, Hess = TRUE)
summary(fit)
exp(-coef(fit))                         # common OR：對照組 vs 甘草組
exp(-rev(confint.default(fit)))         # 95% CI

# proportional odds：每個切點各配一個二元 logistic
summary(glm(I(y >= 1) ~ factor(treat), data = d, family = binomial))
summary(glm(I(y >= 2) ~ factor(treat), data = d, family = binomial))

# 因為只有一個二元共變項，2x3 表可以直接做飽和模型的概似比檢定
tb <- table(d$treat, d$y)
ll_sat <- sum(ifelse(tb > 0, tb * log(tb / rowSums(tb)), 0))
lrt <- 2 * (ll_sat - as.numeric(logLik(fit)))
c(statistic = lrt, p = pchisq(lrt, df = 1, lower.tail = FALSE))

# 格子太空會怎樣：換成 pacu30min_cough
ds <- subset(data.frame(y = lg$pacu30min_cough, treat = lg$treat), !is.na(y))
table(ds$treat, ds$y)
summary(polr(factor(y, ordered = TRUE) ~ factor(treat), data = ds, Hess = TRUE))
summary(glm(I(y >= 2) ~ factor(treat), data = ds, family = binomial))

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.miscmodels.ordinal_model import OrderedModel

lg = sm.datasets.get_rdataset("licorice_gargle", "medicaldata").data
d = lg[["extubation_cough", "treat"]].dropna().rename(
    columns={"extubation_cough": "y"})

print(np.asarray(sm.stats.Table.from_data(d[["treat", "y"]]).table_orig))

mod = OrderedModel(d["y"].astype(int), d[["treat"]], distr="logit").fit(method="bfgs")
print(mod.summary())
print(np.exp(-mod.params["treat"]))     # common OR：對照組 vs 甘草組

# 每個切點各一個二元 logistic
for k in (1, 2):
    d[f"ge{k}"] = (d["y"] >= k).astype(int)
    print(smf.logit(f"ge{k} ~ treat", data=d).fit(disp=0).summary())

假設怎麼檢查

最直覺、也最容易向臨床讀者解釋的檢查方式是：把每一個切點各配一個二元 logistic 迴歸， 把得到的 log OR 並排看。如果比例勝算假設在資料裡站得住，這些數字應該接近。

兩個切點的 log OR 相差 0.780， 換算成 OR 是 2.18 倍。看起來差滿多的——但這個「看起來」值多少錢， 要看不確定性有多大。

因為本頁只有一個二元共變項，這張 2×3 表可以直接做精確的概似比檢定： 飽和模型有四個自由參數、比例勝算模型有三個，所以檢定有 1 個自由度。 結果是 χ² = 2.411，p = 0.121。

格子太空的時候

同一份資料換一個時點：PACU 30 分鐘的咳嗽 pacu30min_cough。這個變項的分布是

全體只有 4 位病人達到中重度， 而且這 4 位全部落在對照組，甘草組是 0 位。這一格是空的。

等級多不等於資訊多

同一份資料的喉嚨痛評分 pacu30min_throatPain 有 7 個等級（0 到 6），比咳嗽的三級多得多。但它的次數分布是：

73% 的人是 0 分，尾端有 2 個等級 的人數少於五位，最少的只有 1 位。等級數是量表設計的性質， 有效資訊量是資料的性質，兩者無關。 這種分布配順序型模型時，尾端的切點基本上是靠假設在支撐， 所以本頁不對它配模型——實務上會先合併尾端等級，或改用不需要比例勝算假設的做法。

違反比例勝算假設時怎麼辦

檢查發現各切點的 OR 明顯不同時，選項由簡到繁：

常見誤用

相關頁面

• 群集與重複測量資料——同一份 licorice_gargle 資料，改看四個時點的喉嚨痛評分
• 邏輯迴歸——二元結果的基礎，以及完全分離長什麼樣子
• 多重比較——為什麼切點要事先指定
• 卡方與 Fisher 精確檢定——列聯表的無母數檢定，不使用順序資訊

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-13-ordinal.R

回報這一頁的問題