【Kaggle入門, R言語】Titanic号の乗客データを用いた生存者予測―ロジスティックモデルでスコア80%を超えることができるか挑戦してみた―

Kaggleに入門してみたので,その内容をまとめてみます.
基本的にはKernelsで優秀な参加者のコードを参考に所々アレンジしているといった内容の記事になります.

Kaggleでは入門者向けにタイタニック号の乗客データを用いた生存者予測を行うためのデータセットを提供しており、今回はそのデータセットを用いて推定モデルの構築及び提出までの一連の流れをまとめます.
今回の挑戦ではランダムフォレストのようなアンサンブル学習器を利用するのではなく,単純なロジスティックモデルでスコア80%を超えることができるかを一つの目標としてチャレンジしました.
結果を先に言うと最高スコアは0.79426(上位26%)で残念ながらあともう一息というところでした.

Titanic: Machine Learning from Disaster | Kaggle


目次

データセット及び変数の定義

予測モデルを構築するのに必要なデータとしてtrainデータ(train.csv)とtestデータ(test.csv)が用意されており,これらをダウンロードして分析を行います.乗客の生死を示すSurvived変数についてはtrainデータのみカラムに追加されており,それ以外の変数は両方のデータセットに含まれています.具体的には以下のような変数がデータに含まれています。

変数名 定義
PassengerId 乗客の識別ID
Survived 生存の場合1の値をとるダミー変数(train.csvのみ)
Pclass 乗客の社会経済的地位を代理するカテゴリカル変数(1st: Upper, 2nd: Middle, 3rd: Lower)
Name 乗客の名前
Sex 性別(maleとfemaleの文字列)
Age 年齢
SibSp 乗船している夫、妻と兄弟姉妹の数
Parch 乗船している両親と子供の数
Ticket チケット番号
Fare 運賃
Cabin 客室番号
Embarked 船場所(C: Cherbourg, Q: Queenstown, S: Southampton)


データの読み込み

# load data
> train_path <- 'C:\\Users\\admin\\kaggle\\Titanic\\train.csv'
> test_path <- 'C:\\Users\\admin\\kaggle\\Titanic\\test.csv'
> train <- read.csv(train_path, stringsAsFactors = F,na.strings=(c("NA", "")))
> test  <- read.csv(test_path, stringsAsFactors = F,na.strings=(c("NA", "")))
# データの結合
> alldata <- bind_rows(train,test) 

read.csv関数を用いてデータを読み込む際に欠損値を空白にするのでなくNAと認識させるためにオプションna.stringsを(c("NA", ""))のように指定しています.この点については以下の記事を参考にしました.

qiita.com

使用するパッケージ

使用するパッケージは以下のとおりです.

> # load packages
> library(dplyr)
> library(magrittr)
> library(ggplot2)
> library(ggthemes)
> library(car)
> library(corrplot)
> library(caret)
> library(makedummies)


データの確認

alldataの中身を確認します.

> glimpse(alldata)
Observations: 1,309
Variables: 12
$ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
$ Survived    <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0...
$ Pclass      <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3...
$ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley ...
$ Sex         <chr> "male", "female", "female", "female", "male", "male", "...
$ Age         <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 1...
$ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1...
$ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0...
$ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", ...
$ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.86...
$ Cabin       <chr> NA, "C85", NA, "C123", NA, NA, "E46", NA, NA, NA, "G6",...
$ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", ...

次に変数の基本統計量の確認をします.

