最近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とデータ前処理、モデルの構築・予測を行いました。
- 予測モデルの構築はロジスティック回帰とランダムフォレストのそれぞれについて行いました。
もっと良い実装があれば教えてください！