井出草平の研究ノート

順序ロジスティック回帰とBrant検定[R]

前回と同じく順序ロジスティック回帰モデルの話。今回はBrant検定を利用したパターン。

ides.hatenablog.com

使用するのはMASSパッケージのpolr関数。

www.rdocumentation.org

データ

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となり、前の仮定が満たされていないことを示す証拠がないことを示している。したがって、順序回帰という方法はサンプルデータセットに適していると判断できる。