# 基本統計量の確認
> summary(alldata)
 PassengerId      Survived          Pclass          Name          
 Min.   :   1   Min.   :0.0000   Min.   :1.000   Length:1309       
 1st Qu.: 328   1st Qu.:0.0000   1st Qu.:2.000   Class :character  
 Median : 655   Median :0.0000   Median :3.000   Mode  :character  
 Mean   : 655   Mean   :0.3838   Mean   :2.295                     
 3rd Qu.: 982   3rd Qu.:1.0000   3rd Qu.:3.000                     
 Max.   :1309   Max.   :1.0000   Max.   :3.000                     
                NA's   :418                                        
     Sex                 Age            SibSp            Parch      
 Length:1309        Min.   : 0.17   Min.   :0.0000   Min.   :0.000  
 Class :character   1st Qu.:21.00   1st Qu.:0.0000   1st Qu.:0.000  
 Mode  :character   Median :28.00   Median :0.0000   Median :0.000  
                    Mean   :29.88   Mean   :0.4989   Mean   :0.385  
                    3rd Qu.:39.00   3rd Qu.:1.0000   3rd Qu.:0.000  
                    Max.   :80.00   Max.   :8.0000   Max.   :9.000  
                    NA's   :263                                     
    Ticket               Fare            Cabin             Embarked        
 Length:1309        Min.   :  0.000   Length:1309        Length:1309       
 Class :character   1st Qu.:  7.896   Class :character   Class :character  
 Mode  :character   Median : 14.454   Mode  :character   Mode  :character  
                    Mean   : 33.295                                        
                    3rd Qu.: 31.275                                        
                    Max.   :512.329                                        
                    NA's   :1                                  
  • Survivedの欠損値はテストデータによるもの
  • Survivedの平均値は約0.38


Survivedとの関係を確認

予測モデルの目的変数に使用するSurvivedと説明変数に使用する変数との単一変量間の関係を確認していきます.

#'Pclass', 'Sex', 'Survived'をファクター型へ変換する
train$Pclass <- as.factor(train$Pclass)
train$Sex <- as.factor(train$Sex)
train$Survived <- factor(train$Survived,levels=c(0,1),labels=c("Died","Survived"))
> # Survived × Pclass
> SP <- table(train$Survived, train$Pclass) 
> print(SP) # クロス集計表
             1   2   3
  Died      80  97 372
  Survived 136  87 119
> round(prop.table(SP,2),digit=2)  # Pclassごとの生存(死亡)率
             1    2    3
  Died     0.37 0.53 0.76
  Survived 0.63 0.47 0.24
  • Pclassが1の値を有する乗客の生存率が0.63と高く,反対にPclassが3の値を有する乗客の生存率は0.24と低い値となっている.このことは,乗客の社会経済的地位が高いほど生存率が高いことを示唆している.
# Survived × Age
> ggplot(train, aes(Age, fill = Survived)) + 
>   geom_histogram() + 
>   theme_igray() +
>   xlab("Age") +
>   scale_fill_discrete(name = "Survived") + 
>   ggtitle("Age vs Survived")

f:id:joure:20170826120438p:plain

# SibSp
> table(train$Survived, train$SibSp)
            0   1   2   3   4   5   8
  Died     398  97  15  12  15   5   7
  Survived 210 112  13   4   3   0   0

> round(prop.table(table(train$Survived,train$SibSp),2),digit=2)          
              0    1    2    3    4    5    8
  Died     0.65 0.46 0.54 0.75 0.83 1.00 1.00
  Survived 0.35 0.54 0.46 0.25 0.17 0.00 0.00
  • 乗船している配偶者(夫、妻)と兄弟姉妹の数は1人のときに生存率が最も高くなりその後減少する
# Parch
> table(train$Survived, train$Parch) 
             0   1   2   3   4   5   6
  Died     445  53  40   2   4   4   1
  Survived 233  65  40   3   0   1   0

> round(prop.table(table(train$Survived, train$Parch),2),digit=2)
              0    1    2    3    4    5    6
  Died     0.66 0.45 0.50 0.40 1.00 0.80 1.00
  Survived 0.34 0.55 0.50 0.60 0.00 0.20 0.00

・乗船している両親(父、母)と子供の数が0人のとき生存率が低く,1人から3人では生存率が高くなり,4人以上だとまた低くなっている.

# Family sizeを表す変数Fsizeの作成:Fsize = SibSp + Parch + 1
> temp_train$Fsize <- temp_train$SibSp + temp_train$Parch + 1
# Survived × Fsize
> table(temp_train$Survived, temp_train$Fsize)
            1   2   3   4   5   6   7   8  11
  Died     374  72  43   8  12  19   8   6   7
  Survived 163  89  59  21   3   3   4   0   0
          
