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

タイタニックに乗ったときの生還確率

はじめに

f:id:silver801:20161208132149j:plain
タイタニックの沈没、世界でいちばん有名な海難事故といっていいんじゃないでしょうか。都市伝説も多いですよね。実は船がすりかえられてたんじゃないかとか。機会がなくて映画はみたことがないんですけど、母がよく「あれは恋愛ものじゃなくてホラーだ」と言っているのでホラー映画なんだと思います。みたことないけど。

で、そのタイタニックの乗員のデータ、つまり

「名前、客室の等級、性別、年齢、一緒に乗っていた兄弟と配偶者の人数、一緒に乗っていた親と子供の人数、支払った運賃、どの港から乗ったか、生還したか、etc」


WebHome < Main < Vanderbilt Biostatistics Wiki

で手に入ることがわかったので、いろいろ解析してみようと思います。ヴァンタービルト大学ありがとう。*1
で、データをいじってみて「タイタニックに乗ったときの生還確率をもとめる計算式」を作り出すことを目標にします。
計算式だけがみたい人は結論のところまで飛ばそう!

解析

統計解析にはR言語を使います。

子供かどうかを付けくわえる

Rstudioで生データを見てみるとこんな感じ。
f:id:silver801:20161208105414p:plain
大人か子供かのデータがないので、12歳以下を子供として「子供かどうか?」の項目isChildを元々のデータに追加する。
Rの使いかたがよくわかってないのでここはEXCELをつかう。いきなり雰囲気が所帯じみてきた。

=IF((セル) = "","",IF((セル)<=12,1,0))

f:id:silver801:20161208111750p:plain
年齢が空欄の人がいるので処理がガチャガチャしてるけど、とりあえずできました。
あとの方で必要になったので、同じようにして、男性か女性かを1,0で表すisMaleという項目も追加しました。行き当たりばったり

それぞれのデータをみていく

生存率と男女比
setwd('/Users/***/Documents/R/Titanic')
data <- read.csv("titanic3.csv")

survivor <- sum(data$survived)
dead <- nrow(data) - survivor
bars <- c(survivor,dead)
deadper <- round((dead/(survivor+dead))*100,digits = 1)
deadper <- paste(as.character(deadper),"% Died")

male <- sum(data$sex == "male")
female <- sum(data$sex == "female") 

sex = c(male,female)

barplot(bars,names.arg = c("survived","dead"),col = c("blue","red"),main = deadper)
barplot(sex,names.arg = c("male","female"),col = c("blue","red"),main = "sex",ylim = c(0,1000))

ディレクトリ名で本名が見えちゃってたのでそこだけ***に変えました。
setwdをしないとcsvが読み込めなかったり、sumに空白を読ませてしまうと上手く動作しなかったりなど問題はあったがとりあえずできた。
f:id:silver801:20161208112547p:plain
90%ぐらい亡くなったような気持ちでいたが死亡率は61.8%。意外と低いとか思ってしまった。もちろんそんなことはないんだけど。
男女比はこんな感じ。
f:id:silver801:20161208114054p:plain

年齢と運賃

年齢と運賃は大事そうなデータですね。高い料金を払っている客や子供は優先的に避難できそうな気がします。
運賃について補足すると、タイタニックの客室は一等客室から三等客室にわかれています。しかし料金は三種類ではなくて、各等級のなかでさらに細かく段階があります。オプションが細かいのかもしれません。ここらへんは詳しい人にお聞きしたい。
では、まず年齢のヒストグラムから。

age <- data$age
h <- hist(age)

f:id:silver801:20161208115454p:plain
横軸は5歳刻み。豪華客船というと壮年の金持ちが乗っているイメージですが、意外と若者が多いんですね。
次は運賃。

fare <-data$fare
f <-hist(fare,breaks = "FD")

f:id:silver801:20161208124323p:plain
運賃の単位はポンド。正確さは保証できませんがインターネットの情報によると、当時の3ポンドが現在の350ドルにあたるらしいです。
5−10ポンドを払って乗船している人がおおいので、10万円ぐらい払っている客がいちばん多かったということでしょうか。飛行機がないことを考えるとまあまあ良心的なのでは?
あと、見にくくて恐縮ですがヒストグラムの右端に500ポンドぐらい払ってるとんでもない方々がいますね。この三人は有名な富豪らしいです。
500ポンドというと現在の貨幣価値でいうと600万円ぐらいになるのかな?時代が違うとはいえすごい。ちなみに500ポンド払った人たちは全員助かってます。リッチマンズワールド。

重回帰分析

じゃあ、さっそく重回帰分析*2をつかって生還確率の計算式を求めていきます。
変数としては
isMale:男かどうか?
pclass:客室の等級は?
age:年齢は?
をつかってみます。いろいろ試したらこれがいちばんマシな結果だった。
コードは2行で済みます。イルマティックな言語だなあ。

dat.lm <- lm(dat$survived ~ isMale+pclass+age,data = dat)
summary(dat.lm)

summaryを表示するとこんな感じ。

Residuals:
     Min       1Q   Median       3Q      Max 
-1.08309 -0.25801 -0.08043  0.21930  0.99758 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.2762817  0.0540773   23.60  < 2e-16 ***
age         -0.0052002  0.0009285   -5.60 2.73e-08 ***
pclass      -0.1827877  0.0160407  -11.39  < 2e-16 ***
isMale      -0.4914873  0.0255501  -19.24  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3913 on 1042 degrees of freedom
  (263 observations deleted due to missingness)
Multiple R-squared:  0.3686,	Adjusted R-squared:  0.3668 
F-statistic: 202.8 on 3 and 1042 DF,  p-value: < 2.2e-16

何がなんだかわからないと思いますが僕もよくわかりません。
Coefficientsというのが係数で、これは大事。age、pclass、isMaleというパラメータに対してこの係数をかけるわけですね。たとえばageの係数はだいたい-0.001なので、年齢が1歳あがるごとに生存率が0.001さがると仮定していることになります。
あとR-squaredというのも大事。これが1に近いほど回帰式が確かだということになります。
0.3686というのはぜんぜん良くないですね。業界用語でいうと”弱い相関”と言うらしい。
まあ、統計的有意ではあるのでよしとしましょう。よし。

結論

生存率 = -0.0052002×年齢 + (-0.1827877)×客室の等級 + (-0.4914873)×男かどうか + 1.2762817
「客室の等級」は一等客室から三等客室までなので1~3のどれかの数字、「男かどうか」には男性なら1、女性なら0が入ります。
いちおう計算フォームを用意してみたので、「年齢→客室の等級→男かどうか」の順に入力すると結果がわかります。
1に近いほど生存率が高いことになります。


×(-0.0052002)+×(-0.1827877)+×(-0.4914873)+1.2762817



ガバガバの分析だったのでもうちょっといろいろ勉強したいなあと思った。
お付き合いいただきありがとうございました!

*1:統計界では定番の練習問題らしくて、みんなやってるので参考URLはたくさんあるんですが、特にこの記事を参考にしました [Python]Pandasでタイタニック号の乗客データを解析する - Qiita

*2:生存確率=a*年齢+b*客室の等級+....みたいな感じの計算式を考えて、その係数a,b,...をもとめるという手法。正直よくわかってない。