最近Rに入門したので、KaggleのTitanicをRで実践してみました。予測モデルはロジスティック回帰とランダムフォレストのそれぞれで構築し、Submissionまで一通り実施しています。

パッケージの読込み

In

library('randomForest')  # ランダムフォレスト
library('na.tools')      # 欠損値の処理
library('Metrics')       # Accuracy算出用

データの読込み

In

train <- read.csv(
  'data/titanic/train.csv', 
  header=T, 
  stringsAsFactors=F, 
  na.strings=(c("NA", "")))

test <- read.csv(
  'data/titanic/test.csv', 
  header=T, 
  stringsAsFactors=F, 
  na.strings=(c("NA", "")))

sample <- read.csv(
  'data/titanic/gender_submission.csv', 
  header=T, 
  stringsAsFactors=F, 
  na.strings=(c("NA", "")))

EDA （探索的データ分析）

いくつかのカラムについて、簡単なEDAを行ってみます。

性別 (Sex)

性別と生存に関連性があるかどうかを調べるため、独立性の検定を行います。

# Sex列のユニーク値を確認
unique(train$Sex)
# Sex列を数値に変換したfemale列を作成
train$female <- ifelse(train$Sex=='female', 1, 0)
test$female <- ifelse(test$Sex=='female', 1, 0)

train_f <- train[train$Sex == 'female', ]
train_m <- train[train$Sex == 'male', ]

N_fa <- nrow(train_f[train_f$Survived == 1, ])
N_fd <- nrow(train_f[train_f$Survived == 0, ])
N_ma <- nrow(train_m[train_m$Survived == 1, ])
N_md <- nrow(train_m[train_m$Survived == 0, ])

# クロス集計表変数の作成
X <- matrix(c(N_fa, N_fd, N_ma, N_md), 2, 2)
X

# 独立性の検定
chisq.test(X)

有意水準を0.05として$\chi^2$検定を行うと、検定結果のp値は有意水準を大きく下回っており、性別と生存には関連性があると言えます。

Out

[1] "male"   "female"
     [,1] [,2]
[1,]  233  109
[2,]   81  468

	Pearson's Chi-squared test with Yates' continuity correction

data:  X
X-squared = 260.72, df = 1, p-value < 2.2e-16

年齢 (Age)

年齢と生存の関連性を調べるため、ヒストグラムを確認します。

In

train_alive <- train[train$Survived == 1, ]
train_dead  <- train[train$Survived == 0, ]

hist_1 <- hist(train_alive$Age, breaks=seq(0, 80, 1), xlab='Age', ylim=c(0,25), col ="#ff00ff40")
hist_2 <- hist(train_dead$Age,  breaks=seq(0, 80, 1), xlab='Age', ylim=c(0,25), col ="#0000ff40", add=T)

傾向として10歳未満の子供は生存確率が高く、10代後半~20代前後の若者が多く死亡していることが分かります。

チケットクラス (Pclass)

チケットクラスと生存の関連性を調べるため、各クラスの生存者数・死亡者数を棒グラフで表示します。また、クロス集計表を作成して独立性の検定を行います。

In

train_pc1 <- train[train$Pclass == 1, ]
train_pc2 <- train[train$Pclass == 2, ]
train_pc3 <- train[train$Pclass == 3, ]

N_pc1_a <- nrow(train_pc1[train_pc1$Survived == 1, ])
N_pc1_d <- nrow(train_pc1[train_pc1$Survived == 0, ])
N_pc2_a <- nrow(train_pc2[train_pc2$Survived == 1, ])
N_pc2_d <- nrow(train_pc2[train_pc2$Survived == 0, ])
N_pc3_a <- nrow(train_pc3[train_pc3$Survived == 1, ])
N_pc3_d <- nrow(train_pc3[train_pc3$Survived == 0, ])

# クロス集計表変数の作成
X <- matrix(c(N_pc1_a, N_pc1_d, N_pc2_a, N_pc2_d, N_pc3_a, N_pc3_d), 2, 3)
# 棒グラフの表示
barplot(X, col = c("#0000ff40", "#ff00ff40"), xlab = 'Pclass', names = c(1,2,3))

# 独立性の検定
chisq.test(X)

有意水準を0.05とするとき、p値は有意水準を大きく下回っており、関連性があるといえます。

Out

	Pearson's Chi-squared test

data:  X
X-squared = 102.89, df = 2, p-value < 2.2e-16

欠損値の処理

予測モデルを構築するにあたり、本記事ではチケットクラス (Pclass), 性別 (female), 年齢 (Age), 兄弟・配偶者の数 (SibSp) を使用します。よってこれらの列に含まれる欠損値を処理していきます。

まず欠損値の有無・数を確認します。

In

na_count_train <- sapply(train, function(y) sum(is.na(y)))
na_count_train
na_count_test <- sapply(test, function(y) sum(is.na(y)))
na_count_test

Out

PassengerId    Survived      Pclass        Name         Sex         Age       SibSp       Parch      Ticket        Fare       Cabin    Embarked      female 
          0           0           0           0           0           0           0           0           0           0         687           2           0 
