「先生！疾患の重症度（軽症/中等症/重症/最重症）を目的変数にして順序ロジスティック回帰モデルを組もうとし、Brant test（ブラント検定）をやったら $P < 0.05$ と判定されてしまいました！比例オッズ仮定が崩れているから順序ロジスティック回帰は諦めて、多項ロジスティック回帰に変更するしかないでしょうか？」
順序カテゴリカル変数を扱う臨床研究において、このような手詰まり感に直面する研究者は非常に多い。
しかし、Brant test などの統計的仮説検定で $P < 0.05$ が出たからといって、すぐに順序ロジスティック回帰を諦める必要は全くない。
統計的仮説検定はサンプルサイズ（$N$）に極めて強く依存するため、大標本データでは臨床的に無視できるほどわずかな傾きのズレでも機械的に $P < 0.05$（仮定破綻）となり、逆に小標本では重大な傾きのズレを見過ごしてしまう。
正攻法は、Cox回帰の比例ハザード性を検定の $P$ 値ではなくSchoenfeld残差プロットで目視確認するのと全く同じように、「目的変数を各閾値で二値化して独立実行した二値ロジスティック回帰の係数プロット（グラフ法）」で視覚的に平行性を確認することである。
本記事では、検定依存の罠からグラフ法による比例オッズ性の正当な評価アプローチ、Rでの具体的な可視化手順、および論文でそのまま使える英文テンプレートまでを徹底解説する。
1. なぜ「Brant検定の P < 0.05」で順序ロジスティック回帰を諦めてはいけないのか？
臨床現場の罠・手詰まり感
改良Rankinスケール（mRS: 0〜6）、病理学的Grade（1〜4）、疾患重症度（軽症/中等症/重症）などの順序カテゴリ変数を目的変数として解析する際、第一選択となるのが順序ロジスティック回帰（Ordinal Logistic Regression / Proportional Odds Model）である。
教科書通りに比例オッズ性の検証として Brant test や Likelihood Ratio test などを実行し、$P < 0.05$ と判定された際に「多項ロジスティック回帰（Multinomial Logistic Regression）」へ回避するケースが後を絶たない。
しかし、多項ロジスティック回帰に変更すると、カットポイント（境界）ごとに別々のオッズ比（OR）が出力され、変数の数が膨れ上がる。結果として臨床的な解釈が極めて難解になり、「どの因子が患者の重症化を抑えるのか」という本質的な結論がぼやけてしまう。
仮説検定（Brant test等）をお勧めしない絶対的理由（サンプルサイズ依存性）
統計的仮説検定に過度に依存してはいけない理由は、検定結果がサンプルサイズ（$N$）に翻弄されるからである。
- 大標本（Large Sample）の罠:$N$ が数千人〜数万人規模のデータベース研究では、各閾値間での回帰係数のズレが臨床的にほぼ無視できる微小なものであっても、検出力が過剰に高いため機械的に $P < 0.05$（仮定棄却）と判定される。
- 小標本（Small Sample）の罠:サンプルサイズが小さい臨床研究では、あからさまに平行性が崩れていてモデルの適用が危うい状態であっても、検出力不足により $P > 0.05$ となり「仮定充足」と誤認してしまう。
このように、$P$ 値のみでモデルの良否を判定することは非常に危険である。
※補足（関連記事）:
Cox回帰においても比例ハザード性を検定の$P$値ではなくSchoenfeld残差プロットで確認する鉄則については Cox回帰での比例ハザード性破綻とRMST を参照されたい。
2. 比例オッズ性（Parallel Slopes）の直感的意味と「グラフ法」による正攻法
比例オッズ性（Parallel Slopes Assumption）とは何か？
順序ロジスティック回帰（比例オッズモデル）は、目的変数の境界（Cutpoint）をどこに設定しても、説明変数（共変量）が及ぼす効果（対数オッズの傾き $\beta$）が一定であるという前提条件に基づいている。
目的変数を $Y \in \{1, 2, \dots, K\}$ としたとき、閾値 $k$（$k = 1, 2, \dots, K-1$）における累積確率の対数オッズ（Logit）は次式で表される。
※記号の意味: ここで使われている「$\in$」は「〜に属する・〜に含まれる」を意味する数学記号である。つまり$Y \in \{1, 2, \dots, K\}$とは、「目的変数$Y$が$1, 2, \dots, K$という離散的なカテゴリ値のいずれかをとる」ことを表している。
$$\ln \left( \frac{P(Y \ge k)}{P(Y < k)} \right) = \alpha_k + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p$$
ここで注目すべきは、切片 $\alpha_k$ は閾値 $k$ ごとに変化するが、回帰係数 $\beta_1, \beta_2, \dots, \beta_p$ には下付き文字 $k$ がついていない点である。
つまり、「カテゴリ1 vs 2以上」「カテゴリ1,2 vs 3以上」といったどの境界で分割しても、共変量 $X$ の及ぼすオッズ比 $\exp(\beta)$（傾き）は共通であると規定している。これが「平行傾き仮定（Parallel Slopes Assumption）」と呼ばれる所以である。
臨床的に理解する可視化アプローチ（グラフ法）の仕組み
比例オッズ性が臨床的に許容できる範囲で保持されているかを判定するベストアプローチは、閾値別の二値ロジスティック回帰係数プロット（グラフ法）である。
手順は以下の通りシンプルである。
- 順序目的変数 $Y$ を各境界（Cutpoint）で分割し、「0 / 1」の二値変数を $K-1$ 個作成する（例: $Y \ge 2$, $Y \ge 3$, $Y \ge 4$）。
- それぞれの二値変数に対して、独立して通常の「二値ロジスティック回帰モデル（
glm）」を実行する。 - 各閾値で得られた説明変数の回帰係数 $\beta$（対数オッズ）またはオッズ比（OR）とその95%信頼区間を並べてプロットする。
- 各変数の係数がカットポイント間で概ね平行・近接しているか（重なり合っているか）を目視評価する。
点推定値（$\beta$）の並びが概ね平坦であり、95%信頼区間が互いに大きく重なり合っていれば、臨床的に比例オッズ仮定は十分満たされていると判断できる。
3. 比例オッズ性が大きく崩れていた場合の臨床的対処法
グラフ法による視覚的確認で、特定の変数の傾きがあからさまに異なっていた（カットポイントによって効果が全く逆転・急変動している）場合、以下の2つの対処法を検討する。
対処法1：部分比例オッズモデル（Partial Proportional Odds Model）の検討
比例オッズ仮定が保たれている変数には共通の係数 $\beta$ を当てはめ、仮定が大きく崩れている特定の変数のみカットポイントごとの異なる係数 $\beta_k$ を許容する柔軟なモデルである。
Rの VGAM パッケージ（vglm 関数）や rms パッケージ（blrm / clm 等）を用いることで構築可能であり、全変数を分離する多項ロジスティック回帰よりも臨床的解釈性を高く保つことができる。
対処法2：臨床的に本質的な二値化への統合
例えば mRS 0〜6（脳卒中の機能予後評価）を解析する際、臨床試験の主要評価項目として標準的に使われる「臨床的良好（mRS 0-2）vs 不良（mRS 3-6）」のように、医学的妥当性・先行研究の標準基準に基づいた単一の閾値で二値ロジスティック回帰に統合するアプローチである。
4. 【実践Rコード】閾値別二値ロジスティック回帰の可視化と R 実装（Step by Step）
標準パッケージ MASS および tidyverse / ggplot2 を用いて、疑似臨床データを用いた順序ロジスティック回帰の構築、閾値別二値ロジスティック回帰の独立実行、係数プロットの可視化までを Step by Step で解説する。
Step 1: 疑似臨床データの作成と確認
患者300人分、4段階の疾患重症度（grade: 1/2/3/4）を目的変数とし、治療群（treatment）、年齢（age）、BMI（bmi）を説明変数とする疑似データを生成する。
library(tidyverse)
library(MASS)
# 乱数シードの設定
set.seed(123)
n <- 300
# 疑似臨床データの作成
df_ordinal <- tibble(
id = 1:n,
treatment = factor(sample(c("Control", "New"), n, replace = TRUE), levels = c("Control", "New")),
age = round(rnorm(n, mean = 65, sd = 10)),
bmi = round(rnorm(n, mean = 24, sd = 3), 1)
) %>%
mutate(
# 順序目的変数の生成（対数オッズに基づいてカテゴリ化）
latent_score = 0.8 * (treatment == "New") + 0.03 * (age - 65) - 0.05 * (bmi - 24) + rlogis(n),
grade = case_when(
latent_score < 0.0 ~ 1,
latent_score < 1.2 ~ 2,
latent_score < 2.5 ~ 3,
TRUE ~ 4
) %>% factor(ordered = TRUE)
)
# データ構造の確認
table(df_ordinal$grade)
Step 2: 順序ロジスティック回帰の構築とP値の算出（MASS::polr）
MASS::polr() 関数を用いて、標準的な順序ロジスティック回帰モデルを構築し、オッズ比（OR）および95%信頼区間を取得する。
# 順序ロジスティック回帰（Proportional Odds Model）の実行
fit_polr <- polr(grade ~ treatment + age + bmi, data = df_ordinal, Hess = TRUE)
# モデル概要の確認
summary(fit_polr)
# 係数・標準誤差・t値および95%信頼区間の取得・算出
ord_coef <- summary(fit_polr)$coefficients
ord_ci <- confint(fit_polr)
# t値（表記はtだが中身はz値＝Wald統計量）から標準正規分布に基づきP値を算出
p_values <- pnorm(abs(ord_coef[, "t value"]), lower.tail = FALSE)
# オッズ比（OR）・95％信頼区間・P値の統合テーブル作成
df_polr_res <- tibble(
term = rownames(ord_ci),
estimate = ord_coef[rownames(ord_ci), "Value"],
std_error = ord_coef[rownames(ord_ci), "Std. Error"],
t_value = ord_coef[rownames(ord_ci), "t value"],
p_value = p_values[rownames(ord_ci)],
or = exp(estimate),
or_lower = exp(ord_ci[, 1]),
or_upper = exp(ord_ci[, 2])
)
print(df_polr_res)
Step 3: 各閾値での二値ロジスティック回帰の独立実行と係数抽出
目的変数 grade（1〜4）の各閾値（grade >= 2, grade >= 3, grade >= 4）で「0 / 1」に二値化し、独立して glm(family = binomial) を実行する。
# カットポイントリスト（2以上、3以上、4）
cutpoints <- c(2, 3, 4)
# 閾値別二値ロジスティック回帰の一括実行と係数抽出
df_binary_res <- map_dfr(cutpoints, function(k) {
# 閾値 k で二値化データを動的作成
df_temp <- df_ordinal %>%
mutate(y_binary = if_else(as.numeric(grade) >= k, 1, 0))
# 二値ロジスティック回帰
fit_bin <- glm(y_binary ~ treatment + age + bmi, data = df_temp, family = binomial)
# 係数と信頼区間の抽出
ci_bin <- confint(fit_bin)
coef_bin <- summary(fit_bin)$coefficients
tibble(
cutpoint = paste0("Y >= ", k),
term = names(coef(fit_bin))[-1], # 切片を除外
estimate = coef_bin[-1, "Estimate"],
std_error = coef_bin[-1, "Std. Error"],
conf_low = ci_bin[-1, 1],
conf_high = ci_bin[-1, 2]
)
})
print(df_binary_res)
Step 4: 係数プロット（ggplot2）による平行性の視覚的評価
各閾値ごとの回帰係数（対数オッズ $\beta$）と95%信頼区間をプロットし、直線・平坦さの確認（視覚的評価）を行う。
# 係数プロット（平行性の視覚的評価）の描画
ggplot(df_binary_res, aes(x = cutpoint, y = estimate, color = term, group = term)) +
geom_point(position = position_dodge(width = 0.3), size = 3) +
geom_errorbar(aes(ymin = conf_low, ymax = conf_high), width = 0.15, position = position_dodge(width = 0.3)) +
geom_line(position = position_dodge(width = 0.3), linetype = "dashed") +
facet_wrap(~ term, scales = "free_y") +
theme_minimal() +
labs(
title = "Parallel Slopes Assessment via Binary Logistic Regressions",
subtitle = "Log-odds coefficients across cumulative cutpoints (Y >= k)",
x = "Cutpoint Threshold",
y = "Regression Coefficient (Log-Odds)",
color = "Variable"
) +
theme(legend.position = "none")
プロット上で、各変数の回帰係数（点）が閾値間で大きく拡散・逆転することなく横ばい（平行）に推移し、95%信頼区間の帯が十分に重なり合っていれば、比例オッズ仮定は妥当であると自信を持って判定できる。
5. 論文でそのまま使える英文テンプレート（Methods & Results）
論文の Methods および Results セクションに記述する標準的な英文表現テンプレートを示す。
Methods セクション英文テンプレート
“An ordinal logistic regression (proportional odds) model was fitted to evaluate the association between baseline covariates and disease severity grades. The proportional odds assumption was assessed visually by estimating separate binary logistic regression models at each cumulative threshold (Grade >= 2, Grade >= 3, and Grade >= 4) and inspecting the graphical parallelism of the estimated regression coefficients, rather than relying solely on sample-size-dependent formal hypothesis tests.”
Results セクション英文テンプレート
“Separate binary logistic regressions across all cumulative thresholds demonstrated consistent effect estimates across cutpoints, confirming the adequacy of the proportional odds assumption for the primary model. In the ordinal logistic regression model, treatment with the new agent was significantly associated with a higher likelihood of achieving lower severity grades (Common Odds Ratio [OR], 2.02; 95% CI, 1.33 to 3.08; P < 0.001).”
6. まとめ＆一括コピペ用Rスクリプト
本記事のポイントおさらい
- Brant検定の $P < 0.05$ で慌ててモデルを破棄しない: 統計的検定はサンプルサイズに依存し、大標本では微小なズレで有意となり、小標本では検出力不足になる。
- 正攻法は「閾値別二値ロジスティック回帰の係数プロット」: 目的変数を分割して各閾値で
glmを実行し、係数が概ね平行・近接しているかを目視確認する。 - 大幅な破綻時は部分比例オッズモデルや臨床的二値化へ: 視覚的評価で大きく傾きが異なっている場合は、特定変数のみ別係数を許容する部分比例オッズモデルや単一の臨床的閾値への統合を検討する。
一括実行用Rスクリプト
以下のコードをコピー＆ペーストすることで、疑似データ作成から順序ロジスティック回帰、閾値別二値ロジスティック回帰、係数プロットの描画までを一気通貫で実行できる。
# ==============================================================================
# 順序ロジスティック回帰：比例オッズ性の視覚的確認（グラフ法）一括実行スクリプト
# パッケージ: MASS, tidyverse
# ==============================================================================
if (!requireNamespace("MASS", quietly = TRUE)) install.packages("MASS")
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
library(tidyverse)
library(MASS)
# 1. 疑似臨床データの作成
set.seed(123)
n <- 300
df_ordinal <- tibble(
id = 1:n,
treatment = factor(sample(c("Control", "New"), n, replace = TRUE), levels = c("Control", "New")),
age = round(rnorm(n, mean = 65, sd = 10)),
bmi = round(rnorm(n, mean = 24, sd = 3), 1)
) %>%
mutate(
latent_score = 0.8 * (treatment == "New") + 0.03 * (age - 65) - 0.05 * (bmi - 24) + rlogis(n),
grade = case_when(
latent_score < 0.0 ~ 1,
latent_score < 1.2 ~ 2,
latent_score < 2.5 ~ 3,
TRUE ~ 4
) %>% factor(ordered = TRUE)
)
# 2. 順序ロジスティック回帰の構築
cat("\n--- 順序ロジスティック回帰モデル構築 ---\n")
fit_polr <- polr(grade ~ treatment + age + bmi, data = df_ordinal, Hess = TRUE)
summary(fit_polr)
# 係数・標準誤差・t値および95%信頼区間の取得・算出
ord_coef <- summary(fit_polr)$coefficients
ord_ci <- confint(fit_polr)
# t値（表記はtだが中身はz値＝Wald統計量）から標準正規分布に基づきP値を算出
p_values <- pnorm(abs(ord_coef[, "t value"]), lower.tail = FALSE)
# オッズ比（OR）・95％信頼区間・P値の統合テーブル作成
df_polr_res <- tibble(
term = rownames(ord_ci),
estimate = ord_coef[rownames(ord_ci), "Value"],
std_error = ord_coef[rownames(ord_ci), "Std. Error"],
t_value = ord_coef[rownames(ord_ci), "t value"],
p_value = p_values[rownames(ord_ci)],
or = exp(estimate),
or_lower = exp(ord_ci[, 1]),
or_upper = exp(ord_ci[, 2])
)
print(df_polr_res)
# 3. 閾値別二値ロジスティック回帰の実行と係数抽出
cutpoints <- c(2, 3, 4)
df_binary_res <- map_dfr(cutpoints, function(k) {
df_temp <- df_ordinal %>%
mutate(y_binary = if_else(as.numeric(grade) >= k, 1, 0))
fit_bin <- glm(y_binary ~ treatment + age + bmi, data = df_temp, family = binomial)
ci_bin <- confint(fit_bin)
coef_bin <- summary(fit_bin)$coefficients
tibble(
cutpoint = paste0("Y >= ", k),
term = names(coef(fit_bin))[-1],
estimate = coef_bin[-1, "Estimate"],
std_error = coef_bin[-1, "Std. Error"],
conf_low = ci_bin[-1, 1],
conf_high = ci_bin[-1, 2]
)
})
cat("\n--- 閾値別二値ロジスティック回帰 係数一覧 ---\n")
print(df_binary_res)
# 4. 係数プロット（平行性の可視化）の作成
p_parallel <- ggplot(df_binary_res, aes(x = cutpoint, y = estimate, color = term, group = term)) +
geom_point(position = position_dodge(width = 0.3), size = 3) +
geom_errorbar(aes(ymin = conf_low, ymax = conf_high), width = 0.15, position = position_dodge(width = 0.3)) +
geom_line(position = position_dodge(width = 0.3), linetype = "dashed") +
facet_wrap(~ term, scales = "free_y") +
theme_minimal() +
labs(
title = "Parallel Slopes Assessment via Binary Logistic Regressions",
subtitle = "Log-odds coefficients across cumulative cutpoints (Y >= k)",
x = "Cutpoint Threshold",
y = "Regression Coefficient (Log-Odds)",
color = "Variable"
) +
theme(legend.position = "none")
# プロットの表示
print(p_parallel)
コメント