> round(prop.table(table(temp_train$Survived, temp_train$Fsize),2),digit=2)
             1    2    3    4    5    6    7    8   11
  Died     0.70 0.45 0.42 0.28 0.80 0.86 0.67 1.00 1.00
  Survived 0.30 0.55 0.58 0.72 0.20 0.14 0.33 0.00 0.00

・乗船している家族の人数は1人のとき生存率は0.3と低いが2人から4人までの生存率は高く,5人以上からは生存率が低くなる傾向にある.このことは,家族がいない乗客の生存率は低い傾向にあるが,家族が多すぎる乗客の生存率も低い傾向にあることを示唆している.

# Fare
ggplot(train, aes(x = Fare, fill = Survived)) +
  geom_histogram() +
  scale_x_continuous() +
  labs(x = 'Fare') +
  theme_igray()

f:id:joure:20170826184850p:plain
・運賃が高いほど生存率が上がる?

# Embarked
> table(train$Survived, train$Embarked)
            C   Q   S
  Died      75  47 427
  Survived  93  30 217
          
> round(prop.table(table(train$Survived, train$Embarked),2),digit=2)
              C    Q    S
  Died     0.45 0.61 0.66
  Survived 0.55 0.39 0.34

・乗船場所がQueenstownとSouthamptonの乗船客の死亡率が高い

説明変数の作成と欠損値の補完

# DFamsize変数の作成(Family sizeをSingle, Small, Largeにカテゴリ化)
> alldata$Famsize <- alldata$SibSp + alldata$Parch + 1 
> alldata$DFamsize[alldata$Famsize==1] <- 'Single'
> alldata$DFamsize[alldata$Famsize>1 & alldata$Famsize<5] <- 'Small'
> alldata$DFamsize[alldata$Famsize>4] <- 'Large'

# ファクター型へ変換
factor_vars <- c('Pclass','Sex','Embarked','Dfamsize','Survived')
alldata[factor_vars] <- lapply(alldata[factor_vars], function(x) as.factor(x))
# 欠損値の補完
> sapply(alldata,function(x) sum(is.na(x)))
PassengerId
0
Survived
418
Pclass
0
Name
0
Sex
0
Age
263
SibSp
0
Parch
0
Ticket
0
Fare
1
Cabin
1014
Embarked
2
Famsize
0
Dfamsize
0

・欠損値は次のとおり.Age:263, Fare:1, Cabin:1014, Embarked:2
・Cabinは使用しないので,Age, Fare, Embarkedの欠損値を補完する.

# Fareの欠損値の補完
> sum(is.na(alldata$Fare)) 
1
> which(is.na(alldata$Fare))
1044
> alldata[1044,]
PassengerId	Survived	Pclass	Name	Sex	Age	SibSp	Parch	Ticket	Fare	Cabin	Embarked	Child	Fsize	famsize	Dfamsize
1044	1044	NA	3	Storey, Mr. Thomas	male	60.5	0	0	3701	NA	NA	S	Adult	1	1	Single

# 利用できそうな変数:Pclass=3、Embarked="S"
# Pclassが3の値を有しており,乗船場所がSouthamptonの乗客の運賃の中央値で補完する
> tempdata<-alldata[which(alldata[,"Pclass"]==3 & alldata[,"Embarked"]== "S"),]
> alldata[1044,]$Fare<-median(tempdata$Fare,na.rm=T)
# Embarkedの欠損値の補完
> sum(is.na(alldata$Embarked)) 
2
> a <- which(is.na(alldata$Embarked)) # 62 and 830
> alldata[a,] # この2人の乗客の共通点はPclassが1とFareが80ということ.
PassengerId	Survived	Pclass	Name	Sex	Age	SibSp	Parch	Ticket	Fare	Cabin	Embarked	Child	Fsize	famsize	Dfamsize
62	62	1	1	Icard, Miss. Amelie	female	38	0	0	113572	80	B28	NA	Adult	1	1	Single
830	830	1	1	Stone, Mrs. George Nelson (Martha Evelyn)	female	62	0	0	113572	80	B28	NA	Adult	1	1	Single

