タイタニック号は、いまから100年以上前の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 で確認してみる
2x2の分割表のデータに対して、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)
4x2の分割表の時はどうしたらよいか？
> 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"
まとめ
100年以上も前に悲劇の大事故を起こしたタイタニック号の生存・死亡データを、ロジスティック回帰モデルで分析してみた。
多変量の結果だけでなく、単変量の結果とも突き合わせながら、実施した。
ロジスティック回帰モデルの分析方法は、glm()を使い、family=binomialを追加する。
今回のように集計値のデータセットの場合、as.data.frame()でデータフレームに変換し、weights=Freq（Freqは頻度変数の名前。任意）を忘れないようにする。
