R でロジスティック回帰を行う方法

タイタニック号は、1912年4月14日の夜、氷山に激突し、北大西洋の底に1,500名以上の命と一緒に沈んだ。

乗客乗員の生存・死亡のデータを用いて、ロジスティック回帰分析を実行してみる。

タイタニック号の生存・死亡データはどこにあるのか？

R はインストールするとデフォルトで使用可能な状態になっている。

datasets パッケージの文字通り Titanic という名前のデータだ。

Titanic の内容を確認する

str()で内容を確認してみる。

str(Titanic)

通常のdata frameとは異なる構造になっている。

> str(Titanic)
'table' num [1:4, 1:2, 1:2, 1:2] 0 0 35 0 0 0 17 0 118 154 ...
- attr(*, "dimnames")=List of 4
..$ Class   : chr [1:4] "1st" "2nd" "3rd" "Crew"
..$ Sex     : chr [1:2] "Male" "Female"
..$ Age     : chr [1:2] "Child" "Adult"
..$ Survived: chr [1:2] "No" "Yes"

データ自体を表示させてみる。

Titanic

Child, AdultとSurvived Yes, Noの4つに分かれたClassとSexの分割表が表示された。

このような形式で格納されている。

> Titanic
, , Age = Child, Survived = No
Sex
Class  Male Female
1st     0      0
2nd     0      0
3rd    35     17
Crew    0      0
, , Age = Adult, Survived = No
Sex
Class  Male Female
1st   118      4
2nd   154     13
3rd   387     89
Crew  670      3
, , Age = Child, Survived = Yes
Sex
Class  Male Female
1st     5      1
2nd    11     13
3rd    13     14
Crew    0      0
, , Age = Adult, Survived = Yes
Sex
Class  Male Female
1st    57    140
2nd    14     80
3rd    75     76
Crew  192     20

通常のデータフレーム状にするには、as.data.frame()を使う。

as.data.frame(Titanic)

4つの変数の組み合わせごとの人数というデータの形に変形された。

> as.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
21   1st Female Child      Yes    1
22   2nd Female Child      Yes   13
23   3rd Female Child      Yes   14
24  Crew Female Child      Yes    0
25   1st   Male Adult      Yes   57
26   2nd   Male Adult      Yes   14
27   3rd   Male Adult      Yes   75
28  Crew   Male Adult      Yes  192
29   1st Female Adult      Yes  140
30   2nd Female Adult      Yes   80
31   3rd Female Adult      Yes   76
32  Crew Female Adult      Yes   20

データフレームに転換したdatで計算していく。

dat <- as.data.frame(Titanic)

性別の相対生存オッズ比をロジスティック回帰モデルで計算する

単変量のロジスティック回帰

性別と生存・死亡を掛け合わせただけのシンプルなオッズ比をロジスティック回帰で計算してみる。

glm()を使い、family=binomialとする。「~」の左に生存・死亡のデータ Survived 、右に性別データ Sexを配置する。

今回のデータは人数の集計データなので、weights=Freqを付け加える。

生データでないときにはこれが必要だ。

logi1 <- glm(Survived~Sex,weights=Freq,data=dat,family=binomial)
summary(logi1)

結果は、logオッズは、男性を0とすると、女性は2.3172と計算された。

統計学的に有意である。

> logi1 <- glm(Survived~Sex,weights=Freq,data=dat,family=binomial)
> summary(logi1)
Call:
glm(formula = Survived ~ Sex, family = binomial, data = dat,
weights = Freq)
Deviance Residuals:
Min       1Q   Median       3Q      Max
-17.869   -3.455    0.000    5.969   24.405
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.3128     0.0588  -22.32   <2e-16 ***
SexFemale     2.3172     0.1196   19.38   <2e-16 ***
---
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 23  degrees of freedom
Residual deviance: 2335.0  on 22  degrees of freedom
AIC: 2339
Number of Fisher Scoring iterations: 5

対数を真数に変換したほうがわかりやすい。

自前の関数 rr.95ci() を使う。

rr.95ci <- function(x){
res <- exp(cbind(coef(x),confint(x)))
colnames(res) <- c("RR","LL(2.5%)","UL(97.5%)")
round(res,2)
}

rr.95ci は Relative Risk and 95% Confidence Intervalの略だ。

rr.95ci(logi1)

女性のオッズ比は、男性を1とすると、10.15と計算された。

それだけ優先されて助けられたという結果だ。

> rr.95ci(logi1)
Waiting for profiling to be done...
RR LL(2.5%) UL(97.5%)
(Intercept)  0.27     0.24      0.30
SexFemale   10.15     8.05     12.86

Fisher Exact Test で確認してみる

2×2の分割表のデータに対して、Fisher Exact Test を実施してみる。

ロジスティック回帰と同様の結果になるか確認してみた。

apply()で、性別と生存・死亡のデータの分割表を、他の要因は合計するという処置をする。

そのテーブルを行列としてFisher Exact Testを実施する。

apply(Titanic,c(2,4),sum)
fisher.test(apply(Titanic,c(2,4),sum))

テーブルの数値を見れば、一目瞭然で、男性に比べ女性のほうが生き残っている。