tempdata <- alldata[-a,]
tempdata <- tempdata[which(tempdata[,"Pclass"]==1),]
table(tempdata$Embarked) # Pclassが1の値を有する乗客の乗船地:Qが圧倒的に少ない
 C   Q   S 
141   3 177 

> summary(tempdata$Fare) # Pclassが1の値を有する乗客の運賃の基本統計量:平均値が87,中央値が60
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   30.70   60.00   87.56  108.90  512.33 

> ggplot(tempdata, aes(x = Embarked, y = Fare)) +
>   geom_boxplot() +
>   geom_hline(yintercept = 80,colour = "red" ,lwd = .5)
# boxplotよりアウトサンプルにおけるPclassが1の値を有する乗客の乗船場所別の運賃の中央値(box内の横線)が80なのはCだと分かる.
# なので欠損値にCを補完する
> alldata$Embarked[a] <- "C"

f:id:joure:20170826193816p:plain

# Ageの欠損値を補完
> sum(is.na(alldata$Age)) #263 missing values for 'Age'
263

# Name変数の情報を利用して補完を行う
> alldata$Title <- gsub('(.*, )|(\\..*)', '', alldata$Name) #Name変数から呼称(Mr, Missなど)部分を抽出して新たな変数Titleとする.
# 使用している正規表現は次のとおり。
# (.*, )|(\\..*)のうち()はグループ化であり,.は任意の一文字。*は直前の一文字を0回以上続ける。
# \\.は単なる「.」文字のことであり\\..は「.の後に任意の一文字」を意味しており,さらに\\..*は「.●●」と.の後の任意の一文字を0回以上続ける,つまり,.の後の文字列を指定していることを意味する.
> table(alldata$Title)
 Capt          Col          Don         Dona           Dr     Jonkheer 
           1            4            1            1            8            1 
        Lady        Major       Master         Miss         Mlle          Mme 
           1            2           61          260            2            1 
          Mr          Mrs           Ms          Rev          Sir the Countess 
         757          197            2            8            1            1 
> officer <- c('Capt', 'Col', 'Don', 'Dr', 'Major', 'Rev')
> royalty <- c('Dona', 'Lady', 'the Countess','Sir', 'Jonkheer')

# Miss, Mrs, Royalty, Officerへ集約する
> alldata$Title[alldata$Title == 'Mlle']        <- 'Miss' 
> alldata$Title[alldata$Title == 'Ms']          <- 'Miss'
> alldata$Title[alldata$Title == 'Mme']         <- 'Mrs' 
> alldata$Title[alldata$Title %in% royalty]  <- 'Royalty'
> alldata$Title[alldata$Title %in% officer]  <- 'Officer'
> alldata$Title<-as.factor(alldata$Title)
# Titleごとの中央値でAgeの欠損値を補完する
> tapply(alldata$Age, alldata$Title,median, na.rm=TRUE)
> title.age <- aggregate(alldata$Age,by = list(alldata$Title), FUN = function(x) median(x, na.rm = T))
> title.age # Titleごとの年齢の中央値
Group.1	x
Master	4
Miss	22
Mr	29
Mrs	35
Officer	49
Royalty	39

> alldata[is.na(alldata$Age), "Age"] <- apply(alldata[is.na(alldata$Age), ] , 1, function(x) title.age[title.age[, 1]==x["Title"], 2])


相関行列

> cormat <- alldata
# ダミー化したい変数をセレクト
> dum <- cordata %>% select(Survived, Pclass,Sex, Embarked, Dfamsize, Title, group, Wom_chd)
# ダミー化しない変数をセレクト
> not_dum <- cordata %>% select(PassengerId, Name, Age, SibSp, Parch, Ticket, Fare, Cabin, famsize)
# makedummies()を使用してダミー変数を作成
> dummy_var <- makedummies(dum, basal_level = FALSE)
#作成したダミー変数の名前を修正
dummy_var %<>% dplyr::rename(Survived=res, Sex=res.1, Group=res.2, Wom_chd=res.3)
# 結合する
> cordata <- cbind(dummy_var, not_dum)      


