R でロジスティック回帰とオッズ比の算出 - glm, MCMClogit
以前、glm・MCMCmetrop1R 関数でロジスティック回帰を試みましたが、今回はその時に利用を断念した MCMCpack の MCMClogit 関数を使ってロジスティック回帰を行います。
題材は、書籍 「 データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus) 」の p.38 と同様の Titanic データセットを使ったロジスティック回帰とオッズ比の算出です。
以前使ったデータは主に数値でしたが、今回は主に因子（factor）データを扱っている点が異なっています。
今回のソースは http://github.com/fits/try_samples/tree/master/blog/20140302/
はじめに
MCMCpack パッケージを R へインストールしておきます。
install.packages("MCMCpack")
Titanic データセット
Titanic のデータセットは R にデフォルトで用意されており、data.frame(Titanic) すると下記のようになります。
data.frame(Titanic) 結果
> data.frame(Titanic) Class Sex Age Survived Freq 1 1st Male Child No 0 2 2nd Male Child No 0 3 3rd Male Child No 35 4 Crew Male Child No 0 5 1st Female Child No 0 6 2nd Female Child No 0 7 3rd Female Child No 17 8 Crew Female Child No 0 9 1st Male Adult No 118 10 2nd Male Adult No 154 11 3rd Male Adult No 387 12 Crew Male Adult No 670 13 1st Female Adult No 4 14 2nd Female Adult No 13 15 3rd Female Adult No 89 16 Crew Female Adult No 3 17 1st Male Child Yes 5 18 2nd Male Child Yes 11 19 3rd Male Child Yes 13 20 Crew Male Child Yes 0 ・・・
内容は以下の通りです。
| Class（船室等級） | Sex（性別） | Age（年齢層） | Survived（生存可否） | Freq（人数） |
|---|---|---|---|---|
| 1st, 2nd, 3rd, Crew | Male, Female | Child, Adult | No, Yes | 数値 |
今回は、ロジスティック回帰の結果から船室等級・性別・年齢層毎の生存率に対するオッズ比を算出します。
(1) glm を使ったロジスティック回帰とオッズ比
書籍の内容とほとんど同じですが、
まずは glm 関数を使ったロジスティック回帰です。
Survived~. としているので Class・Sex・Age を説明変数としたロジスティック回帰を実施します。
オッズ比は exp(<推定値>) で算出できるので、書籍のような epicalc パッケージは使わず、glm 結果の $coefficients を exp 関数へ渡してオッズ比を算出しました。
logiMcmcglmm.R
d <- data.frame(Titanic) d.data <- data.frame( Class = rep(d$Class, d$Freq), Sex = rep(d$Sex, d$Freq), Age = rep(d$Age, d$Freq), Survived = rep(d$Survived, d$Freq) ) d.res <- glm(Survived~., data = d.data, family = binomial) summary(d.res) # オッズ比 exp(d.res$coefficients)
なお、上記では data.frame(Titanic) の Freq（人数） を展開したデータフレーム（下記）を使ってロジスティック回帰を実施しています。
d.data の内容
Class Sex Age Survived 1 3rd Male Child No 2 3rd Male Child No ・・・ 35 3rd Male Child No 36 3rd Female Child No ・・・
実行結果は下記の通りです。
実行結果
・・・
> d.res <- glm(Survived~., data = d.data, family = binomial)
> summary(d.res)
Call:
glm(formula = Survived ~ ., family = binomial, data = d.data)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0812 -0.7149 -0.6656 0.6858 2.1278
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6853 0.2730 2.510 0.0121 *
Class2nd -1.0181 0.1960 -5.194 2.05e-07 ***
Class3rd -1.7778 0.1716 -10.362 < 2e-16 ***
ClassCrew -0.8577 0.1573 -5.451 5.00e-08 ***
SexFemale 2.4201 0.1404 17.236 < 2e-16 ***
AgeAdult -1.0615 0.2440 -4.350 1.36e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2769.5 on 2200 degrees of freedom
Residual deviance: 2210.1 on 2195 degrees of freedom
AIC: 2222.1
Number of Fisher Scoring iterations: 4
>
> # オッズ比
> exp(d.res$coefficients)
(Intercept) Class2nd Class3rd ClassCrew SexFemale AgeAdult
1.9844057 0.3612825 0.1690159 0.4241466 11.2465380 0.3459219
オッズ比から、女性（Female）は男性（Male）の 11倍（11.2465380）、大人（Adult）は子供（Child）の 3分の1 （0.3459219）の生存率という結果になりました。
(2) MCMClogit を使ったロジスティック回帰とオッズ比
それでは、MCMCpack の MCMClogit 関数を使って同様の処理を実施してみます。
glm 関数の場合とほとんど同じですが、MCMClogit 関数は応答変数（下記の Survived）に因子型（factor）を扱えないようなので、Survived の No と Yes を 0（= No） と 1（= Yes）へ変換しています。
（as.numeric(d$Survived) - 1 の箇所）
また、MCMClogit 関数の結果も glm とは違って値の分布となるので summary した結果の平均値（期待値） Mean を使ってオッズ比を算出しています。
logistic_odds_mcmclogit.R
library(MCMCpack) d <- data.frame(Titanic) d.data <- data.frame( Class = rep(d$Class, d$Freq), Sex = rep(d$Sex, d$Freq), Age = rep(d$Age, d$Freq), Survived = rep(as.numeric(d$Survived) - 1, d$Freq) ) d.res <- MCMClogit(Survived~., data = d.data) d.summ <- summary(d.res) d.summ #オッズ比 exp(d.summ$statistics[, "Mean"])
実行結果は下記の通りです。
実行結果
・・・
> d.res <- MCMClogit(Survived~., data = d.data)
>
> d.summ <- summary(d.res)
>
> d.summ
Iterations = 1001:11000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept) 0.6775 0.2682 0.002682 0.012059
Class2nd -1.0367 0.1961 0.001961 0.008828
Class3rd -1.7935 0.1768 0.001768 0.008054
ClassCrew -0.8607 0.1596 0.001596 0.006971
SexFemale 2.4486 0.1384 0.001384 0.006279
AgeAdult -1.0540 0.2392 0.002392 0.010882
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
(Intercept) 0.1455 0.4992 0.6831 0.8630 1.2028
Class2nd -1.4185 -1.1683 -1.0331 -0.9027 -0.6626
Class3rd -2.1481 -1.9130 -1.7926 -1.6760 -1.4476
ClassCrew -1.1764 -0.9661 -0.8651 -0.7535 -0.5488
SexFemale 2.1543 2.3600 2.4525 2.5442 2.7110
AgeAdult -1.5275 -1.2144 -1.0573 -0.8904 -0.5907
>
> #オッズ比
> exp(d.summ$statistics[, "Mean"])
(Intercept) Class2nd Class3rd ClassCrew SexFemale AgeAdult
1.9690188 0.3546137 0.1663746 0.4228566 11.5715841 0.3485381
glm と同じような結果となりました。
最後に、それぞれの値の分布は下記のようになりました。
