quantifyinghealth.com
注：動作しないコードを書き換え、不足しているものに関して注という形で補完している。
順序ロジスティック回帰は、1つまたは複数の予測変数（数値またはカテゴリ）と順序結果の間の関係をモデルする回帰分析の一種だ。順序結果は、次のような論理的順序を持つ2つ以上のカテゴリを持つ変数である：
- 癌のステージ：0、1、2、3、4
- 所得レベル：低、中、高
このチュートリアルでは、GitHubのリンクからダウンロードできる(https://github.com/datasciencedojo/datasets/blob/master/titanic.csv)。
注：ここではtitanicパッケージのtitanic_trainデータを使用する。
titanicデータセットで順序ロジスティック回帰を使い、予測変数Pclass（乗客クラス）を予測する：Pclassは、3つの値（1: first class, 2: second class, 3: third class）のうちの1つを取ることができる： 予測変数：年齢（Age）、性別（Sex）、生存（Survived）。
library(titanic)
titanic.raw <-titanic_train
## サブセットを作成
titanic <- subset(titanic.raw, select = c(Pclass, Age, Sex, Survived))
## Pclassにレベルをつける
titanic$Pclass <- factor(titanic$Pclass, levels = c(3, 2, 1),
labels = c("1st", "2nd", "3rd"), ordered = TRUE)
titanic$Sex <- factor(titanic$Sex)
titanic$Survived <- factor(titanic$Survived,
labels = c("Died", "Survived"))
str(titanic)
'data.frame': 891 obs. of 4 variables:
$ Pclass : Ord.factor w/ 3 levels "1st"<"2nd"<"3rd": 1 3 1 3 1 1 3 1 1 2 ...
$ Age : num 22 38 26 35 35 NA 54 2 27 14 ...
$ Sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
$ Survived: Factor w/ 2 levels "Died","Survived": 1 2 2 2 1 1 1 1 2 2 ...
2. 順序ロジスティック回帰モデルのフィッティング
我々は、順序ロジスティック回帰の最も一般的なタイプである比例オッズ・ロジスティック回帰を使用する。
このモデルは、ordinal パッケージのclm()関数を用いて、下記のように実行できる：
library(ordinal)
model <- clm(Pclass ~ Age + Sex + Survived, data = titanic, na.action = na.omit)
summary(model)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
Age 0.063208 0.005857 10.792 <2e-16 ***
Sexmale 0.315159 0.199921 1.576 0.115
SurvivedSurvived 1.969870 0.204830 9.617 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Threshold coefficients:
Estimate Std. Error z value
1st|2nd 2.8120 0.2838 9.91
2nd|3rd 4.2024 0.3091 13.60
注
na.actionはリストワイズを意味するna.omit、na.excludeは（予測値、標準誤差、区間限界において）残差値NAと表示される。
各予測変数と結果の関係の解釈
予測変数 X1 に関連づけられた順序ロジスティック回帰係数 β1 の一般的な解釈は，次式である．
予測変数X1が1単位増加すると，結果Yがより高いカテゴリーにあることの対数オッズがβ1変化する．
予測変数の係数は対数オッズ比なので、オッズ比を得るためには指数化しなければならない。
Rでは、係数とp値を抽出するために、モデル・フィットの関数tidy()を、引数exponentiate = TRUEとconf.int = TRUE（95%信頼区間を表示する）で呼び出すことができる。そして、表を小さくするために、出力から標準誤差とzスコアを取り除く：
library(tidymodels)
tidy(model, exponentiate = TRUE, conf.int = TRUE) |>
select(-c(std.error, statistic))
term estimate p.value conf.low conf.high coef.type
<chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 1st|2nd 16.6 3.77e-23 NA NA intercept
2 2nd|3rd 66.8 4.17e-42 NA NA intercept
3 Age 1.07 3.77e-27 1.05 1.08 location
4 Sexmale 1.37 1.15e- 1 0.930 2.04 location
5 SurvivedSurvived 7.17 6.77e-22 4.83 10.8 location
我々は、予測変数の係数（最後の3行）のみを見て、切片は見ない。以下は、その解釈方法 である：
年齢の係数の解釈：
年齢に関する指数係数（オッズ比）は1.07である。年齢が高いことは、より高い客室であることと正の関係があり、解釈は次のようになる：
年齢が高い乗客は、より高いクラスのチケットを持つ傾向がある（p < 0.001）。具体的には、1歳上の乗客は、1.07倍高い確率で、1クラス高い客室のチケットを持つ。
Sexの係数の解釈：
男性のSexに関する指数係数（オッズ比）は統計的に有意ではない（p = 0.115）。
生存率の係数の解釈
Survivedの指数係数（オッズ比）7.17である。
生存した乗客は、より高いクラスのチケットを持っている傾向がある（p < 0.001）。具体的には、生存した乗客は、生存しなかった乗客よりも7.17倍の高い確率で1クラス上の航空券を持っていた。
モデル適合のチェック
rcompanionパッケージの関数nagelkerke()は、擬似R2乗を計算し、尤度比検定を実行する：
library(rcompanion)
nagelkerke(model)
$Models
Model: "clm, Pclass ~ Age + Sex + Survived, titanic"
Null: "clm, Pclass ~ 1, titanic"
$Pseudo.R.squared.for.model.vs.null
Pseudo.R.squared
McFadden 0.155862
Cox and Snell (ML) 0.277186
Nagelkerke (Cragg and Uhler) 0.316640
$Likelihood.ratio.test
Df.diff LogLik.diff Chisq p.value
-3 -115.88 231.77 5.738e-50
$Number.of.observations
Model: 714
Null: 714
$Messages
[1] "Note: For models fit with REML, these statistics are based on refitting with ML"
$Warnings
[1] "None"
特に重要なのは、
NagelkerkeのR2乗: ロジスティック回帰モデルの適合度を測定する0から1の間の数値である。
尤度比検定: フル・モデル（すべての予測変数を含むモデル）がヌル・モデル（変数を含まないモデル）よりもよくデータに適合するかどうかを検定する。我々の事例では、LogLik.diffは-115.88で、p < 0.001であり、これは予測変数を追加することが予測変数を含まない帰無モデルよりもよいことを意味する。
注
nagelkerkeの疑似R二乗を用いた尤度比検定は行わない方が良いと思われる。一つはnagelkerkeの疑似R二乗が通常のR二乗とは異なるということ、もう一つはR二乗はがモデルのフィッティングを示す指標としては不適切であるためである。
参考：
https://ides.hatenablog.com/entry/2022/05/06/165738
https://ides.hatenablog.com/entry/2022/06/19/030458
また、lipsitz検定を実行して適合度をチェックすることもできる：
library(generalhoslem)
lipsitz.test(model)
Lipsitz goodness of fit test for ordinal response models
data: formula: Pclass ~ Age + Sex + Survived
LR statistic = 5.2796, df = 9, p-value = 0.8093
帰無仮説は、指定されたモデルがデータによく適合する。5%水準以下の小さなP値は、モデルが不完全であることを示唆する。つまり、我々のモデルはデータに適合している。
注
リプシッツ検定は、順序対応ロジスティック回帰モデルの適合度検定である。これは、順序応答スコアに基づいて、観察されたデータを等しい大きさのg群に分割することを含む。このスコアは、各結果水準での各被験者の予測確率に等間隔の整数の重みを掛けて合計することにより計算される。ユーザーは、デフォルトでは10であるgに整数値を代入して、グループの数を指定できる。
データのこの分割が与えられると、各グループについて、被験者が領域gにあればI = 1、なければI = 0となるようなダミー変数Iが導かれる。そして、モデルをこれらのダミー変数で再フィットする。Lipsitz ら (1996) は、尤度比検定、Wald 検定、スコア検定が使用できることを示唆している; lipsitz.test は、自由度 g-1 の尤度比検定を使用するだけである。
モデルを実行する前に、結果変数を因子に変換しなければならないことに注意してほしい。モデル関数内で as.factor() を使用すると、 lipsitz.test がモデルを再適合するために update() 関数を使用する方法のため、エラーが発生する。
Lipsitz 検定は、順序Hosmer-Lemeshow 検定（logitgof）とPulkstenis-Robinson 検定（pulkrob.chisq と pulkrob.deviance）と一緒に実行することが推奨さ れる（Fagerland and Hosmer, 2016）。
https://www.rdocumentation.org/packages/generalhoslem/versions/1.3.4/topics/lipsitz.test
注：Pulkstenis-Robinson 検定
帰無仮説は、指定されたモデルがデータによく適合する。5%水準以下の小さなP値は、モデルが不完全であることを示唆する。
c("Sex","Survived")のところはカテゴリカル変数を指定する必要がある。
pulkrob.chisq(model, c("Sex","Survived"))
Pulkstenis-Robinson chi-squared test
data: formula: Pclass ~ Age + Sex + Survived
X-squared = 28.637, df = 11, p-value = 0.002583
pulkrob.deviance(model, c("Sex","Survived"))
Pulkstenis-Robinson deviance test
data: formula: Pclass ~ Age + Sex + Survived
Deviance-squared = 28.448, df = 11, p-value = 0.002763
注：Hosmer-Lemeshow検定
帰無仮説は、指定されたモデルがデータによく適合する。5%水準以下の小さなP値は、モデルが不完全であることを示唆する。
## 欠損地を削除してtitanic.na.omitに格納
titanic.na.omit <- na.omit(titanic)
## 従属変数のみのデータをpredprobに格納
predprob <- data.frame(Age = titanic.na.omit$Age, Sex = titanic.na.omit$Sex, Survived = titanic.na.omit$Survived)
## フィットさせたオブジェクトmodelから必要なものをfvに格納
fv <- predict(model, newdata = predprob, type = "prob")$fit
## Hosmer-Lemeshow検定を実行，従属変数、フィットさせたオブジェクト、従属変数のカテゴリー数、順序変数かどうかを指定する。
logitgof(titanic.na.omit$Pclass, fv, g = 3, ord = TRUE)
Hosmer and Lemeshow test (ordinal model)
data: titanic.na.omit$Pclass, fv
X-squared = 3.5784, df = 3, p-value = 0.3107
順序ロジスティック回帰モデルの精度
次に、3つのステップでモデルの精度を計算する：
ステップ#1: フィット値を取得し、predsに保存する：
注：ここではbroomパッケージのaugment()を用いる。
preds <- augment(model, type = "class", na.action="na.omit")
preds
# A tibble: 714 × 6
.rownames Pclass Age Sex Survived .fitted
<chr> <ord> <dbl> <fct> <fct> <fct>
1 1 1st 22 male Died 1st
2 2 3rd 38 female Survived 3rd
3 3 1st 26 female Survived 3rd
4 4 3rd 35 female Survived 3rd
5 5 1st 35 male Died 1st
6 7 3rd 54 male Died 3rd
7 8 1st 2 male Died 1st
8 9 1st 27 female Survived 3rd
9 10 2nd 14 female Survived 1st
10 11 1st 4 female Survived 1st
ステップ2：混同行列を見る：
注：yardstickのconf_mat()を用いる。
conf_mat(preds, truth = Pclass, estimate = .fitted)
Truth
Prediction 1st 2nd 3rd
1st 306 114 73
2nd 0 0 0
3rd 49 59 113
このモデルは、2等席のチケットを持っている人を分類していないことに注意してほしい。
ここで停止して、混同行列を用いて手動で精度を計算することもできますが、yardstickのaccuracy()関数を用いるとより速く計算できるでしょう。
ステップ3: モデルの精度を計算する：
yardstick::accuracy(preds, truth = Pclass, estimate = .fitted)
.metric .estimator .estimate
<chr> <chr> <dbl>
1 accuracy multiclass 0.587
つまり、年齢、性別、生存率に基づいて、モデルは乗客のクラスを58.7％の確率で正しく推測していることになる。
比例オッズの仮定のチェック
上記の順序ロジスティック回帰モデルの特定のタイプ（比例オッズ・ロジスティック回帰）は、結果での各予測変数の効果が結果のすべてのレベルについて同じであると仮定している。たとえば、2等席のチケット vs 1等席のチケットを持つオッズでの年齢の1歳増加の効果は、3等席のチケット vs 2等席のチケットを持つオッズでの年齢の1歳増加の効果と同じといったものだ。
ここで帰無仮説は、比例オッズの仮定が成り立つというものです。帰無仮説は，比例オッズの仮定が成り立つというものである．オムニバス検定でp < 0.05と変数の少なくとも1つが重なると，仮定が破られたとみなされる[source： McNulty K. Handbook of Regression Modeling in People Analytics： With Examples in R and Python. 第1版。Chapman and Hall/CRC; 2021.].
library(gofcat)
brant.test(model)
Brant Test:
chi-sq df pr(>chi)
Omnibus 3.699 3 0.30
Age 2.027 1 0.15
Sexmale 1.554 1 0.21
SurvivedSurvived 0.679 1 0.41
H0: Proportional odds assumption holds
ここでは、すべてのp値が0.05を超えているので、比例オッズの仮定が成り立つ。