# 相関行列の表示
> cordata_ver1 <- cordata %>% dplyr::select(Pclass_2,Pclass_3, Sex, Age, Fare, Embarked_Q,Embarked_S,
>                                           Dfamsize_Single, Dfamsize_Small, Title_Miss, Title_Mr,
>                                            Title_Mrs, Title_Officer, Title_Royalty)
> factor_vars <- c('Pclass_2','Pclass_3', 'Sex', 'Age', 'Fare', 'Embarked_Q', 'Embarked_S',
>                  'Dfamsize_Single', 'Dfamsize_Small', 'Title_Miss', 'Title_Mr', 'Title_Mrs', 
>                  'Title_Officer', 'Title_Royalty')
> cordata_ver1[factor_vars] <- lapply(cordata_ver1[factor_vars], function(x) as.numeric(x))
> cormat1 <- cor(cordata_ver1)
> corrplot(cormat1,method="circle",,numbers=T)

f:id:joure:20170830002942p:plain


データセットの分割

# alldataをもとのtrainデータとtestデータに分割
> train <- alldata[1:891,]
> test <- alldata[892:1309,]


ロジスティック回帰モデルの推定(1回目)

> # 説明変数formulaの作成
> n <- names(train)
> n <- n[-14][-15][-14:-15][-15:-17][-16:-18] # 推定に使用しないカラム名を取り除く
> formula_train <- as.formula(paste("Survived~",paste(n[!n%in%c("Survived")],collapse="+")))
> formula_train
Survived ~ Pclass_2 + Pclass_3 + Sex + Embarked_Q + Embarked_S + 
    Dfamsize_Single + Dfamsize_Small + Title_Miss + Title_Mr + 
    Title_Mrs + Title_Officer + Title_Royalty + Age + Fare
# VIFの確認
> data_vif_ver1 <- glm(formula_train,family=binomial(link='logit'),data=train)
> car::vif(data_vif_ver1) 
Pclass_2
2.0690831639864
Pclass_3
2.99470227203914
Sex
33889618.0577909
Embarked_Q
1.64928964379448
Embarked_S
1.56332970492066
Dfamsize_Single
6.98253855750122
Dfamsize_Small
5.86125762723412
Title_Miss
27259220.7276544
Title_Mr
10.1725086745876
Title_Mrs
16400975.6961539
Title_Officer
2.2800413412682
Title_Royalty
2.00000137763095
Age
1.82769059866997
Fare
1.80256557574874

VIFが10を超えている変数はSex, Title_Miss, Title_Mr, Title_Mrsである.とりあえず,Sex変数を除いてもう一度VIFを計算する.

> # 説明変数formulaの作成(Sex変数を除くバージョン)
> n <- names(train)
> n <- n[-14][-15][-14:-15][-15:-17][-16:-18][-4] # 不要なカラム名を取り除く
> formula_train <- as.formula(paste("Survived~",paste(n[!n%in%c("Survived")],collapse="+")))
> formula_train
Survived ~ Pclass_2 + Pclass_3 + Embarked_Q + Embarked_S + Dfamsize_Single + 
    Dfamsize_Small + Title_Miss + Title_Mr + Title_Mrs + Title_Officer + 
    Title_Royalty + Age + Fare
# VIFの確認
> data_vif_ver1 <- glm(formula_train,family=binomial(link='logit'),data=train)
> car::vif(data_vif_ver1) 
Pclass_2
2.01492616765143
Pclass_3
2.89688819262399
Embarked_Q
1.59598916928248
Embarked_S
1.52675468475323
Dfamsize_Single
7.04456490757105
Dfamsize_Small
5.80350635227861
Title_Miss
6.28385241950577
Title_Mr
9.39757461218912
Title_Mrs
4.80463078750706
Title_Officer
2.26189845814524
Title_Royalty
1.2724980598487
Age
1.89565761614302
Fare
1.71653693562228

Title_Mrが限りなく10に近く怪しいがすべての変数が10以下なので,この説明変数を使用してモデルを推定する(VIFが10以下はあくまで目安?).

