読者です 読者をやめる 読者になる 読者になる

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 結果の $coefficientsexp 関数へ渡してオッズ比を算出しました。

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 と同じような結果となりました。

最後に、それぞれの値の分布は下記のようになりました。

f:id:fits:20140302174258p:plain

f:id:fits:20140302174314p:plain