Fisher Exact Testの結果、性別と生存・死亡は、統計学的有意に独立で、オッズ比は10.13であった。

これは先ほどのロジスティック回帰の結果と同様だ。

> apply(Titanic,c(2,4),sum)
Survived
Sex        No Yes
Male   1364 367
Female  126 344
> fisher.test(apply(Titanic,c(2,4),sum))
Fisher's Exact Test for Count Data
data:  apply(Titanic, c(2, 4), sum)
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  7.97665 12.92916
sample estimates:
odds ratio 
   10.1319 

客室・性別・年齢を加味した多重ロジスティック回帰モデルで分析する

客室 Class、性別 Sex、年齢 Ageを同時に投入した多変量モデルで分析してみる。

glm()に、family=binomialを加えるのは同様。「~」の左は同様にSurvived、右はClass, Sex, Ageを「+」でつないで投入。

logi2 <- glm(Survived ~ Class+Sex+Age, weights=Freq, family=binomial, data=dat)
summary(logi2)
rr.95ci(logi2)

結果として全部の説明変数が統計学的有意であった。

一等客室の乗客に比べ、二等、三等、乗務員は誰もが生き残れなかった。

女性のほうが生き残ったのは単変量の結果と同じ。

子供に比べ大人のほうが生き残れなかった。

子供を優先したということだろう。

> summary(logi2)
Call:
glm(formula = Survived ~ Class + Sex + Age, family = binomial,
data = dat, weights = Freq)
Deviance Residuals:
Min       1Q   Median       3Q      Max
-18.505   -4.247    0.000    4.747   23.915
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 23  degrees of freedom
Residual deviance: 2210.1  on 18  degrees of freedom
AIC: 2222.1
Number of Fisher Scoring iterations: 5
> rr.95ci(logi2)
Waiting for profiling to be done...
RR LL(2.5%) UL(97.5%)
(Intercept)  1.98     1.16      3.39
Class2nd     0.36     0.25      0.53
Class3rd     0.17     0.12      0.24
ClassCrew    0.42     0.31      0.58
SexFemale   11.25     8.57     14.87
AgeAdult     0.35     0.21      0.56

Fisher Exact Testで単変量解析もやってみる

大人・子供と生存・死亡の関係をFisher Exact Testでも実行してみる。

apply(Titanic,c(3,4),sum)
fisher.test(apply(Titanic,c(3,4),sum))

子供のほうが生存した割合が高そうだ。

> apply(Titanic,c(3,4),sum)
Survived
Age       No Yes
Child   52  57
Adult 1438 654

大人のほうが統計学的有意に生き残れなく、オッズ比は0.415であった。

> fisher.test(apply(Titanic,c(3,4),sum))
Fisher's Exact Test for Count Data
data:  apply(Titanic, c(3, 4), sum)
p-value = 1.234e-05
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.2760923 0.6228448
sample estimates:
odds ratio 
 0.4150843 

客室クラスと生存・死亡の単変量解析はどうやったいいか？

Class変数は、客室クラスと乗務員を合わせて、4つに分かれている。

apply(Titanic,c(1,4),sum)

4×2の分割表の時はどうしたらよいか？

> apply(Titanic,c(1,4),sum)
Survived
Class   No Yes
1st  122 203
2nd  167 118
3rd  528 178
Crew 673 212

epitools パッケージを使う。

Epidemiology (疫学) のツールが詰まったパッケージである。

最初に一回だけインストールする。

install.packages("epitools")

library()で呼び出す。

library(epitools)

oddsratio()を使うと、一行目をオッズ1として、二行目以降のオッズ比を計算してくれる。

oddsratio(apply(Titanic,c(1,4),sum))

$measureの項に結果があり、一等客室の乗客に比べ、二等は0.425、三等は0.203、乗務員は0.190と計算されている。

多変量モデル中の結果とおおむね一致する。

> oddsratio(apply(Titanic,c(1,4),sum))
$data
Survived
Class     No Yes Total
1st    122 203   325
2nd    167 118   285
3rd    528 178   706
Crew   673 212   885
Total 1490 711  2201
$measure
odds ratio with 95% C.I.
Class   estimate     lower     upper
1st  1.0000000        NA        NA
2nd  0.4254346 0.3064855 0.5884387
3rd  0.2030991 0.1528567 0.2686306
Crew 0.1897539 0.1441302 0.2487196
$p.value
two-sided
Class    midp.exact fisher.exact   chi.square
1st            NA           NA           NA
2nd  2.029144e-07 2.777190e-07 2.026267e-07
3rd  0.000000e+00 3.675395e-30 1.141202e-30
Crew 0.000000e+00 1.811851e-34 6.880615e-36
$correction
[1] FALSE
attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"

まとめ

悲劇の大事故を起こしたタイタニック号の生存・死亡データを、ロジスティック回帰モデルで分析してみた。

多変量の結果だけでなく、単変量の結果とも突き合わせながら、実施した。

ロジスティック回帰モデルの分析方法は、glm()を使い、family=binomialを追加する。

今回のように集計値のデータセットの場合、as.data.frame()でデータフレームに変換し、weights=Freq（Freqは頻度変数の名前）を忘れないようにする。