医学研究において、順序カテゴリ（ordinal category）を持つアウトカムは頻繁に遭遇します。疾患の重症度（軽度、中等度、重度）、患者の自覚的健康状態、疼痛スコアなど、カテゴリ間に自然な順序が存在するデータです。通常のロジスティック回帰は2値アウトカムのみを扱うため、このような順序カテゴリには適用できません。順序ロジスティック回帰（ordered logistic regression）は、こうした順序カテゴリデータを適切に解析するための手法です。 ## 順序ロジスティック回帰とは 順序ロジスティック回帰は、従属変数が3つ以上の順序カテゴリを持つ場合に用いる回帰モデルです。カテゴリ間の順序情報を考慮した解析が可能になります。 ### 通常のロジスティック回帰との違い 通常のロジスティック回帰は、アウトカムが2値（例：発症した/しなかった）の場合に使用します。一方、順序ロジスティック回帰は以下の特徴を持ちます。 - アウトカムが3つ以上の順序カテゴリを持つ - カテゴリ間の順序関係を考慮する - 比例オッズ仮定（proportional odds assumption）を用いる ### 順序ロジスティック回帰の考え方 順序ロジスティック回帰では、「あるカテゴリ以下である確率」と「それより上のカテゴリである確率」を比較します。 例えば、血栓の程度が「なし」「軽度」「中等度」の3段階ある場合を考えてみます。このとき、以下の2つの比較を行います。 - なし vs. 軽度／中等度 - なし／軽度 vs. 中等度 順序ロジスティック回帰の特徴は、説明変数（例：年齢、血流量など）の関連が、これらすべての比較で共通であると仮定する点です。これを「比例オッズ仮定」と呼びます。 具体的には、「年齢が1歳高いと血栓が起きやすい」という関連性が、上記のどの比較でも同じ強さで成り立つと仮定します。この仮定により、1つのモデルで複数のカテゴリ間の関係を効率的に推定できます。 生物統計における仮定の考え方については、[[生物統計における仮定とは]]を参照してください。 ## 実践：Rによる解析 ### 必要なパッケージのインストールと読み込み 順序ロジスティック回帰には`MASS`パッケージの`polr()`関数を使用します。`polr`は"proportional odds logistic regression"の略です。また、比例オッズ仮定を検定するために`brant`パッケージも使用します。 ```r # パッケージのインストール（初回のみ） install.packages("MASS") # 順序ロジスティック回帰の実行用 install.packages("brant") # 比例オッズ仮定の検定用 # パッケージの読み込み library(MASS) library(brant) ``` ### データセットの準備 透析終了時の回路における血栓の程度を例として、順序ロジスティック回帰を実行してみます。アウトカムは血栓の程度（1=なし、2=軽度、3=中等度）、説明変数は血流量、抗血小板薬使用の有無、年齢とします。 ```r # データフレームの作成 dat <- data.frame( id = 1:100, thrombosis = c(1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,1,1,1,1,1, 1,1,1,1,1,2,1,2,1,1,1,1,1,1,2,1,1,2,1,2,1,2,1,1,1, 1,1,1,2,2,1,1,1,1,2,1,1,1,1,3,1,1,1,1,1,1,1,1,1,2, 1,1,1,1,1,1,1,1,1,2,1,2,1,1,1,1,1,1,1,1,3,1,1,1,1), blood_flow = c(220,209,170,196,223,226,211,197,176,178,167,209,217,206,155, 197,181,209,160,155,223,206,221,196,208,169,198,214,159,198, 165,193,165,224,190,208,207,190,189,175,183,200,172,219,165, 203,168,213,213,219,189,185,187,179,165,170,218,195,220,160, 212,217,211,204,145,197,160,214,210,158,223,176,206,202,155, 194,208,214,179,200,175,220,165,165,203,222,196,195,195,216, 175,170,186,214,211,159,177,176,199,155), antiplatelet = c(0,1,0,1,0,1,0,1,1,0,0,0,1,1,0,1,1,0,1,1,1,0,0,1,1, 0,1,0,1,1,1,1,0,1,0,0,1,1,0,0,1,1,1,1,1,1,1,0,1,0, 0,1,1,0,1,1,1,1,0,1,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0, 1,0,1,1,0,0,0,1,0,1,1,0,0,0,1,0,0,0,0,0,1,0,0,0,1), age = c(46,59,89,75,51,25,64,56,51,64,81,60,60,79,55,52,62,39,75,86, 56,71,54,62,65,74,75,62,69,55,68,63,68,52,55,69,72,56,52,66, 53,47,67,73,75,67,75,66,88,38,60,66,51,62,72,48,51,62,71,65, 52,67,53,27,68,65,68,66,53,73,52,78,62,59,60,63,59,66,57,67, 76,68,57,56,60,61,59,54,61,49,63,73,76,74,49,83,53,63,47,62) ) # 順序変数としてのfactor型に変換 dat$thrombosis <- factor(dat$thrombosis, levels = c(1, 2, 3), labels = c("なし", "軽度", "中等度"), ordered = TRUE) # 二値変数をfactor型に変換 dat$antiplatelet <- factor(dat$antiplatelet, levels = c(0, 1), labels = c("なし", "あり")) ``` ### データの確認 まず、データセットの構造を確認します。 ```r # データ構造の確認 str(dat) # 出力: # 'data.frame': 100 obs. of 5 variables: # $ id : int 1 2 3 4 5 6 7 8 9 10 ... # $ thrombosis : Ord.factor w/ 3 levels "なし"<"軽度"<..: 1 1 2 1 1 1 1 1 1 1 ... # $ blood_flow : num 220 209 198 196 223 226 211 197 176 178 ... # $ antiplatelet: Factor w/ 2 levels "なし","あり": 1 2 1 2 1 2 1 2 2 1 ... # $ age : num 46 59 89 75 51 25 64 56 51 64 ... ``` 次に、記述統計を確認します。 ```r # 記述統計 summary(dat) # 出力: # id thrombosis blood_flow antiplatelet age # Min. : 1.0 なし :83 Min. :145.0 なし:46 Min. :25.00 # 1st Qu.: 25.8 軽度 :15 1st Qu.:175.0 あり:54 1st Qu.:54.75 # Median : 50.5 中等度: 2 Median :196.0 Median :62.00 # Mean : 50.5 Mean :192.4 Mean :62.09 # 3rd Qu.: 75.2 3rd Qu.:211.0 3rd Qu.:68.25 # Max. :100.0 Max. :226.0 Max. :89.00 ``` 血栓の程度は、なし83例（83%）、軽度15例（15%）、中等度2例（2%）と分布しています。血流量は平均192.4 mL/分、抗血小板薬使用者は54例（54%）、年齢は平均62.1歳です。 アウトカムの分布を詳しく確認します。 ```r # 頻度表 table(dat$thrombosis) # 出力: # なし 軽度 中等度 # 83 15 2 # 割合表 prop.table(table(dat$thrombosis)) # 出力: # なし 軽度 中等度 # 0.83 0.15 0.02 ``` ### 単変量解析 まず、血流量単独で順序ロジスティック回帰を実行します。`MASS`パッケージの`polr()`関数を使用します。 #### 血流量と血栓の程度 ```r # 単変量モデル：血流量 model_blood <- polr(thrombosis ~ blood_flow, data = dat, Hess = TRUE) summary(model_blood) # 出力: # Call: # polr(formula = thrombosis ~ blood_flow, data = dat, Hess = TRUE) # # Coefficients: # Value Std. Error t value # blood_flow -0.09293 0.02188 -4.246 # # Intercepts: # Value Std. Error t value # なし|軽度 -15.2399 3.7900 -4.0211 # 軽度|中等度 -12.3992 3.7026 -3.3488 # # Residual Deviance: 71.65286 # AIC: 77.65286 ``` 係数は-0.0929で、血流量が1 mL/分増加すると、より高い血栓カテゴリに属するオッズの対数が0.0929減少します。統計的に有意です（t = -4.246）。オッズ比に変換してみます。 ```r # オッズ比の計算 exp(coef(model_blood)) # 出力: # blood_flow # 0.911261 # 信頼区間の計算 exp(confint(model_blood)) # 出力: # 2.5 % 97.5 % # 0.8666100 0.9467987 ``` 血流量のオッズ比は0.911（95%CI: 0.867-0.947）で、血流量が1 mL/分増加すると、より高い血栓カテゴリに属するオッズが0.911倍（約9%減少）になります。信頼区間が1を含まないため、統計的に有意です。 出力に表示される`なし|軽度`と`軽度|中等度`はカットポイント（threshold coefficient）と呼ばれ、カテゴリ間の閾値を表します。これらは通常、解釈には使用しません。 ### 多変量順序ロジスティック回帰 次に、すべての説明変数を同時に投入した多変量モデルを構築します。 ```r # 多変量モデル model_multi <- polr(thrombosis ~ blood_flow + antiplatelet + age, data = dat, Hess = TRUE) summary(model_multi) # 出力: # Call: # polr(formula = thrombosis ~ blood_flow + antiplatelet + age, data = dat, Hess = TRUE) # # Coefficients: # Value Std. Error t value # blood_flow -0.08508 0.02201 -3.866 # antiplateletあり 0.85288 0.68605 1.243 # age 0.06964 0.03533 1.971 # # Intercepts: # Value Std. Error t value # なし|軽度 -8.7530 4.6500 -1.8824 # 軽度|中等度 -5.5523 4.6491 -1.1943 # # Residual Deviance: 66.02672 # AIC: 76.02672 ``` オッズ比と信頼区間を計算します。 ```r # オッズ比の計算 exp(coef(model_multi)) # 出力: # blood_flow antiplateletあり age # 0.9184418 2.3464013 1.0721268 # 信頼区間の計算 exp(confint(model_multi)) # 出力: # 2.5 % 97.5 % # blood_flow 0.8734015 0.9547687 # antiplateletあり 0.6323051 9.6603791 # age 1.0023841 1.1533789 ``` オッズ比は、血流量0.918、抗血小板薬2.346、年齢1.072です。信頼区間を見ると、血流量（95%CI: 0.873-0.955）と年齢（95%CI: 1.002-1.153）の信頼区間は1を含まないため統計的に有意です。一方、抗血小板薬（95%CI: 0.632-9.660）の信頼区間は1を含むため統計的に有意ではありません。 ### 結果のまとめ表示 オッズ比と信頼区間をまとめて表示する関数を作成すると便利です。 ```r # 係数表の取得（説明変数のみ） ctable <- coef(summary(model_multi)) # 説明変数のみを抽出（切片を除外） coef_only <- ctable[1:3, ] # 最初の3行が説明変数 # p値の計算（t値から） p_values <- pnorm(abs(coef_only[, "t value"]), lower.tail = FALSE) * 2 # 信頼区間の計算 CI <- exp(confint(model_multi)) # 結果をまとめる result <- data.frame( Variable = rownames(coef_only), OR = round(exp(coef_only[, "Value"]), 2), CI_low = round(CI[, 1], 2), CI_high = round(CI[, 2], 2), p_value = round(p_values, 3) ) print(result) # 出力: # Variable OR CI_low CI_high p_value # blood_flow blood_flow 0.92 0.87 0.95 0.000 # antiplateletあり antiplateletあり 2.35 0.63 9.66 0.214 # age age 1.07 1.00 1.15 0.049 ``` - **血流量**: オッズ比は0.92（95%CI: 0.87-0.95）で、血流量が1 mL/分増加すると、より高い血栓カテゴリに属するオッズが0.92倍（約8%減少）になることを示します（p < 0.001）。血流量が低いほど血栓形成が多くなる傾向があります。 - **抗血小板薬使用**: オッズ比は2.35（95%CI: 0.63-9.66）で、抗血小板薬を使用している場合、より高い血栓カテゴリに属するオッズが2.35倍になることを示しますが、統計的に有意ではありません（p = 0.214）。 - **年齢**: オッズ比は1.07（95%CI: 1.00-1.15）で、年齢が1歳増加すると、より高い血栓カテゴリに属するオッズが1.07倍（約7%増加）になることを示します（p = 0.049）。高齢になるほど血栓形成が多くなる傾向があります。 ## 比例オッズ仮定の検定 順序ロジスティック回帰の重要な仮定は、比例オッズ仮定（proportional odds assumption、平行回帰の仮定とも呼ばれる）です。この仮定は、各説明変数の効果がすべてのカテゴリ間の比較で一定であることを意味します。 例えば、血栓の程度が「なし」「軽度」「中等度」の3カテゴリの場合、以下の2つの比較があります。 - なし vs. 軽度/中等度 - なし/軽度 vs. 中等度 比例オッズ仮定では、血流量の効果（オッズ比）がこれら両方の比較で同じであると仮定します。もしこの仮定が成立しない場合、順序ロジスティック回帰の結果は適切ではなく、より柔軟なモデルが必要になります。 ### 比例オッズ仮定の検定（Brant検定） 先ほど推定した多変量モデル`model_multi`を使って、Brant検定で比例オッズ仮定を検定します。Brant検定は比例オッズ仮定が成立しているかを統計的に検証する標準的な方法です。 ```r # Brant検定の実行 brant(model_multi) # 出力: # ---------------------------------------------------- # Test for X2 df probability # ---------------------------------------------------- # Omnibus 1.46 3 0.69 # blood_flow 1.46 1 0.23 # antiplateletあり 0.00 1 1.00 # age 0.14 1 0.71 # ---------------------------------------------------- # # H0: Parallel Regression Assumption holds ``` ### 結果の解釈 Brant検定の帰無仮説は「比例オッズ仮定（平行回帰の仮定）が成立している」です。 - **Omnibus（全体検定）**: カイ二乗値 = 1.46、p = 0.69 - 全変数をまとめた検定結果 - p > 0.05のため、全体として比例オッズ仮定は成立していると判断 - **個別変数の検定**: - **blood_flow**: p = 0.23（仮定は成立） - **antiplateletあり**: p = 1.00（仮定は成立） - **age**: p = 0.71（仮定は成立） すべてのp値が0.05より大きいため、比例オッズ仮定は成立していると判断できます。 重要な注意点として、**有意なp値（p < 0.05）が出た場合は比例オッズ仮定違反の証拠となります**。その場合は、多項ロジスティック回帰（multinomial logistic regression）の使用を検討する必要があります。 ### 比例オッズ仮定違反への対処 もしBrant検定で仮定違反が検出された場合、以下の対処法があります： 1. **多項ロジスティック回帰を使用**: カテゴリ間の順序を考慮しないモデルを推定 2. **カテゴリの統合**: アウトカムのカテゴリ数を減らして再解析 3. **変数の変換**: 問題となっている説明変数を変換（例：連続変数をカテゴリ化） 多項ロジスティック回帰は順序を考慮しないため、解釈は複雑になりますが、比例オッズ仮定に縛られません。 ## 留意点 ### 比例オッズ仮定 順序ロジスティック回帰の最も重要な仮定です。この仮定が成立しない場合、推定結果にバイアスが生じる可能性があります。Brant検定などで仮定を確認することが推奨されます。 ### サンプルサイズ 順序ロジスティック回帰は、各カテゴリに十分な観測数が必要です。カテゴリ数が多い場合や、一部のカテゴリに観測数が少ない場合、推定が不安定になる可能性があります。 ### 完全分離（Complete Separation） 説明変数がアウトカムを完全に予測する場合、モデルの推定に問題が生じます。このような場合、標準誤差が非常に大きくなり、推定結果の信頼性が低下します。データの分布を確認し、必要に応じて説明変数の調整やカテゴリの統合を検討してください。 ### カテゴリの順序 順序ロジスティック回帰は、カテゴリ間の順序が意味を持つ場合にのみ適用できます。順序に意味がない場合（例：血液型）は、多項ロジスティック回帰（multinomial logistic regression）を使用してください。 ### 変数のスケーリング 説明変数のスケールが大きく異なる場合、「Model is nearly unidentifiable: large eigenvalue ratio」という警告が表示されることがあります。この場合、変数を標準化（平均0、標準偏差1に変換）することで警告を解消できますが、解釈が複雑になります。警告が出ても、推定結果自体は通常問題ありません。 ## 参考コード（完全版） 以下は、データ作成から多変量解析、比例オッズ仮定の検定までを一度に実行できる完全なコードです。 ```r # 透析血栓データの作成と順序ロジスティック回帰の実行 # パッケージの読み込み library(MASS) library(brant) # データフレームの作成 dat <- data.frame( id = 1:100, thrombosis = c(1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,1,1,1,1,1, 1,1,1,1,1,2,1,2,1,1,1,1,1,1,2,1,1,2,1,2,1,2,1,1,1, 1,1,1,2,2,1,1,1,1,2,1,1,1,1,3,1,1,1,1,1,1,1,1,1,2, 1,1,1,1,1,1,1,1,1,2,1,2,1,1,1,1,1,1,1,1,3,1,1,1,1), blood_flow = c(220,209,170,196,223,226,211,197,176,178,167,209,217,206,155, 197,181,209,160,155,223,206,221,196,208,169,198,214,159,198, 165,193,165,224,190,208,207,190,189,175,183,200,172,219,165, 203,168,213,213,219,189,185,187,179,165,170,218,195,220,160, 212,217,211,204,145,197,160,214,210,158,223,176,206,202,155, 194,208,214,179,200,175,220,165,165,203,222,196,195,195,216, 175,170,186,214,211,159,177,176,199,155), antiplatelet = c(0,1,0,1,0,1,0,1,1,0,0,0,1,1,0,1,1,0,1,1,1,0,0,1,1, 0,1,0,1,1,1,1,0,1,0,0,1,1,0,0,1,1,1,1,1,1,1,0,1,0, 0,1,1,0,1,1,1,1,0,1,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0, 1,0,1,1,0,0,0,1,0,1,1,0,0,0,1,0,0,0,0,0,1,0,0,0,1), age = c(46,59,89,75,51,25,64,56,51,64,81,60,60,79,55,52,62,39,75,86, 56,71,54,62,65,74,75,62,69,55,68,63,68,52,55,69,72,56,52,66, 53,47,67,73,75,67,75,66,88,38,60,66,51,62,72,48,51,62,71,65, 52,67,53,27,68,65,68,66,53,73,52,78,62,59,60,63,59,66,57,67, 76,68,57,56,60,61,59,54,61,49,63,73,76,74,49,83,53,63,47,62) ) # 順序変数としてのfactor型に変換 dat$thrombosis <- factor(dat$thrombosis, levels = c(1, 2, 3), labels = c("なし", "軽度", "中等度"), ordered = TRUE) # 二値変数をfactor型に変換 dat$antiplatelet <- factor(dat$antiplatelet, levels = c(0, 1), labels = c("なし", "あり")) # データの確認 str(dat) summary(dat) # 頻度表 table(dat$thrombosis) prop.table(table(dat$thrombosis)) # クロス集計 table(dat$thrombosis, dat$antiplatelet) # 単変量解析：血流量 model_blood <- polr(thrombosis ~ blood_flow, data = dat, Hess = TRUE) summary(model_blood) exp(coef(model_blood)) # オッズ比 exp(confint(model_blood)) # 信頼区間 # 多変量順序ロジスティック回帰 model_multi <- polr(thrombosis ~ blood_flow + antiplatelet + age, data = dat, Hess = TRUE) summary(model_multi) # オッズ比の計算 exp(coef(model_multi)) # 信頼区間の計算 exp(confint(model_multi)) # 結果のまとめ表示 ctable <- coef(summary(model_multi)) coef_only <- ctable[1:3, ] # 説明変数のみ抽出 p_values <- pnorm(abs(coef_only[, "t value"]), lower.tail = FALSE) * 2 CI <- exp(confint(model_multi)) result <- data.frame( Variable = rownames(coef_only), OR = round(exp(coef_only[, "Value"]), 2), CI_low = round(CI[, 1], 2), CI_high = round(CI[, 2], 2), p_value = round(p_values, 3) ) print(result) # 比例オッズ仮定の検定（Brant検定） brant(model_multi) ```