kaggleのtitanicで0.81340を出した話

はじめに

みんな大好き（？）kaggleのtitanicコンペ

Titanic: Machine Learning from Disaster | Kaggle

で0.81340のスコアを出したので色々まとめてみましたという記事です。

このコンペはtitanic号が沈没したときのデータを使用して、乗客が生存したか死亡したかを予測するモデル構築し、そのモデルの予測精度を競うコンペです。
一般に80%の精度を超えると良いとされているコンペのようですので、この記事がそこを目指している方々への参考となればよいかと思います。

ここで今回使うパッケージを列挙しておきます。予めインストールしておくといいかもしれません。

{dplyr}
{tidyr}
{dummies}
{ranger}
{e1071}

dplyrとtidyrはデータ整形で使用します。
dummiesは整形したデータをデミー変数化するときに使用します。
rangerはランダムフォレストを実行できるパッケージで、欠損値の補完や変数選択で使用します。
e1071はSVMを実行できるパッケージで、予測モデルを構築するときに使用します。

それでは、レッツデータ解析((´^ω^)

集計関数の準備

これからデータ解析をするわけですが、データをみていくにあたって便利な関数を作っておきます。

# Survivedと変数xのクロス集計の右側に
# 変数x内のカテゴリー比率を加えたものを表示する関数
S_table <- function(x, survived){
  tbl <- table(x, survived)
  row_sum <- apply(tbl, 1, sum)
  s_ratio <- tbl[,2]/row_sum
  
  return(cbind(tbl, s_ratio))
}

# クロス割合を表示する関数
# (クロス集計表の要素が割合)
rate_table <- function(x, y){
  tbl <- table(x, y)
  row_sum <- apply(tbl, 1, sum)
  s_ratio <- tbl/row_sum
  
  return(s_ratio)
}

S_tableはクロス集計をする際に、行についての割合もテーブルの右側に表示する関数です。
今回はtitanicの生存・死亡を表すカラム「Survived」と他のカテゴリカル変数とのクロス集計でこの関数を用いるので、S_tableと名付けました。

rate_tableはクロス集計のテーブルの要素が割合になったものです。

データの読み込み

データを読み込みます

# データ読み込み
d <- read.csv("train.csv")
d_t <- read.csv("test.csv")

データの前処理

titanicのデータ次のような形をしています：

ひとつひとつ見ていきましょう。

PassengerId: 乗客のID
Survived: 生存(1)、死亡(0)
Pclass: 乗客の階級
Name: 乗客の名前
Sex: 性別
Age: 年齢
SibSp: タイタニック号に乗っていた兄弟、姉妹、義兄弟、義姉妹、夫、妻の数（自分を除く）
Parch: タイタニック号に乗っていた母親、父親、息子、娘の数
Ticket: チケットナンバー
Fare: 乗船料金
Cabin: キャビン番号
Embarked: 乗船場

まずNameに着目すると

「名字, 敬称. 名前」

といった順にデータが入っています。
後に確認するAgeには多くの欠損値を含んでいるのですが、この敬称はこのAgeの欠損値を補完するのに役立ちますのでNameのカラムから敬称だけを取り出して新たな変数とします。
例えば敬称「Master」は子供につけられる敬称です。ここからAgeの欠損が補完できそうです。

またある方がブログで「自分以外の家族の生死がその人の生死に影響を与えた」という仮説を立てていて、ふむなるほどなと思いました。ただ自分以外の家族の生死を変数とするために、どのひとがどの人家族関係にあるのかを知る必要があります。そこで名字が同じ人は家族だということにすればとりあえず良いことにして、Nameのカラムから名字だけを取り出して新たな変数とします。

つぎにCabinを見てみます。
Cabinは部屋の情報を表した変数ですが、凄まじいほどの欠損値があります。ただなんとなく生データを眺めてみると、Cabinの情報がある乗客はSurvived=1となっている傾向がある、つまり生存しやすい傾向がありそうな気がします。このことは後に確認してみるとして、Cabinを使える変数にするために、欠損していないものについては先頭のアルファベットを抽出し、欠損しているものについてはUnknownの頭文字である「U」を
入れてあげて新たな変数とすることにします。

以上の操作をするコードが以下になります：

# Cabinの処理
train <- d %>% 
  separate("Cabin", into=c("Cabin1", "Cabin2", "Cabin3", "Cabin4"), sep=" ") %>% 
  mutate(Cabin1 = substr(Cabin1, 1, 1),
         Cabin2 = substr(Cabin2, 1, 1),
         Cabin3 = substr(Cabin3, 1, 1),
         Cabin4 = substr(Cabin4, 1, 1)) %>% 
  mutate(Cabin1 = if_else(Cabin1=="", "U", Cabin1),
         Cabin2 = if_else(is.na(Cabin2), "U", Cabin2),
         Cabin3 = if_else(is.na(Cabin3), "U", Cabin3),
         Cabin4 = if_else(is.na(Cabin4), "U", Cabin4))

# Nameの処理
train <- train %>% 
  separate("Name", into=c("Last_Name", "Title"), sep=",") %>% 
  mutate(Title = gsub("\\..+$", "", Title)) %>% 
  mutate(Title = gsub(" ", "", Title))

ここでCabinについて、例えば以下のように乗客一人について最大で4つのキャビン情報が付与されているデータがあります：

C○○○　CXXX　CYYYY　CZZZZ

ですのでこの時点ではそれぞれ前からCabin1, Cabin2, Cabin3, Cabin4と名前をつけて、カラムを追加しています。
（ちなみに、後の分析で最初のCabin情報、つまりCabin1しかSurvivedの予測に影響を与えないことがわかりましたので、最終的にはCabin1しか使いません。）

基礎集計と分布の確認

なんといってもまずは基礎集計と分布の確認です。

###################################
# 基礎集計や分布の確認
###################################
# 全体の生存率
# 0.3838384
sum(train$Survived)/nrow(train)

基本的にこの全体の生存率0.3838384が、ある変数が生存・死亡に影響を与えているのかを判断する数字になります。

次から、この記事の最初で作ったS_table関数やggplotの可視化を使って基礎集計や分布の確認を行っていきます。
以下がそのコードです。

library(ggplot2)


# 2変数の関係の確認
# 対 Survived
# Pclass
S_table(train$Pclass, train$Survived)

# Sex
S_table(train$Sex, train$Survived)

# Title
S_table(train$Title, train$Survived)

# Embarked
S_table(train$Embarked, train$Survived)

# SibSp
S_table(train$SibSp, train$Survived)

# Parch
S_table(train$Parch, train$Survived)

# Cabin1
S_table(train$Cabin1, train$Survived)
S_table(train$Family_size, train$Survived)

# 気になった変数
table(train$Cabin1, train$Pclass)
rate_table(train$Embarked, train$Pclass)
rate_table(train$Embarked, train$Cabin1)

#Age (散布図)
plot(train$Age, train$Survived+runif(nrow(train), -0.3, 0.3))

# Age (密度トレース)
dat <- train %>% 
  select(Survived, Age) %>% 
  mutate(Survived = as.factor(Survived))
g <- ggplot(dat, aes(Age, colour=Survived, fill=Survived, alpha=0.5)) +
  geom_density()
plot(g)


# Fare (散布図)
plot(train$Fare, train$Survived+runif(nrow(train), -0.3, 0.3))

# Fare (密度トレース)
dat <- train %>% 
  select(Survived, Fare) %>% 
  mutate(Survived = as.factor(Survived))
g <- ggplot(dat, aes(Fare, colour=Survived, fill=Survived, alpha=0.5)) +
  geom_density()
plot(g)


# PclassとFareの密度トレース
dat <- train %>% 
  select(Pclass, Fare) %>% 
  mutate(Pclass = as.factor(Pclass))
g <- ggplot(dat, aes(Fare, colour=Pclass, fill=Pclass, alpha=0.5)) +
  geom_density()
plot(g)




# 3変数の関係の確認
# SurvivedとFareとPclass
plot(train$Fare, train$Pclass+runif(nrow(train), -0.3, 0.3), 
     col=ifelse(train$Survived == 1, "red", "blue"))

# SurvivedとFareとSex
plot(train$Fare, train$Sex+runif(nrow(train), -0.3, 0.3), 
     col=ifelse(train$Survived == 1, "red", "blue"))

# SurvivedとFareとPclass
dat <- train %>% 
  select(Survived, Pclass, Fare) %>% 
  mutate(Survived = as.factor(Survived),
         Pclass = as.factor(Pclass))
g <- ggplot(dat, aes(Fare, colour=Pclass, fill=Pclass, alpha=0.5)) +
  geom_density() +
  facet_wrap(~Survived, nrow=2) 
plot(g)

例えば最初の集計である

S_table(train$Pclass, train$Survived)

の結果は次のようになります：

縦軸がPclass、横軸がSurvivedを表しています。
この結果をみると、Pclassが1、2、3の順で生存率が高い、つまり乗客の階級が高い人ほど生存率が高い傾向が読みとれます。

また、生存・死亡別のAgeの分布を確認してみます。

# Age (密度トレース)
dat <- train %>% 
  select(Survived, Age) %>% 
  mutate(Survived = as.factor(Survived))
g <- ggplot(dat, aes(Age, colour=Survived, fill=Survived, alpha=0.5)) +
  geom_density()
plot(g)

0〜10歳の部分でSurvived=0の分布よりもSurvived=1の分布の方が密度か高くなっているのがわかります。
どうやら子供は優先的に救助されたみたいですね。

他のコードも実行した結果、以上の基礎集計から次のことがわかりそうです：

・Pclassが良いほど生存する傾向

・女性の方が男性に比べて圧倒的に生存する傾向

・SibSp=0は平均的な生存率
　SibSp=1,2の順に生存率が高い
　それ以降は低い生存率
　
・Parch=0は平均的な生存率
　Parch=1,2の順に生存率が高い
　Parch=3は生存率は高いがサンプル数が少ない（5）
　それ以降は低い生存率(サンプル数も少ない)
　
・Family_size = SibSp + Parch + 1　を本人を含めた家族数としたときに
　Family_size = 1の人、つまり一人で乗っていたひとは平均よりやや低い生存率
　Family_size = 2,3,4 の人たちは高い生存率
　それ以降では低い生存率
　

　家族の規模が生存に影響しているかも
　

・Embarkedについて、Cから出発した人の生存率が高い
　それ以外は平均的
　
　疑問
　なんでただ出発地点が生存に影響を与えているのか？
　
　⇓
　
　PclassとEmbarkedのクロス集計をすると、Cから出発した人はPclass=1の方が多い
　お偉いさんはCから乗った？
　
・Ageについて、子供ほど生存しやすい傾向

・Fareについて、高いお金を払っている人の方が生存しやすい傾向
　Fare=30~50以上でその傾向が出始めている感じ
　
・Cabinについて、
　 欠損値多い
　 そもそもCabinの情報が残っている人の生存率が高い
　 PclassとCabinをクロス集計すると、Cabinの情報が残っている人はPclass=1の人が多い。お偉いさんのCabin情報だけは予めとってて、他の乗客のCabin情報はとってなかった？
　

基礎集計と分布の確認を踏まえて、データを整形

基礎集計と分布の確認を踏まえて、データを整形を整形していきます。

Ageの欠損値処理について

PclassとTitle(敬称)ごとの平均値で補完をしました。
それでも補完できなかった部分が１つできてしまい、そのデータがSibSp=0でParch=0の人だったため、SibSp=0でParch=0の人の平均値でもって補完をしました。

Fareの欠損値処理について

Pclass=3の人が欠損をしていたので、おそらく低い乗船料だったであろうということで0.0で補完をしました。
その後正規化（平均0、分散1に調整）をしています。

Embarkedの欠損値処理について

2つ欠損値がありましたが、一番数の多かった「S」を入れています。

自分以外の家族の生存率について

説明が面倒なのでコードをみて理解してください。fam_suv_rateという名前をつけています。
全家族分の生存率を出せるわけではないので、欠損した部分についてはランダムフォレストでfam_suv_rateの予測モデルを作り、そのモデルから得られる予測値で欠損値を埋めています。

それ以外の説明はしません。汚いコードで申し訳ないですが、頑張って解読してください。

※今回は予測モデルを構築する際にSVMを使用しますのでデータを予めダミー変数化します。その際に昔の記事

kefism.hatenablog.com

で紹介したデータフレームを渡すとすべてのカラムについてダミー変数化してくれる関数、convDummiesを使用します。
ここでconvDummiesのコードを再掲します：

convDummies <- function(data, is.drop = FALSE){
  library(dummies)
  
  N <- ncol(data)
  row_names <- names(data)
  
  names_list <- c()
  new_data <- rep(NA, nrow(data))
  for(n in 1:N){
    unique_value <- sort(unique(data[,n]))
    dummied_data <- dummy(data[,n])
    
    if(is.drop == TRUE){
      new_data <- cbind(new_data, dummied_data[,-ncol(dummied_data)])
      names_list <- c(names_list, 
                      paste(row_names[n], unique_value, sep = ".")[-ncol(dummied_data)])
    } else {
      new_data <- cbind(new_data, dummied_data)
      names_list <- c(names_list, paste(row_names[n], unique_value, sep = "."))
    }
  }
  
  new_data <- as.data.frame(new_data)
  names(new_data) <- c("temp", names_list)
  
  return(new_data[,-1])
}

以下からデータの処理になります。

# データ整形関数
createData <- function(){
  library(dplyr)
  library(tidyr)
  library(ranger)
  
  train <- read.csv("train.csv")
  test <- read.csv("test.csv")
  
  train$Embarked[train$Embarked == ""] <- "S"
  
  Survived <- train$Survived
  
  X <- train %>% 
    select(-Survived) %>% 
    bind_rows(test) %>% 
    select(-PassengerId, -Ticket)
  
  # Cabinの処理
  X <- X %>% 
    separate("Cabin", into=c("Cabin1", "Cabin2", "Cabin3", "Cabin4"), sep=" ") %>% 
    mutate(Cabin1 = substr(Cabin1, 1, 1),
           Cabin2 = substr(Cabin2, 1, 1),
           Cabin3 = substr(Cabin3, 1, 1),
           Cabin4 = substr(Cabin4, 1, 1)) %>% 
    mutate(Cabin1 = if_else(Cabin1=="", "U", Cabin1),
           Cabin2 = if_else(is.na(Cabin2), "U", Cabin2),
           Cabin3 = if_else(is.na(Cabin3), "U", Cabin3),
           Cabin4 = if_else(is.na(Cabin4), "U", Cabin4)) %>% 
    select(-c(Cabin2, Cabin3, Cabin4))
  
  X <- X %>% 
    mutate(is_Cabin = if_else(Cabin1 == "U", "No", "Yes"))
  
  # Nameの処理
  Name_list <- c("Master", "Miss", "Mr", "Mrs", "Rev")
  X <- X %>% 
    separate("Name", into=c("Last_Name", "Title"), sep=",") %>% 
    mutate(Title = gsub("\\..+$", "", Title)) %>% 
    mutate(Title = gsub(" ", "", Title)) %>% 
    mutate(Title = if_else(Title %in% Name_list, Title, "Otherwise"))
  
  # Family_sizeの追加
  X <- X %>% 
    mutate(Family_size = SibSp + Parch +1)
  
  # Familly_sizeで家族をカテゴリー分け
  X <- X %>% 
    mutate(Family_type = if_else(Family_size == 1, "singleton",
                                 if_else(2 <= Family_size & Family_size <= 4, "middle", "large")))
  
  # Ageの欠損値補完
  age_for_na <- 
    na.omit(X) %>% 
    group_by(Pclass, Title) %>% 
    summarise(age_ave = mean(Age))
  
  X <- X %>% 
    left_join(age_for_na) %>% 
    mutate(Age = if_else(is.na(Age), age_ave, Age)) %>% 
    select(-age_ave)
  
  X$Age[is.na(X$Age)] <- mean(X$Age[!is.na(X$Age) & X$SibSp == 0 & X$Parch == 0])
  
  X <- X %>%
    mutate(Age_desc = if_else(0 <= Age & Age <= 6, "0_6",
                              if_else(7 <= Age & Age <=10, "7_10",
                                      if_else(11 <= Age & Age <= 15, "11_15",
                                              if_else(16 <= Age & 20 <= Age, "16_20",
                                                      if_else(21 <= Age & Age <= 30, "21_30", "30_"))))))
  
  # Fareの欠損値補完と正規化
  X$Fare[is.na(X$Fare)] <- 0.0
  X$Fare <- (X$Fare - mean(X$Fare))/sd(X$Fare)
  
  # カテゴリー変数のダミー化
  df_category <- X %>% 
    select(Sex, Title, Pclass, Cabin1, is_Cabin, Family_type, Age_desc, Embarked)
  X_continuous <- X %>% 
    select(-c(Sex, Title, Pclass, Cabin1, is_Cabin, Family_type, Age_desc, Embarked))
  
  X_dummy <- convDummies(df_category, is.drop = TRUE)
  
  X <- cbind(X_continuous, X_dummy)
  
  
  # trainデータとtestデータに分割
  train_new <- X[1:891,]
  train_new <- data.frame(Survived = Survived) %>% 
    cbind(train_new)
  test_new <- X[892:1309,]
  
  # 自分以外の家族の生存割合を追加
  fam_suv <- train_new %>% 
    group_by(Last_Name) %>% 
    summarise(f_n = n(), s_n = sum(Survived))
  
  train_new <- train_new %>% 
    left_join(fam_suv) %>% 
    mutate(fam_suv_rate = if_else(f_n==1, 0, (s_n-Survived)/(f_n-1))) %>% 
    select(-f_n, -s_n, -Last_Name)
  
  test_new <- test_new %>% 
    left_join(fam_suv) %>% 
    mutate(fam_suv_rate = if_else(Family_size==1, 0, s_n/f_n)) %>% 
    select(-f_n, -s_n, -Last_Name)
  
  fam_suv_df <- train_new %>% 
    select(-Survived)
  fam_suv_df_na <- test_new[is.na(test_new$fam_suv_rate),] %>% 
    select(-fam_suv_rate)
  fam_suv.fit <- ranger(fam_suv_rate~., fam_suv_df[,29:31])
  pred_fam_suv <- predict(fam_suv.fit, fam_suv_df_na)
  test_new$fam_suv_rate[is.na(test_new$fam_suv_rate)] <- pred_fam_suv$predictions
  
  
  
  
  train_test <- list()
  train_test[[1]] <- train_new
  train_test[[2]] <- test_new
  
  return(train_test)
}

ここで、あらたに追加したfam_suv_rateが乗客の生存に影響を与えていそうかをちょっと確認してみましょう。

# 先程作成したデータ前処理関数で前処理をする。
data <- createData()
train <- data[[1]]
test <- data[[2]]


# Surviveとfam_sv_rateの散布図
plot(train$fam_suv_rate+runif(nrow(train), -0.02, 0.02),
     train$Survived+runif(nrow(train), -0.3, 0.3))

# Surviveとfam_sv_rateの密度トレース
dat <- train %>% 
  select(Survived, fam_suv_rate) %>% 
  mutate(Survived = as.factor(Survived))
g <- ggplot(dat, aes(fam_suv_rate, colour=Survived, fill=Survived, alpha=0.5)) +
  geom_density()
plot(g)

密度トレースの方のプロットは以下のようになります：

Survived=1のときのfam_suv_rateの密度とSurvived=0のときのfam_suv_rateを比べてみますと、fam_suv_rateが高い部分でSurvived=1の場合の密度のほうがSurvived=0の場合の密度よりも高くなっています。fam_suv_rateはSurvivedの予測に影響を与えるかもしれません。

予測モデルの構築

さあ、いよいよ予測モデルの構築です。
今回はSVMを使用します。

変数選択

まず、色々作った変数の選択をするために、ランダムフォレストを使用して変数の重要度をみてみます。

# SVMモデルを構築するためのライブラリ
library(e1071)

# 先程作成したデータ前処理関数で前処理をする。
data <- createData()
train <- data[[1]]
test <- data[[2]]


# 変数の重要度をみる
rf.fit <- ranger(as.factor(Survived)~., train, num.trees=4000,
                 importance = 'impurity')

# 変数の重要度を可視化
pd <- data.frame(Variable = names(rf.fit$variable.importance),
                 Importance = as.numeric(rf.fit$variable.importance)) %>% 
  arrange(Importance)
p <- ggplot(pd, aes(x=factor(Variable, levels=unique(Variable)), y=Importance)) +
  geom_bar(stat="identity") +
  xlab("Variables") + 
  coord_flip()
plot(p)

ふむ、なるほど...

プロットをみて、Sex.female〜Parchまでの変数を使用することにしました。

ハイパーパラメータのチューニング

SVMのハイパーパラメータcostとgammaをチューニングします。

# パラメータチューニング
# Cross Validation
svm.df <- train[(c("Survived", as.character(pd$Variable[17:31])))]
svm.df_t <- test[as.character(pd$Variable[17:31])]


tune_res <- tune.svm(as.factor(Survived)~., data=svm.df)
tune_res$best.model

cost: 1
gamma: 0.06666667

が良いとのことです。

予測

cost: 1
gamma: 0.06666667

のSVMで予測モデルを構築し、テストデータで予測値を出します。
その後、kaggleに投稿できる形にデータを整え、csvファイルとして書き出します。

svm.fit <- svm(as.factor(Survived)~., data=svm.df, cost=1, gamma=0.06666667)

# 予測
pred <- predict(svm.fit, svm.df_t)
prediction <- data.frame(PassengerId = 892:1309,
                         Survived = pred)

write.csv(prediction, "prediction.csv", row.names = FALSE)

kaggleへ投稿

やったぜ。