PassengerId      Pclass        Name         Sex         Age       SibSp       Parch      Ticket        Fare       Cabin    Embarked      female 
          0           0           0           0           0           0           0           0           1         327           0           0

使用する列のうち年齢 (Age) のみ欠損値を含んでいるので、これを平均値で埋めることとします。

In

age_mean = mean(na.rm(train$Age))
train$Age[is.na(train$Age)] <- age_mean
test$Age[is.na(test$Age)] <- age_mean

モデルの構築・予測

目的変数の型変換

生存・死亡 {0, 1} を予測する2値分類問題のため、目的変数の型を数値型からfactor型に変換します。これを行わない場合、ランダムフォレストにおいて回帰問題として扱われてしまいます。

In

train$Survived <- as.factor(train$Survived)

学習用データの分割

学習用データを、予測モデルの学習用と検証用に分割します。ここでは学習用を全体の80%、検証用を残り20%としています。

In

N <- nrow(train)
idx <- sample(N, N*0.8)
titanic.train.data <- train[idx, ]
titanic.valid.data <- train[-idx, ]
titanic.test.data <- test

ロジスティック回帰

1. モデルの構築

stats パッケージに含まれる glm() 関数を使用し、分布族には二項分布 (binomial) を指定します。

In

titanic.lr.model <- glm(Survived ~ Pclass+female+Age+SibSp, data=titanic.train.data, family=binomial)
summary(titanic.lr.model)

Out

Call:
glm(formula = Survived ~ Pclass + female + Age + SibSp, family = binomial, 
    data = titanic.train.data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0662  -0.5698  -0.4180   0.5833   2.4552  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.570297   0.496878   5.173  2.3e-07 ***
Pclass      -1.186926   0.137520  -8.631  < 2e-16 ***
female       2.817840   0.220202  12.797  < 2e-16 ***
Age         -0.043851   0.009061  -4.839  1.3e-06 ***
SibSp       -0.276739   0.120190  -2.303   0.0213 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 946.06  on 711  degrees of freedom
Residual deviance: 615.81  on 707  degrees of freedom
AIC: 625.81

Number of Fisher Scoring iterations: 5

2. モデルの検証

In

titanic.lr.valid.pred1 <- predict(titanic.lr.model, newdata = titanic.valid.data)
titanic.lr.valid.pred2 <- 1/(1+exp(-titanic.lr.valid.pred1))
titanic.lr.valid.pred3 <- array(ifelse(titanic.lr.valid.pred2 >= 0.5, 1, 0))

accuracy(titanic.valid.data$Survived, titanic.lr.valid.pred3)
# [1] 0.7597765

3. テストデータによる予測

In

titanic.lr.test.pred1 <- predict(titanic.lr.model, newdata = titanic.test.data)
titanic.lr.test.pred2 <- 1/(1+exp(-titanic.lr.test.pred1))
titanic.lr.test.pred3 <- array(ifelse(titanic.lr.test.pred2 >= 0.5, 1, 0))

titanic.lr.result <- data.frame(titanic.test.data$PassengerId, titanic.lr.test.pred3)
names(titanic.lr.result) <- c('PassengerId', 'Survived')
write.csv(titanic.lr.result, 'resul_logistic-regression.csv', row.names=FALSE)

Public Score は 0.75358 でした。

ランダムフォレスト

1. モデルの構築

In

titanic.rf.model <- randomForest(Survived ~ Pclass+female+Age+SibSp, data=titanic.train.data)
summary(titanic.rf.model)

Out

                Length Class  Mode     
call               3   -none- call     
type               1   -none- character
predicted        712   factor numeric  
err.rate        1500   -none- numeric  
confusion          6   -none- numeric  
votes           1424   matrix numeric  
oob.times        712   -none- numeric  
classes            2   -none- character
importance         4   -none- numeric  
importanceSD       0   -none- NULL     
localImportance    0   -none- NULL     
proximity          0   -none- NULL     
ntree              1   -none- numeric  
mtry               1   -none- numeric  
forest            14   -none- list     
y                712   factor numeric  
test               0   -none- NULL     
inbag              0   -none- NULL     
terms              3   terms  call

2. モデルの検証

In

titanic.rf.valid.pred <- array(predict(titanic.rf.model, newdata = titanic.valid.data))
accuracy(titanic.valid.data$Survived, titanic.rf.valid.pred)
# [1] 0.7932961

3. テストデータによる予測

In

titanic.rf.test.pred <- array(predict(titanic.rf.model, newdata = titanic.test.data))

titanic.rf.result <- data.frame(titanic.test.data$PassengerId, titanic.rf.test.pred)
names(titanic.rf.result) <- c('PassengerId', 'Survived')
write.csv(titanic.rf.result, 'resul_random-forest.csv', row.names=FALSE)

Public Score は 0.76794 でした。

まとめ

KaggleのTitanicをRで実装しました。
簡単なEDAとデータ前処理、モデルの構築・予測を行いました。
予測モデルの構築はロジスティック回帰とランダムフォレストのそれぞれについて行いました。

もっと良い実装があれば教えてください！