# ロジスティック回帰モデルの推定(1回目)
> train$Survived <- as.factor(train$Survived) #Survived変数がnumeric型になっていたのでfactor型に変換する
> acc_data_ver1 <- train(
>  data = train,
>  formula_train,
>  method = "glm",
>  family = binomial())
> summary(acc_data_ver1)
Call:
NULL

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6752  -0.5159  -0.3964   0.5543   2.4282  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.979935   0.665638   1.472 0.140974    
Pclass_2        -1.110466   0.324996  -3.417 0.000633 ***
Pclass_3        -2.072061   0.322692  -6.421 1.35e-10 ***
Embarked_Q      -0.011785   0.399375  -0.030 0.976459    
Embarked_S      -0.314052   0.254656  -1.233 0.217485    
Dfamsize_Single  3.110923   0.512692   6.068 1.30e-09 ***
Dfamsize_Small   2.749036   0.481291   5.712 1.12e-08 ***
Title_Miss      -0.608128   0.545230  -1.115 0.264696    
Title_Mr        -3.525318   0.582137  -6.056 1.40e-09 ***
Title_Mrs        0.127135   0.616027   0.206 0.836495    
Title_Officer   -3.447805   0.827004  -4.169 3.06e-05 ***
Title_Royalty   -1.794382   1.313209  -1.366 0.171810    
Age             -0.024596   0.009658  -2.547 0.010871 *  
Fare             0.004111   0.002681   1.533 0.125187    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1186.66  on 890  degrees of freedom
Residual deviance:  716.11  on 877  degrees of freedom
AIC: 744.11

Number of Fisher Scoring iterations: 5
# ブートストラップ法による予測精度の検証(caret::train()を使用)
> acc_data_ver1
Generalized Linear Model 

891 samples
 13 predictor
  2 classes: '0', '1' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 891, 891, 891, 891, 891, 891, ... 
Resampling results:

  Accuracy   Kappa    
  0.8246645  0.6273165

Accuracyが約0.825である.

# 提出
> logit_ver1 <- glm(formula_train,family=binomial(link='logit'),data=train)
> 
> fitted_ver1 <- predict(logit_ver1, newdata=test,type='response')
> fitted_ver1 <- ifelse(fitted_ver1 > 0.5,1,0)
> solution_ver1 <- data.frame(PassengerId = test$PassengerId, Survived = fitted_ver1)
> write.csv(solution_ver1, file = 'submit_ver1.csv', row.names = F)
> # # Your submission scored 0.77990

提出した結果,スコア(Accuracy)は0.77990であった.


説明変数の追加

# 女性もしくは18歳以下ならばYes, それ以外はNoの値とるWom_chd変数の作成
> alldata$Wom_chd <- "No"
> alldata$Wom_chd[which(alldata$Sex == "female" | alldata$Age < 18)] = "Yes"
> alldata$Wom_chd <- as.factor(alldata$Wom_chd)
# Ticketを利用した変数の作成
> head(table(alldata$Ticket),20)
110152 110413 110465 110469 110489 110564 110813 111163 111240 111320 111361 
     3      3      2      1      1      1      2      1      1      1      2 
111369 111426 111427 111428 112050 112051 112052 112053 112058 
     1      1      1      1      1      1      1      1      2 

・同じチケット番号ならば行動を共にしていたと考えられる.
・チケット番号ごとに集計して,その数が2より大きければYes, それ以外ならばNoの値をとる変数を作成したい

# Ticketk変数を番号ごとに集計した結果を新たな変数groupTKTとする.
> groupTKT <- alldata %>%
>         group_by(Ticket) %>%
>         summarise(N = n()) %>%
>         filter(N > 2) %>%
>         arrange(desc(N))
> head(groupTKT, 5)
Ticket	N
CA. 2343 11
1601	8
CA 2144	8
3101295	7
347077	7

# groupTKTの値が2より大きければYes,それ以外ならNoの値をとるGroup変数の作成
> alldata$Group = "No"
> alldata$Group[which(alldata$Ticket %in% groupTKT$Ticket)] = "Yes"
> alldata$Group <- as.factor(alldata$Group)



ロジスティック回帰モデルの推定(2回目)

