前回と同じく順序ロジスティック回帰モデルの話。今回はBrant検定を利用したパターン。
使用するのはMASSパッケージのpolr関数。
データ
library(MASS) data(housing) dat<- housing head(dat)
Sat Infl Type Cont Freq 1 Low Low Tower Low 21 2 Medium Low Tower Low 21 3 High Low Tower Low 28 4 Low Medium Tower Low 34 5 Medium Medium Tower Low 22 6 High Medium Tower Low 36
https://rdrr.io/cran/MASS/man/housing.html
Sat: 住宅所有者の現在の住宅環境に対する満足度（高、中、低の順に順序が設定されている）。
Infl: 所有者が不動産の管理に与える影響の度合い（「高い」「中程度」「低い」の順に回答。
Type: 賃貸住宅のタイプ（タワー型、アトリウム型、アパート型、テラス型）。
Cont: 居住者が他の居住者と接する機会があるか（低い、高い）。
Freq: 各クラスの住人の数。
カテゴリカル変数の順序を設定する
もし順序が設定されていなかった場合には順序を設定しておく。今回はしなくてよい。
dat$Sat<- factor(dat$Infl, ordered = T,
levels = c("low", "High"))
順序の確認はstr(dat)で確認できる。
'data.frame': 72 obs. of 5 variables: $ Sat : Ord.factor w/ 2 levels "low"<"High": NA NA NA NA NA NA 2 2 2 NA ... $ Infl: Factor w/ 3 levels "Low","Medium",..: 1 1 1 2 2 2 3 3 3 1 ... $ Type: Factor w/ 4 levels "Tower","Apartment",..: 1 1 1 1 1 1 1 1 1 2 ... $ Cont: Factor w/ 2 levels "Low","High": 1 1 1 1 1 1 1 1 1 1 ... $ Freq: int 21 21 28 34 22 36 10 11 36 61 ...
データの概観
library(summarytools) dfSummary(dat)
-----------------------------------------------------------------------------------------------------------------
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
---- ------------------- -------------------------- -------------------- ------------------- ---------- ---------
1 Sat 1. Low 24 (33.3%) IIIIII 72 0
[ordered, factor] 2. Medium 24 (33.3%) IIIIII (100.0%) (0.0%)
3. High 24 (33.3%) IIIIII
2 Infl 1. High 24 (33.3%) IIIIII 72 0
[ordered, factor] 2. Medium 24 (33.3%) IIIIII (100.0%) (0.0%)
3. Low 24 (33.3%) IIIIII
3 Type 1. Tower 18 (25.0%) IIIII 72 0
[factor] 2. Apartment 18 (25.0%) IIIII (100.0%) (0.0%)
3. Atrium 18 (25.0%) IIIII
4. Terrace 18 (25.0%) IIIII
4 Cont 1. High 36 (50.0%) IIIIIIIIII 72 0
[ordered, factor] 2. Low 36 (50.0%) IIIIIIIIII (100.0%) (0.0%)
5 Freq Mean (sd) : 23.3 (17.7) 39 distinct values : : 72 0
[integer] min < med < max: : : : (100.0%) (0.0%)
3 < 19.5 < 86 : : :
IQR (CV) : 21.8 (0.8) : : : : .
: : : : : . . . .
-----------------------------------------------------------------------------------------------------------------
モデルの実行
Hess = TRUEで最適化過程で得られる数値近似値が得られる。
weights = Freqはスタック形式のデータであるため。1ケース1行のデータセットの場合、ここは不要。
順序変数の末尾が参照カテゴリとなる。
res01 <- polr(formula = Sat ~ Infl + Type + Cont,
data = housing, weights = Freq,
Hess = TRUE, method = "logistic")
summary(res01)
Call:
polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq,
Hess = TRUE, method = "logistic")
Coefficients:
Value Std. Error t value
InflMedium 0.5664 0.10465 5.412
InflHigh 1.2888 0.12716 10.136
TypeApartment -0.5724 0.11924 -4.800
TypeAtrium -0.3662 0.15517 -2.360
TypeTerrace -1.0910 0.15149 -7.202
ContHigh 0.3603 0.09554 3.771
Intercepts:
Value Std. Error t value
Low|Medium -0.4961 0.1248 -3.9739
Medium|High 0.6907 0.1255 5.5049
Residual Deviance: 3479.149
AIC: 3495.149
オッズ比の表示
odds01 <- exp(coef(res01)) odds01
InflMedium InflHigh TypeApartment TypeAtrium TypeTerrace ContHigh
1.7619017 3.6284990 0.5641979 0.6933734 0.3358754 1.4337368
P値の表示
ctable <- coef(summary(res01)) p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2 ctable <- cbind(ctable, "p-value" = p) ctable
Value Std. Error t value p-value InflMedium 0.5663937 0.1046528 5.412123 6.228192e-08 InflHigh 1.2888191 0.1271561 10.135720 3.835394e-24 TypeApartment -0.5723501 0.1192380 -4.800064 1.586146e-06 TypeAtrium -0.3661866 0.1551733 -2.359855 1.828208e-02 TypeTerrace -1.0910149 0.1514860 -7.202083 5.929964e-13 ContHigh 0.3602841 0.0955358 3.771195 1.624675e-04 Low|Medium -0.4961353 0.1248472 -3.973939 7.069368e-05 Medium|High 0.6907083 0.1254719 5.504882 3.694147e-08
平行性の検定
平行性の仮定または比例オッズの仮定は、順序付きカテゴリ変数の順序ロジスティック回帰モデルを適用するために必要になる。そうしないのであれば、多項式モデルを使用する必要がある。Brant検定を用いて検定することができ，Brantパッケージのbrant関数で利用できる。
library(brant) brant(res01)
-------------------------------------------- Test for X2 df probability -------------------------------------------- Omnibus 0 6 1 InflMedium 0 1 1 InflHigh 0 1 1 TypeApartment 0 1 1 TypeAtrium 0 1 1 TypeTerrace 0 1 1 ContHigh 0 1 1 -------------------------------------------- H0: Parallel Regression Assumption holds
P値が0.05より大きいと、平行性の仮定が成り立つという帰無仮説が棄却されないか、または、結果変数の異なるカットポイントで係数が異ならないということになる。出力によるとP値が1となり、前の仮定が満たされていないことを示す証拠がないことを示している。したがって、順序回帰という方法はサンプルデータセットに適していると判断できる。
