kaggleでデータ分析の基礎 #3【titanic, tidyr, dplyr,ggplot】

[#1はこちら](https://qiita.com/H_Y_J/items/71ddb72e0968e4f1d900)
[前回#2はこちら](https://qiita.com/H_Y_J/items/de59bfc34f3a47134e07)

#2.3 欠損値

いろいろ変数について調べました。学んだことから欠損値についていくつか考えてみます。

乗船場所のNAは補間してもいいのでは？
部屋番号Cabinには欠損値が多すぎます。
年齢は下記のアプローチをしてみます。
運賃のひとつかけている部分も値を入れてみましょう。

```R
combine %>% filter(is.na(Embarked))
## # A tibble: 2 x 12
##   PassengerId Survived Pclass Name    Sex     Age SibSp Parch Ticket  Fare
##         <dbl> <fct>    <fct>  <chr>   <fct> <dbl> <dbl> <dbl> <chr>  <dbl>
## 1        62.0 1        1      Icard,… fema…  38.0     0     0 113572  80.0
## 2       830   1        1      Stone,… fema…  62.0     0     0 113572  80.0
## # ... with 2 more variables: Cabin <chr>, Embarked <fct>
```
乗船位置が欠損している二人の情報を見てみましょう。
どちらも一等客室の女性で(38,62歳)家族は連れていません。

```R
combine %>%
  filter(Embarked != "Q" & Pclass == 1 & Sex == "female") %>%
  group_by(Embarked, Pclass, Sex, Pclass, Parch, SibSp) %>%
  summarise(count = n())
## # A tibble: 16 x 6
## # Groups:   Embarked, Pclass, Sex, Parch [?]
##    Embarked Pclass Sex    Parch SibSp count
##    <fct>    <fct>  <fct>  <dbl> <dbl> <int>
##  1 C        1      female  0     0       30
##  2 C        1      female  0     1.00    20
##  3 C        1      female  1.00  0       10
##  4 C        1      female  1.00  1.00     6
##  5 C        1      female  2.00  0        2
##  6 C        1      female  2.00  2.00     2
##  7 C        1      female  3.00  1.00     1
##  8 S        1      female  0     0       20
##  9 S        1      female  0     1.00    20
## 10 S        1      female  0     2.00     3
## 11 S        1      female  1.00  0        7
## 12 S        1      female  1.00  1.00     6
## 13 S        1      female  2.00  0        4
## 14 S        1      female  2.00  1.00     5
## 15 S        1      female  2.00  3.00     3
## 16 S        1      female  4.00  1.00     1
```
家族を持たず、一等客室に乗っている女性を見てみると、Cで乗船したひとが最も多く30人、Sで乗船した女性が最も多いグループで20人という結果になっています。
ここからSかCだと決めつけていいのでしょうか？
さらに乗船位置で合計してみると、Cで乗ってきた女性は71人、Sでのってきた女性は69人と拮抗しています。
他のカーネルを見てみると"S"を入力していたりします。
でも筆者はCを入れたいと思います。理由はありますが、あなた自身で調査・決断し、入力することをオススメします。

（実際にはどれほど重要なのでしょうか。全体像では2人の乗客しかいないため、モデルの精度への影響はそれほど大きくありません。ただし、この課題の主なポイントはデータ分析の実践です。もう少し詳しく疑問に思ったことを追いかけるために時間をかけてください。）

dplyrのmutateとcase_whenで値を置き換えることができます。


```R
combine <- combine %>%
  mutate(Embarked = as.character(Embarked)) %>%
  mutate(Embarked = case_when(
    is.na(Embarked) ~ "C",
    TRUE ~ Embarked
  )) %>%
  mutate(Embarked = as.factor(Embarked))
```


運賃も欠損していました。


```R
print(filter(combine, is.na(Fare)), width = Inf)
## # A tibble: 1 x 12
##   PassengerId Survived Pclass Name               Sex     Age SibSp Parch
##         <dbl> <fct>    <fct>  <chr>              <fct> <dbl> <dbl> <dbl>
## 1        1044 <NA>     3      Storey, Mr. Thomas male   60.5     0     0
##   Ticket  Fare Cabin Embarked
##   <chr>  <dbl> <chr> <fct>   
## 1 3701      NA <NA>  S
```

家族のいない60代で、三等客室に乗っている。
三等客室は運賃がほとんど決まっているので、中央値に基づいてNAを補間しましょう。

```R
med_fare_3 <- combine %>%
  filter(!is.na(Fare)) %>%
  group_by(Pclass) %>% 
  summarise(med_fare = median(Fare)) %>%
  filter(Pclass == 3) %>%
  .$med_fare

combine <- combine %>%
  mutate(Fare = case_when(
    is.na(Fare) ~ med_fare_3,
    TRUE ~ Fare
  ))
```

これでNAを補間できました。
年齢は別の方法で扱ってみましょう。

##3 変数の追加(特徴エンジニアリング)

もとからある変数よりも、信頼でき、生存を予測するための変数を作ってみましょう。
この分析操作を特徴エンジニアリングと呼びます。定義した変数はどこかにリストをつくりメモしておきましょう。
どのように定義したのか見直し、さらに検討できるように。

dplyrのmutateを使い、新しい変数を追加します。
forcatsパッケージのfct_lumpを使います。

欠損値と同じく、これらの変数を追加し、trainとtestに分割します。

```R
combine <- mutate(combine,
       fclass = factor(log10(Fare+1) %/% 1),
       age_known = factor(!is.na(Age)),
       cabin_known = factor(!is.na(Cabin)),
       title_orig = factor(str_extract(Name, "[A-Z][a-z]*\\.")),
       young = factor( if_else(Age<=30, 1, 0, missing = 0) | (title_orig %in% c('Master.','Miss.','Mlle.')) ),
       child = Age<10,
       family = SibSp + Parch,
       alone = (SibSp == 0) & (Parch == 0),
       large_family = (SibSp > 2) | (Parch > 3),
       deck = if_else(is.na(Cabin),"U",str_sub(Cabin,1,1)),
       ttype = str_sub(Ticket,1,1),
       bad_ticket = ttype %in% c('1', '5', '6', '7', '8', 'A', 'F', 'W')
       )

tgroup <- combine %>%
  group_by(Ticket) %>%
  summarise(ticket_group = n()) %>%
  ungroup

combine <- left_join(combine, tgroup, by = "Ticket") %>%
    mutate(shared_ticket = ticket_group > 1)

combine <- combine %>%
  mutate(fare_eff = Fare/ticket_group,
         title = fct_lump(title_orig, n=4),
         )

train <- combine %>% filter(!is.na(Survived))
test <- combine %>% filter(is.na(Survived))
```

新しい変数を定義します。
・fclass 運賃を対数にする。運賃の低中高を区分しやすくする。
・age_known 年齢がわかっているか(NAかどうか)
・cabin_known 上に同じ
・title_orig
・young 30よりも若いこと
・child 10よりも若いこと
・family 親子や兄弟の数を足した
・alone 親子も兄弟も0であること
・large_family 兄弟が2以上もしくは親子が3以上
・deck 客室番号がUである
・ttype チケットの1文字目を格納
・bad_ticket チケットの一文字目がc('1', '5', '6', '7', '8', 'A', 'F', 'W')である。




#3.1 年齢の影響度 年齢不明、若い、子供である
生存と、新しく追加したパラメータがどのように関係しているか確認してみましょう。



```R
p1 <- train %>%
  ggplot(aes(age_known, fill = Survived)) +
  geom_bar(position = "fill")

p2 <- train %>%
  ggplot(aes(child, fill = Survived)) +
  geom_bar(position = "fill")

p3 <- train %>%
  ggplot(aes(young, fill = Survived)) +
  geom_bar(position = "fill")

p4 <- train %>%
  ggplot(aes(Age, fill = young)) +
  geom_density(alpha = 0.5)
```


![p1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d8617a8b-20b8-b96f-527c-cae504bf6ac8.png)

![p2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5a81ae1c-d5e5-14d6-55b4-601cf4027418.png)

![p3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e2bb273d-3976-5665-7d40-36947a8203d0.png)

![p4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2c0c4fa3-8766-5bc1-58e3-1673f10c5081.png)



・child,youngという変数をみると、どちらも生存率をたかめてくれるようです。
・ただ、childという変数には欠損値が含まれています。
・youngという変数は使いやすそうな変数です。30を基準にageをわけたyoungは、密度プロットを見てわかるようにきれいに分割することが出来ています。
・年齢不明age_knownでは、年齢がわかっている乗客ほど生存率が上がりそうです。



youngと客室や性別との関係性も見てみましょう。
2変数間にfacetを使います。
書き直したりしなくていいように、plot_bar_fill_gridという関数を定義して、変数の関係を簡単に比較できるようにしよう。

```R
plot_bar_fill_grid <- function(barx, filly, gridx, gridy){
  train %>%
    ggplot(aes_string(barx, fill = filly)) +
    geom_bar(position = "fill") +
    facet_grid(reformulate(gridy,gridx))
}
```
関数をよびだして使ってみます。


```R
plot_bar_fill_grid("young", "Survived", "Sex", "Pclass")
```

![19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/14f4dc2a-4d03-8736-32b8-722af7834864.png)

・youngは女性に対しては大きな違いを生まなさそうです。
・しかし、男性では、客室とも関係してyoungが生存率について影響のある説明をしてくれそうです。(つまり男性の生存を予想するのにyoungという変数は有効？？)

さきほど作成した関数でage_knowについてもplotします。

![20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d67d663e-cf01-ea3a-cba9-65097d92139a.png)

・女性のage_knownは、客室と関係して、生存率を低下させる変数でした。(年齢がわかっていると生存率が低下)
・さらに、客室が一等、三等の男性は、年齢がわかっていると生存率があがるように見えますが、二等客室では逆のことが起こっています。

**結論：youngとage_knownは新しい変数として採用したいと思います。**

#3.2 家族関係 family alone large_family

家族関係に関するパラメータを見てみます。

```R
p1 <- train %>%
  mutate(family = as.factor(family)) %>%
  ggplot(aes(family, fill = family)) +
  geom_bar() +
  theme(legend.position = "none")

p2 <- train %>%
  ggplot(aes(alone, fill = Survived)) +
  geom_bar(position = "fill")

p3 <- train %>%
  mutate(family = as.factor(family)) %>%
  ggplot(aes(family, fill = Survived)) +
  geom_bar(position = "fill") +
  theme(legend.position = "none")

p4 <- train %>%
  ggplot(aes(large_family, fill = Survived)) +
  geom_bar(position = "fill")
```

![21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/641a9657-b8a2-ea75-7473-36cb29973626.png)

![22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/97dd1bd3-fceb-c793-f677-1a9c7edeea47.png)

![23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/ac00f911-34fd-e037-d774-1ca37004649c.png)

![24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0c361039-b9f6-2b2e-2aaa-18900d9cf7b9.png)

・一人旅が大半です。
・一人旅aloneと、大家族large_familyは生存に悪影響をしているようです。
・familyをみると1~3人の家族は生存の機会を得ることができそうです。

性別と客室と今回の変数を比べてみます。
今までのplotではなく、積み上げ棒グラフで、絶対値も確認しながら見てみます。

```R
p1 <- train %>%
  ggplot(aes(alone, fill = Survived)) +
  geom_bar(position = "stack") +
  facet_grid(Pclass ~ Sex) +
  theme(legend.position = "none")

p2 <- train %>%
  ggplot(aes(large_family, fill = Survived)) +
  geom_bar(position = "stack") +
  facet_grid(Pclass ~ Sex) +
  theme(legend.position = "none")
```

![25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b88a0023-b2e5-b421-bd02-70157ea8db86.png)

![26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b4212d47-1f86-eb8c-9676-a316c2aed507.png)

・一人旅aloneは最大派閥で、特に三等客室に多く、生存率は低い。
・興味深いことに、一等、三等の女性一人旅は生存率が高いです。まず全体的に高いですが。
・一等、二等客室の男性は生存率が低いです。

・large_familyをみると、大家族は三等客室に乗っていたようで、生存率もわずかでした。aloneと似ています。

三等客室に注目してみます。

```R
p1 <- train %>%
  filter(Pclass == 3) %>%
  ggplot(aes(alone, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Sex)

p2 <- train %>%
  filter(Pclass == 3) %>%
  ggplot(aes(large_family, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Sex)
```

![27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/3be51eea-603e-8902-74c4-5ecbb80aab04.png)

![28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e3740e13-8d2a-a366-a10c-a90c3c6fceca.png)

・男性はaloneでないこと、large_familyでないことが生存にとって有利でした。

#3.3客室番号とチケットの関係 deck, chabin_known, ttype, bad_ticket, ticet_group, shared_ticket

deckについて

```R
p1 <- train %>%
  filter(deck != "U") %>%
  ggplot(aes(deck, fill = Pclass)) +
  geom_bar(position = "dodge") +
  coord_polar() +
  #theme(legend.position = "none") +
  scale_y_log10()

p2 <- train %>%
  filter(deck != "U") %>%
  ggplot(aes(deck, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Pclass, nrow = 3)
```

![29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2ac969d2-0e44-b7e1-1474-93c93c406fc8.png)

![30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/bfabc47d-3c41-a19b-c5aa-c0684a2b7b28.png)

deckにはcabinの頭文字をいれました。

・円グラフの背景は対数半径軸になっていることに注意です。
・deckがB,Cに最も人が居ます。
・deckとPclassの関係がわかります。一等客室はA,B,Cに存在しています
・deck=BはA,Cよりも生存率が高いです。
・三等客室でも、deck=Eの人は全員生存しましたが、Eにいたのは3人だけでした。
・二等客室の旅行者は、D,E,Fにいましたが、生存率に差はありませんでした。ただ、どのdeckでも二等客室は生存率が高いです。

cabin_knownについて

```R
p1 <- train %>%
  mutate(cabin_known = fct_recode(cabin_known, F = "FALSE", T = "TRUE")) %>%
  ggplot(aes(cabin_known, fill = Survived)) +
  geom_bar(position = "dodge") +
  facet_grid(Sex ~ Pclass) +
  scale_y_log10() +
  theme(legend.position = "none")

p2 <- train %>%
  mutate(cabin_known = fct_recode(cabin_known, F = "FALSE", T = "TRUE")) %>%
  ggplot(aes(cabin_known, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_grid(Sex ~ Pclass) +
  theme(legend.position = "bottom")
```

![31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/8aca94c1-66bf-ca4a-5462-06d5ce05f0b2.png)

![32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/54b2371a-2233-ea0c-4329-3aea4a61880e.png)

・三等客室の男性はわずかながら生存率を上げることができています。
・cabin_knowとの関係は、一等客室であり、女性であることが生存に影響を与えるようです。

ttypeについて

```R
p1 <- train %>%
  ggplot(aes(ttype, fill = ttype)) +
  geom_bar() +
  theme(legend.position = "none") +
  facet_wrap(~ Pclass, nrow=3)

av_surv <- train %>%
  group_by(Pclass, Survived) %>%
  count() %>%
  spread(key = Survived, value = n) %>%
  mutate(frac = `1`/(`0`+`1`))

p2 <- train %>%
  ggplot(aes(ttype, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Pclass, nrow = 3) +
  geom_hline(data = av_surv, aes(yintercept = frac), linetype=2)
```

![33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0e06a6c6-d6ad-7670-8664-ee5398001c6d.png)

![34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/bcb0e0e1-77c3-27cc-6a82-63c46312f6e7.png)

・ttypeには16種類があり、最も出現頻度の高いものは客室階級と強い関係を持っています。
・1とPは一等客室と、2は二等客室と三等客室に強く関係している。
・C,Sは二等客室と三等客室に同じ程度の強さで関係を持っている。
・最も頻度の高いttypeは客室の生存率と一致します(点線で示しています)
・生存率の高いttype、例えば一等客室の2, 三等客室の9などは、偶然助かったケースと考えられるかもしれません。
・例外として、三等客室のttypeが2であった人は、50件確認されており、生存率の上昇について確かに関係している値かもしれません。
・三等客室のAは件数がおおいけれど生存率は低そうです。

**結論：ttypeを影響力のある変数として採用します。**


bad_ticketについて
前回のttypeについて、bad_ticketというものを定義します。
生存確率の低かった1,5,6,7,8,A,F,Wです。
この"生存率の低さ"というのは、小さな集団に分けてしまって出てきた統計量です。
単なるばらつきである可能性が高いです。
ですが、チケット番号と生存をうまく結び付けようと思ってもこれ以上うまく結びつかないかもしれません。
もちろん、チケットが生存について悪い影響を与えるわけでなく、チケットに基づいて割り振られた船内の場所が生存率に影響を与えていたりします。
二項誤差を計算して表示します。
95%信頼を計算するために短い関数をつかいます。
**get_binCIがカーネル内に見当たりませんでしたので、以下のものを代用しています**

```R
get_binCI <- function(x,n) {
    rbind(setNames(c(binom.test(x,n)$conf.int),c("lwr","upr")))
}

train %>%
  mutate(bad_ticket = factor(bad_ticket)) %>%
  group_by(bad_ticket, Survived) %>%
  count() %>%
  spread(Survived, n, fill = 0) %>%
  mutate(frac_surv = `1`/(`1`+`0`)*100,
         lwr = get_binCI(`1`,(`1`+`0`))[[1]]*100,
         upr = get_binCI(`1`,(`1`+`0`))[[2]]*100,
         ) %>%
  ggplot(aes(bad_ticket, frac_surv, fill = bad_ticket)) +
  geom_col() +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.5, size = 1) +
  labs(y = "Survival fraction") +
  theme(legend.position = "none")
```

![35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/43c28ae0-15ea-b32d-d715-b133fea6a3fe.png)


・bad_ticketで分類した人たちは20~30%程度の生存率であったのに対し、
そのほかのグループでは50%ほどが生存した。
・ttypeのときのように、チケット番号にはモデル作成に使えそうなものがありそうだ。

**結論：ttypeと同じようにbad_tiketがどう動くのか確認します**

ticket_groupとshared_ticketについて


同じチケットで乗船している異なる人がいることがわかります。
チケットによってグループ化し、チケットをシェアしている人を割り当てます。ticket_group
さらに、乗客が自分のチケットを誰かと共有しているか判別するshared_ticketをつくります。

```R
train %>%
  arrange(Ticket) %>%
  select(Ticket, ticket_group, shared_ticket, Name) %>%
  head(9) %>%
  tail(-3)

## # A tibble: 6 x 4
##   Ticket ticket_group shared_ticket Name                                  
##   <chr>         <int> <lgl>         <chr>                                 
## 1 110413            3 T             Taussig, Mr. Emil                     
## 2 110413            3 T             Taussig, Mrs. Emil (Tillie Mandelbaum)
## 3 110413            3 T             Taussig, Miss. Ruth                   
## 4 110465            2 T             Porter, Mr. Walter Chamberlain        
## 5 110465            2 T             Clifford, Mr. George Quincy           
## 6 110564            1 F             Bjornstrom-Steffansson, Mr. Mauritz H…
```

チケット「110413」は3人の「Taussigs」によって共有され、チケット「110465」は「Mr. PorterとMr. Clifford」によって共有されています。チケット「110564」は「Mr. Bjornstrom-Steffansson」によってのみ使用されました。

生存率とのplotをしましょう。


```R
p1 <- train %>%
  group_by(Survived, shared_ticket) %>%
  count() %>%
  ggplot(aes(shared_ticket, n, fill = Survived)) +
  geom_col(position = "dodge") +
  geom_label(aes(label = n), position = position_dodge(width = 1)) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

p2 <- train %>%
  ggplot(aes(ticket_group, fill = Survived)) +
  geom_bar(position = "dodge") +
  theme(legend.position = "none") +
  scale_y_log10()

p3 <- train %>%
  ggplot(aes(shared_ticket, fill = Survived)) +
  geom_bar(position = "fill") +
  facet_wrap(~ Pclass) +
  theme(legend.position = "none")

p4 <- train %>%
  ggplot(aes(ticket_group, fill = Fare)) +
  stat_summary(fun.data = "mean_cl_boot", col = "red")
```

![36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/94cc6189-682b-9cfd-dcc4-51faddfe5684.png)

![37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/978ad664-eb06-1ba8-c112-5fb154f87207.png)

![38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/739ef650-7d4e-d0ac-5857-0a6ef38ccf26.png)

P4がplotできない
library(Hmisc)が抜けているみたい。

```R
library(Hmisc)
p4 <- train %>%
  ggplot(aes(ticket_group, Fare)) +
  stat_summary(fun.data = "mean_cl_boot", col = "red")
```
![39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/911f96f1-5614-128f-0985-2ea9ee2a64a0.png)

・チケットのシェアは生存率にいい影響を与えます。
・一番上のplotでは縦軸をなくし、それぞれにカウントを記入したデザインにしています。
・チケット共有人数は2~4人で分けているグループが生存率を高めていそうです。
・一等客室の乗客は三等客室と比較して、チケットシェアによる生存率が高くなりそうです。
・ticket_groupと運賃を見てみると興味深いことがわかります。
・ticket_groupの人数が多いほど、運賃が上昇していることがわかります。
・ただし、4人がピークであり、4人以上では発生件数が少ないために不確実性が大きくなります。
・6人を迎えると運賃は下がり始めます。

つづく↓


[kaggleでデータ分析の基礎 #4【titanic, tidyr, dplyr,ggplot】](https://qiita.com/H_Y_J/items/ae6f5c30fb58b7d78240)