# 再度データを分割
train <- alldata[1:891,]
test <- alldata[892:1309,]
> # 説明変数formulaの作成
> n <- names(train)
> n <- n[-7:-8][-14:-15][-15:-17][-16][-4][-7] # 不要なカラム名を取り除く
> n
> formula_train_ver2 <- as.formula(paste("Survived~",paste(n[!n%in%c("Survived")],collapse="+")))
> formula_train_ver2
Survived ~ Pclass_2 + Pclass_3 + Embarked_Q + Embarked_S + Title_Miss + 
    Title_Mrs + Title_Officer + Title_Royalty + Group + Wom_chd + 
    Age + Fare + famsize
# VIFの確認
> data_vif_ver2 <- glm(formula_train_ver2,family=binomial(link='logit'),data=train)
> car::vif(data_vif_ver2) 
Pclass_2
2.10544001506684
Pclass_3
2.9097090314934
Embarked_Q
1.55142003971765
Embarked_S
1.48615277528479
Title_Miss
3.0353767265435
Title_Mrs
2.87137916357893
Title_Officer
1.0837337860472
Title_Royalty
1.03124236723129
Group
2.18353274788386
Wom_chd
4.52577456143758
Age
2.03649236513551
Fare
1.79003630590504
famsize
2.19602987015212
# ロジスティック回帰モデルの推定(2回目)
> train$Survived <- as.factor(train$Survived)
> 
> train_acc_ver2 <- train(
>  data = train,
>  formula_train_ver2,
>  method = "glm",
>  family = binomial())
> summary(train_acc_ver2)
Call:
NULL

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4420  -0.6190  -0.3769   0.5642   2.5076  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    1.4962097  0.4964748   3.014 0.002581 ** 
Pclass_2      -1.1683218  0.3197397  -3.654 0.000258 ***
Pclass_3      -2.3262555  0.3146478  -7.393 1.43e-13 ***
Embarked_Q     0.0335398  0.3878851   0.086 0.931094    
Embarked_S    -0.4141321  0.2466684  -1.679 0.093171 .  
Title_Miss     0.8721582  0.3677343   2.372 0.017706 *  
Title_Mrs      1.7756888  0.4532482   3.918 8.94e-05 ***
Title_Officer -0.0172768  0.5949102  -0.029 0.976832    
Title_Royalty  1.0444122  1.2946644   0.807 0.419836    
Group          0.7245855  0.3108476   2.331 0.019753 *  
Wom_chd        1.9237960  0.3923751   4.903 9.44e-07 ***
Age           -0.0307301  0.0096504  -3.184 0.001451 ** 
Fare           0.0006951  0.0026086   0.266 0.789879    
famsize       -0.4784253  0.0942377  -5.077 3.84e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1186.66  on 890  degrees of freedom
Residual deviance:  743.42  on 877  degrees of freedom
AIC: 771.42

Number of Fisher Scoring iterations: 5
# ブートストラップ法による予測精度の検証(caret::train()を使用)
> train_acc_ver2
Generalized Linear Model 

891 samples
 13 predictor
  2 classes: '0', '1' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 891, 891, 891, 891, 891, 891, ... 
Resampling results:

  Accuracy   Kappa   
  0.8139758  0.602856
# 提出2
> logit_ver2 <- glm(formula_train_ver2,family=binomial(link='logit'),data=train)
> 
> fitted_ver2 <- predict(logit_ver2, newdata=test,type='response')
> fitted_ver2 <- ifelse(fitted_ver2 > 0.5,1,0)
> solution_ver2 <- data.frame(PassengerId = test$PassengerId, Survived = fitted_ver2)
> write.csv(solution_ver2, file = 'submit_ver2.csv', row.names = F)
# Your submission scored 0.79426

f:id:joure:20170827001447p:plain


おわりに

Kaggleでは書いたコードをKernelsで公開している人がたくさんいます.そのコードを見るだけでとても勉強になるなーというのが一番の感想です.あとロジスティックモデルを使用した人のKernelsをみていましたが,VIFを確認するなど多重共線性をケアしている人があまりいないんだなーというのも印象に残りました*1


*1:アンサンブル学習でケアしないのはわかりますが・・・