KaggleのTitanicにbase Rと glmnetで挑戦してみました。
quartoでhtmlを作って、下に貼り付けました。
何かのご参考になれば幸いです。
Kaggle Titanic with base R and glmnet
はじめに
この Quarto ドキュメントは、Kaggle の Titanic の competition に
base R と glmnet で挑む記録を記したものです。
ライブラリーの読み込み
はじめに glmnet パッケージの読み込みをします。
library(glmnet)Loading required package: MatrixLoaded glmnet 4.1-10CSVファイルの読み込み
トレーニング用データとテスト用のデータのCSVファイルを読み込みます。
read.csv()関数を使います。
df_train <- read.csv("C:/R_folder/kaggle_titanic_train.csv")
df_test <- read.csv("C:/R_folder/kaggle_titanic_test.csv")正しく読み込めたか、str()関数で確認してみます。
str(df_train)'data.frame': 891 obs. of 12 variables:
$ PassengerId: int 1 2 3 4 5 6 7 8 9 10 ...
$ Survived : int 0 1 1 1 0 0 0 0 1 1 ...
$ Pclass : int 3 1 3 1 3 3 1 3 3 2 ...
$ Name : chr "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
$ Sex : chr "male" "female" "female" "female" ...
$ Age : num 22 38 26 35 35 NA 54 2 27 14 ...
$ SibSp : int 1 1 0 1 0 0 0 3 0 1 ...
$ Parch : int 0 0 0 0 0 0 0 1 2 0 ...
$ Ticket : chr "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...
$ Fare : num 7.25 71.28 7.92 53.1 8.05 ...
$ Cabin : chr "" "C85" "" "C123" ...
$ Embarked : chr "S" "C" "S" "S" ...str(df_test)'data.frame': 418 obs. of 11 variables:
$ PassengerId: int 892 893 894 895 896 897 898 899 900 901 ...
$ Pclass : int 3 3 2 3 3 3 3 2 3 3 ...
$ Name : chr "Kelly, Mr. James" "Wilkes, Mrs. James (Ellen Needs)" "Myles, Mr. Thomas Francis" "Wirz, Mr. Albert" ...
$ Sex : chr "male" "female" "male" "male" ...
$ Age : num 34.5 47 62 27 22 14 30 26 18 21 ...
$ SibSp : int 0 1 0 0 1 0 0 1 0 2 ...
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ Ticket : chr "330911" "363272" "240276" "315154" ...
$ Fare : num 7.83 7 9.69 8.66 12.29 ...
$ Cabin : chr "" "" "" "" ...
$ Embarked : chr "Q" "S" "Q" "S" ...トレーニング用データとテスト用のデータの結合
データの前処理をするためには、
トレーニング用データとテスト用のデータを結合すると便利です。
テスト用データには、Survivedが無いので、Survivedを加えます
df_test$Survived <- NA2つのデータフレームを縦方向につなげるときは、rbind()関数を使いますが
列の順番が同じでないとだめなので、工夫が必要です。
まず、df_trainの列の順番を保存します。
col_order <- names(df_train)
col_order [1] "PassengerId" "Survived" "Pclass" "Name" "Sex"
[6] "Age" "SibSp" "Parch" "Ticket" "Fare"
[11] "Cabin" "Embarked" そうしたら、保存したcol_orderを使って、df_testの列の順番を変更します。
変更前の列の順番はこうなります。
str(df_test)'data.frame': 418 obs. of 12 variables:
$ PassengerId: int 892 893 894 895 896 897 898 899 900 901 ...
$ Pclass : int 3 3 2 3 3 3 3 2 3 3 ...
$ Name : chr "Kelly, Mr. James" "Wilkes, Mrs. James (Ellen Needs)" "Myles, Mr. Thomas Francis" "Wirz, Mr. Albert" ...
$ Sex : chr "male" "female" "male" "male" ...
$ Age : num 34.5 47 62 27 22 14 30 26 18 21 ...
$ SibSp : int 0 1 0 0 1 0 0 1 0 2 ...
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ Ticket : chr "330911" "363272" "240276" "315154" ...
$ Fare : num 7.83 7 9.69 8.66 12.29 ...
$ Cabin : chr "" "" "" "" ...
$ Embarked : chr "Q" "S" "Q" "S" ...
$ Survived : logi NA NA NA NA NA NA ...Survivedが一番最後にあることがわかります。
これをdf_trainと同じ順番にします。
df_test <- df_test[ , col_order]
str(df_test)'data.frame': 418 obs. of 12 variables:
$ PassengerId: int 892 893 894 895 896 897 898 899 900 901 ...
$ Survived : logi NA NA NA NA NA NA ...
$ Pclass : int 3 3 2 3 3 3 3 2 3 3 ...
$ Name : chr "Kelly, Mr. James" "Wilkes, Mrs. James (Ellen Needs)" "Myles, Mr. Thomas Francis" "Wirz, Mr. Albert" ...
$ Sex : chr "male" "female" "male" "male" ...
$ Age : num 34.5 47 62 27 22 14 30 26 18 21 ...
$ SibSp : int 0 1 0 0 1 0 0 1 0 2 ...
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ Ticket : chr "330911" "363272" "240276" "315154" ...
$ Fare : num 7.83 7 9.69 8.66 12.29 ...
$ Cabin : chr "" "" "" "" ...
$ Embarked : chr "Q" "S" "Q" "S" ...SurvivedがPassengerIdとPclassの間にあることがわかります。 そうしたら、rbind()関数でdf_trainとdf_testを結合します。
df_all <- rbind(df_train, df_test)
head(df_all) PassengerId Survived Pclass
1 1 0 3
2 2 1 1
3 3 1 3
4 4 1 1
5 5 0 3
6 6 0 3
Name Sex Age SibSp Parch
1 Braund, Mr. Owen Harris male 22 1 0
2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38 1 0
3 Heikkinen, Miss. Laina female 26 0 0
4 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1 0
5 Allen, Mr. William Henry male 35 0 0
6 Moran, Mr. James male NA 0 0
Ticket Fare Cabin Embarked
1 A/5 21171 7.2500 S
2 PC 17599 71.2833 C85 C
3 STON/O2. 3101282 7.9250 S
4 113803 53.1000 C123 S
5 373450 8.0500 S
6 330877 8.4583 Qtail(df_all) PassengerId Survived Pclass Name Sex Age
1304 1304 NA 3 Henriksson, Miss. Jenny Lovisa female 28.0
1305 1305 NA 3 Spector, Mr. Woolf male NA
1306 1306 NA 1 Oliva y Ocana, Dona. Fermina female 39.0
1307 1307 NA 3 Saether, Mr. Simon Sivertsen male 38.5
1308 1308 NA 3 Ware, Mr. Frederick male NA
1309 1309 NA 3 Peter, Master. Michael J male NA
SibSp Parch Ticket Fare Cabin Embarked
1304 0 0 347086 7.7750 S
1305 0 0 A.5. 3236 8.0500 S
1306 0 0 PC 17758 108.9000 C105 C
1307 0 0 SOTON/O.Q. 3101262 7.2500 S
1308 0 0 359309 8.0500 S
1309 1 1 2668 22.3583 C結合できました。
EDAと前処理
それでは、EDA (Exploratory Data Analysis) と前処理をやっていきます。
Survived
Survivedは0なら死亡、1なら生存です。
生存率を計算しましょう。
mean(df_all$Survived, na.rm = TRUE)[1] 0.3838384生存率は、38.4%でした。
Pclass
次は、Pclassです。これは、1、2、3に分かれています。
まず、NAが無いかみてみます。
sum(is.na(df_all$Pclass))[1] 0PclassにはNAは無いようです。
Pclassは、1, 2, 3と数値データですが、これは本当は、
1等客室、2等客室、3等客室のことなのでファクター型に変更します。
df_all$Pclass <- as.factor(df_all$Pclass)
summary(df_all$Pclass) 1 2 3
323 277 709 3等客室の乗客が709人で一番多いことがわかります。
客室の等級によって生存率に差があるでしょうか？
バーチャートで視覚化してみましょう。
colors <- c("red", "green", "yellow")
barplot(table(df_all$Pclass, df_all$Survived),
col = colors, legend = FALSE,
main = "1 = Survived / 0 = Not-Survived")
legend("topright", legend = c("Pclass 1", "Pclass 2", "Pclass 3"),
fill = colors, border = "black")
0, 死亡者のほうが、1, 生存者より多いこと、
Pclass 1の人は生存者のほうが多いことなどがわかります。
PclassがSurvivedに有意に関連しているかどうか、chisq.test()関数で検定します。
chisq.test( table(df_all$Pclass, df_all$Survived) )
Pearson's Chi-squared test
data: table(df_all$Pclass, df_all$Survived)
X-squared = 102.89, df = 2, p-value < 2.2e-16p-value < 2.2e-16 とほとんど0ですので、
PclassはSurvivedに関連していることが確認できました。
Name
Nameは文字通り乗客の名前です。
いろいろ調べると、このNameの中からMr.やMiss.やMaster.などの肩書を
取得するのがTitancのキーポイントらしいです。
とりあえず、数行表示してみます。
head(df_all[ , "Name"], 20) [1] "Braund, Mr. Owen Harris"
[2] "Cumings, Mrs. John Bradley (Florence Briggs Thayer)"
[3] "Heikkinen, Miss. Laina"
[4] "Futrelle, Mrs. Jacques Heath (Lily May Peel)"
[5] "Allen, Mr. William Henry"
[6] "Moran, Mr. James"
[7] "McCarthy, Mr. Timothy J"
[8] "Palsson, Master. Gosta Leonard"
[9] "Johnson, Mrs. Oscar W (Elisabeth Vilhelmina Berg)"
[10] "Nasser, Mrs. Nicholas (Adele Achem)"
[11] "Sandstrom, Miss. Marguerite Rut"
[12] "Bonnell, Miss. Elizabeth"
[13] "Saundercock, Mr. William Henry"
[14] "Andersson, Mr. Anders Johan"
[15] "Vestrom, Miss. Hulda Amanda Adolfina"
[16] "Hewlett, Mrs. (Mary D Kingcome) "
[17] "Rice, Master. Eugene"
[18] "Williams, Mr. Charles Eugene"
[19] "Vander Planke, Mrs. Julius (Emelia Maria Vandemoortele)"
[20] "Masselmani, Mrs. Fatima" 観察すると、はじめのカンマの後にスペースがあって、
次にMr.やMrs.などの肩書が来るようです。
この肩書を取り出すにはどうしたらいいでしょうか？
はじめのカンマの位置と2番目の空白(スペース)の位置がわかれば
なんとかなりそうです。
例えば、
1234, Mr. ab,c
という文字列があるとします。
そうすると、はじめのカンマの位置は、5です。
2番目の空白の位置は、10です。
すると、Mr. は、(5+2)から(10-1)までになります。
はじめのカンマの位置を取得する関数は、regexpr()関数を使います 例:
full_name <- "1234, Mr. ab,c"
comma <- regexpr(",", full_name)[1]
comma[1] 52番目の空白の位置を取得するのは、gregexpr()関数を使います。
regexpr()関数は、はじめに出てくる位置しか返しませんが、
gregexpr()関数は、すべての位置を返します。
例:
full_name <- "1234, Mr. ab,c"
spaces <- gregexpr(" ", full_name)[[1]]
space2 <- spaces[2]
space2[1] 10こうして、前後の位置がわかりましたので、Mr.が取り出せます。
これは、substr()関数で取り出せます。
今、commaには、5が入っていて、space2には10が入っています。
例:
full_name <- "1234, Mr. ab,c"
comma <- regexpr(",", full_name)[1]
spaces <- gregexpr(" ", full_name)[[1]]
space2 <- spaces[2]
substr(full_name, comma + 2, space2 - 1)[1] "Mr."上手く取り出せました。 いままでのことを利用して、肩書を取り出す関数式をつくりましょう。
get_katagaki <- function(string) {
comma <- regexpr(",", string)[1] # カンマの位置
spaces <- gregexpr(" ", string)[[1]] # 空白の位置(複数ある)
space2 <- spaces[2] # 2番目の空白の位置
substr(string, comma + 2, space2 - 1) #肩書を取り出す
}このget_katagaki関数がうまく動くか試してみます。
get_katagaki(df_all$Name[1])[1] "Mr."get_katagaki(df_all$Name[2])[1] "Mrs."get_katagaki(df_all$Name[3])[1] "Miss."get_katagaki(df_all$Name[4])[1] "Mrs."get_katagaki(df_all$Name[5])[1] "Mr."うまくいきましたね！
このget_katagaki()関数をdf_all$Nameに使うためには、sapply()関数を使います。
このようにします。
df_all$title <- sapply(df_all$Name, get_katagaki)
head(df_all[ , "title"])[1] "Mr." "Mrs." "Miss." "Mrs." "Mr." "Mr." titleの頻度をみてみます。
rev(sort(table(df_all$title)))
Mr. Miss. Mrs. Master. Rev. Dr. Col.
736 256 191 59 36 8 8 4
Ms. Mlle. Major. the Mme. Jonkheer. Don. Capt.
2 2 2 1 1 1 1 1 Mr.が736人、Miss.が256人、Mrs.が191人、Master.が59人で、空白が36人です。
その他は8人以下ですね。
titleは、Mr. Miss. Mrs. Master. その他は othersとして5つにまとめましょう。
# 4つのうちのどれかなら、TRUE
good_katagaki <- df_all$title %in% c("Mr.", "Miss.", "Mrs.", "Master.")
# good_katagakiでないものをothersにする
df_all$title[!good_katagaki] <- "others"
# うまくできたか確認
table(df_all$title)
Master. Miss. Mr. Mrs. others
59 256 736 191 67 うまくできました。
では、titleとSurvivedの関係を見てみます。
table(df_all$title, df_all$Survived)
0 1
Master. 17 23
Miss. 53 126
Mr. 423 79
Mrs. 24 97
others 32 17Mr.は死亡者が多くて、Miss.やMrs.は生存者が多いですね。
比率のテーブルにしてみます。prop.table()関数を使います。
prop.table(table(df_all$title, df_all$Survived), margin = 1)
0 1
Master. 0.4250000 0.5750000
Miss. 0.2960894 0.7039106
Mr. 0.8426295 0.1573705
Mrs. 0.1983471 0.8016529
others 0.6530612 0.3469388margin = 1 を使うと、横方法の合計が1.0になる比率のテーブルになります。
Master.は58%, Miss.は70%, Mr.は16%, Mrs.は80%, othersは35%の生存率です。
グラフも描いてみます。
colors <- c("yellow", "red", "green", "pink", "gray")
barplot(prop.table(table(df_all$title, df_all$Survived)),
legend = FALSE, col = colors,
main = "Title and Survived (1 = Survuved)")
legend("topright", legend = c("Master.", "Miss.", "Mr.", "Mrs.", "others"),
fill = colors, border = "black")
Mr.は多くが死んでしまったことがわかります。
Nameの情報をtitleとして5つの肩書に集約しました。
titleはファクター型にして、
Nameは、もう必要が無いので、df_allから削除します。
df_all$title <- as.factor(df_all$title)
df_all$Name <- NULLSex
次はSexを調べます。まず、NAがあるかどうかを確認します。
sum(is.na(df_all$Sex))[1] 0結果は0でしたので、ファクター型にします。
df_all$Sex <- as.factor(df_all$Sex)SexとSurvivedの関係をテーブルにしてみます。
table(df_all$Sex, df_all$Survived)
0 1
female 81 233
male 468 109女性は生存が多く、男性は生存は少ないですよね。
これは、titleでもMiss., Mrs.が生存が多くて、Mr.が生存が少なかったので
予想された結果です。
SurvivedとSexの関係もグラフにします。
barplot(table(df_all$Sex, df_all$Survived),
legend = TRUE,
main = "Survuved(= 1) and Sex")
femaleは多くが生存しています。
次に気になるのは、titleとSexです。
テーブルにしてみます。
table(df_all$title, df_all$Sex)
female male
Master. 0 59
Miss. 256 0
Mr. 0 736
Mrs. 191 0
others 19 48others だけ男女ともにあります。Master.とMr.は男性だけ、
Miss.とMrs.は女性だけです。
Age
Ageを調べます。まず、NAがあるかどうかを確認します。
sum(is.na(df_all$Age))[1] 263Ageは、263もNAがあります。このNAを何かの値で置換しないといけないです。 その前に、AgeがNAだったら1で、そうでなければ0のダミー変数を作っておきます。
df_all$agena <- ifelse(is.na(df_all$Age), 1, 0)
table(df_all$agena)
0 1
1046 263 1の値を持つ観測データが263個です。
AgeはPclassとtitleで中央値や平均値が違うと思います。
確認してみましょう。
# 中央値
tapply(df_all$Age, list(df_all$Pclass, df_all$title),
function(x) median(x, na.rm = TRUE)) Master. Miss. Mr. Mrs. others
1 6 30 42 45 45
2 2 20 30 31 30
3 6 18 26 31 30# 平均値
tapply(df_all$Age, list(df_all$Pclass, df_all$title),
function(x) round(mean(x, na.rm = TRUE),2)) Master. Miss. Mr. Mrs. others
1 6.98 30.34 41.58 43.61 42.18
2 2.76 20.38 32.40 33.70 35.65
3 5.94 17.43 28.22 32.04 28.98# 観測数
tapply(df_all$Age, list(df_all$Pclass, df_all$title),
function(x) sum(!is.na(x))) Master. Miss. Mr. Mrs. others
1 5 59 129 66 25
2 11 46 134 53 17
3 36 101 299 45 20Master. は子どもなんですね。知らなかったです。
AgeのNAをどうするか、今回はPclass(3通り) X title(5通り) = 15通りの
組み合わせで、NAは中央値に置き換えてみたいと思います。
まず、NAに置き換える自作関数を作ります。
replace_na <- function(x) { # xはベクトル
med <- median(x, na.rm = TRUE) # xの中央値を計算
x[is.na(x)] <- med # NAのところをmed(中央値)に置換
return(x) # x を返す
}replace_na()が上手く動くか試します。
test <- c(1, 2, NA, 4, 5)
test[1] 1 2 NA 4 5replace_na(test)[1] 1 2 3 4 5NAが中央値の3に置き換わっています。
次は、Pclass X title = 15とおりにこのreplace_na()関数を適合することを
考えます。
これは、ave()関数を使います。
ave()関数は、ave(x, g1, g2, FUN)のように使い、
x を g1 X g2 のグループにして、それぞれにFUNという関数を適用する。
というものです。
試しにこんなデータフレームを作ってみます。
df_tameshi <- data.frame(
x = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
g1 = c("a", "a", "a", "b", "b", "b", "c", "c", "c", "c"),
g2 = c("d", "e", "d", "e", "d", "e", "d", "e", "d", "e")
)
df_tameshi x g1 g2
1 1 a d
2 2 a e
3 3 a d
4 4 b e
5 5 b d
6 6 b e
7 7 c d
8 8 c e
9 9 c d
10 10 c eこのデータフレームで、g1(3とおり) X g2(2とおり) = 6つのグループでの
平均値をみてみます。
ave(df_tameshi$x, df_tameshi$g1, df_tameshi$g2, FUN = mean) [1] 2 2 2 5 5 5 8 9 8 9この結果が正しいか検証しやすいように、tapply()関数で
それぞれのグループの平均値を算出してみます。
tapply(df_tameshi$x, list(df_tameshi$g1, df_tameshi$g2), mean) d e
a 2 2
b 5 5
c 8 9平均値は、ad = 2, ae = 2, bd = 5, be = 5, cd = 8, ce = 9 です。
xは、ad, ae, ad, be, bd, be, cd, ce, cd, ce の順番になっています。
平均値に置き換えると、
2, 2, 2, 5, 5, 5, 8, 9, 8, 9ですので、ave()関数の結果と同じです。
ave()関数の働きを確認しましたので、実際にAgeのNAを
それぞれの中央値に置き換えます。
df_all$Age <- ave(df_all$Age,
df_all$Pclass, df_all$title,
FUN = replace_na)
summary(df_all$Age) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.17 21.00 26.00 29.30 36.50 80.00 replace_na()関数は、NAを中央値に置き換える関数なので、
summary()関数の結果をみると、NAがなくなっていることがわかります。
これで、Ageの前処理は終わりました。
SurvivalとAgeの関係を見てみます。
hist()関数でヒストグラムで、Survived == 1 と Survived == 0 で
わけて描いてみます。
par(mfrow = c(2, 1))
kugiri <- seq(0, 80, by = 10)
hist(df_all$Age[df_all$Survived == 1], main = "生存", breaks = kugiri,
xlab = "年齢")
hist(df_all$Age[df_all$Survived == 0], main = "死亡", breaks = kugiri,
xlab = "年齢")
par(mfrow = c(1, 1))生存者のほうが年齢の幅が広いように見えます。
SibSpとParch
次は、SibSpとParchです。
これは、合計して、familyとします。
df_all$family = df_all$SibSp + df_all$Parch
summary(df_all$family) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.0000 0.0000 0.8839 1.0000 10.0000 familyは一緒に乗船している家族の人数です。0は一人でタイタニックに
乗った人です。
最大が10ですから、11人の家族でのった人もいますね。
SibSpとParchはもう必要ないので、削除します。
df_all$SibSp <- NULL
df_all$Parch <- NULLfamilyの分布をみてみます。 table()関数を使います。
table(df_all$family)
0 1 2 3 4 5 6 7 10
790 235 159 43 22 25 16 8 11 0が790人、1が235人、2が159人でこの3つで大半を占めていますね。
familyとSurvivedの関係をみてみます。
familyごとにSurvivedの平均値を計算して、生存率をみます。
tapply(df_all$Survived, df_all$family,
FUN = function(x) round(mean(x, na.rm = TRUE), 2)) 0 1 2 3 4 5 6 7 10
0.30 0.55 0.58 0.72 0.20 0.14 0.33 0.00 0.00 0は30%, 1は55%, 2は58%, 3は72%です。
グラフも描いてみます。
# グラフのための表を生成
object <- tapply(df_all$Survived, df_all$family,
FUN = function(x) round(mean(x, na.rm = TRUE), 2))
# 高さに比例した色（生存率が高いほど明るい）
cols <- gray(object)
# グラフ描画
barplot(object, col = cols, main = "Survival Rate by Family")
こうしてグラフを見ると、familyは、0, 1~3, 4~10 の3つのグループに
カテゴリー化したほうがよさそうです。
どうやりましょうか？
cut()関数を使います。 cut()関数は数値型のベクトルをファクター型に関数です。 cut(ベクトル, breaks = 区切り, labels = ファクターの水準) という構造です。
breaks = c(-1, 0, 3, Inf) とすると、
(-1, 0]のグループ(-1より大きくて0以下)
(0, 3]のグループ(0より大きくて3以下)
(3, Inf]のグループ(3より大きくてInf以下)
に分かれます。
df_all$family <- cut(df_all$family,
breaks = c(-1, 0, 3, Inf),
labels = c("alone", "small", "large"))
summary(df_all$family)alone small large
790 437 82 できました。
Ticket
Ticketをみてみましょう。これは文字列ですね。とりあえず、20個くらいを
みてみます。
head(df_all[ , c("PassengerId", "Ticket")], 20) PassengerId Ticket
1 1 A/5 21171
2 2 PC 17599
3 3 STON/O2. 3101282
4 4 113803
5 5 373450
6 6 330877
7 7 17463
8 8 349909
9 9 347742
10 10 237736
11 11 PP 9549
12 12 113783
13 13 A/5. 2151
14 14 347082
15 15 350406
16 16 248706
17 17 382652
18 18 244373
19 19 345763
20 20 2649数字だけのものもあれば、アルファベットも混じったものもありますね。
NAはあるでしょうか？
sum(is.na(df_all$Ticlet))[1] 0NAは無いです。
どうしましょうか？ はじめの1文字だけにしてみましょう。
substr()関数を使います。
df_all$Ticket <- substr(df_all$Ticket, 1, 1)
head(df_all$Ticket, 20) [1] "A" "P" "S" "1" "3" "3" "1" "3" "3" "2" "P" "1" "A" "3" "3" "2" "3" "2" "3"
[20] "2"できました。
頻度と、Survivedとのクロス表を作成しましょう。
# 頻度
table(df_all$Ticket)
1 2 3 4 5 6 7 8 9 A C F L P S W
210 278 429 11 3 9 13 2 2 42 77 13 5 98 98 19 # 数のクロス表
table(df_all$Survived, df_all$Ticket)
1 2 3 4 5 6 7 8 9 A C F L P S W
0 54 98 229 8 3 5 8 2 0 27 31 3 3 23 44 11
1 92 85 72 2 0 1 1 0 1 2 16 4 1 42 21 2# 比率のクロス表
round(prop.table(table(df_all$Survived, df_all$Ticket), margin = 2), 2)
1 2 3 4 5 6 7 8 9 A C F L P S
0 0.37 0.54 0.76 0.80 1.00 0.83 0.89 1.00 0.00 0.93 0.66 0.43 0.75 0.35 0.68
1 0.63 0.46 0.24 0.20 0.00 0.17 0.11 0.00 1.00 0.07 0.34 0.57 0.25 0.65 0.32
W
0 0.85
1 0.15ある程度数が多くて、生存率が高いのは、1(63%), 2(46%), P(65%)ですね。
1, 2, P をgt(good ticketの意味)、それ以外をbt(bad ticket)にしてみます。
df_all$Ticket[df_all$Ticket %in% c("1", "2", "P")] <- "gt"
df_all$Ticket[df_all$Ticket != "gt"] <- "bt"
table(df_all$Ticket)
bt gt
723 586 btが723人、gtが586人です。1(210人) + 2(278人) + P(98人) = 586人なので正しく
分類されました。
ファクター型にしておきます。
df_all$Ticket <- as.factor(df_all$Ticket)Ticket別のSurvivedをグラフにしてみます。
colors <- c("blue", "red")
barplot(prop.table(table(df_all$Survived, df_all$Ticket)),
legend = FALSE, col = colors,
main = "Survived by Ticket")
legend("topright", legend = c("Not Survived", "Survived"),
fill = colors)
gtのチケットを持った人は生存者が多いことがわかります。
Fare
Fareをみてみましょう。これは乗船料ですね。
NAがあるかどうかを確認します。
is.na()関数で、NAならばTRUEを返して、sum()関数でTRUEの数を合計します。
Rでは、TRUEは1、FALSEは0としてsum()関数は処理します。
sum(is.na(df_all$Fare))[1] 11つNAがありました。これは、中央値に置き換えます。
df_all$Fare[is.na(df_all$Fare)] <- median(df_all$Fare, na.rm = TRUE)Fareの分布をヒストグラムにしてみます。
hist(df_all$Fare)
あらら、これは完全に右側の裾野が広い対数正規分布のような分布ですね。
summary()関数で最小値、平均値、中央値、最大値などと確認します。
summary(df_all$Fare) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 7.896 14.454 33.281 31.275 512.329 最小値が0なので、そのまま対数変換はできないですね。
1を足してから対数変換してみようかと思います。
ただ、その前にヒストグラムで、対数変換したらどのような分布になるかを
確認します。
hist(log(df_all$Fare + 1, base = 10))
う～ん、ちょっと気に入らないです。
平方根ではどうでしょうか
hist(sqrt(df_all$Fare))
3乗根はどうでしょうか?
hist(df_all$Fare^(1 / 3))
悩みますね。。。
そもそも、なんでこんなに乗船料がいろいろなんでしょうね。
新幹線なら自由席、指定席、グリーン席の3種類だし、
飛行機ならエコノミー、プレミアムエコノミー、ビジネス、ファーストクラスの
4種類くらいですよね。。。
ここは思い切って、4分位にわけてカテゴリー化してしまいましょう。
まず、4分位の値を求めます。quantile()関数を使います。
yonbuni <- quantile(df_all$Fare,
probs = c(0, 0.25, 0.5, 0.75, 1))
yonbuni 0% 25% 50% 75% 100%
0.0000 7.8958 14.4542 31.2750 512.3292 このyonbuniに基づいてFareをカテゴリー化します。cut()関数を使います。
df_all$Fare <- cut(df_all$Fare,
breaks = yonbuni,
include.lowest = TRUE,
labels = c("economy", "pre_eco", "biz", "first"))
table(df_all$Fare)
economy pre_eco biz first
337 321 328 323 飛行機のように、エコノミー、プレミアムエコノミー、ビジネス、ファーストクラスに
してみました。
Survivedとの関係を見てみます。
barplot(table(df_all$Survived, df_all$Fare),
legend = TRUE)
やっぱり、エコノミークラスの乗客は生存者は少なく、
ファーストクラスの乗客は生存者は多いですね。
Cabin
次はCabinです。NAの有無を確認します。
sum(is.na(df_all$Cabin))[1] 0NAはありませんでした。
Cabinも適当に20個くらいをみてみます。sample()関数で適当に20個表示します。
set.seed(3)
sample(df_all$Cabin, 20, replace = FALSE) [1] "E77" "" "" "" "" "" "" "" "" "B18" "" "D37"
[13] "" "" "" "" "" "" "" "" “” というのが多いですね。
“” が何個あるか数えてみます。
sum(df_all$Cabin == "")[1] 10141014もあります。df_allは全部で1309の観測数ですので、大半が “” です。
Cabinは、“” とそれ以外 のダミー変数にしましょう。
df_all$Cabin <- ifelse(df_all$Cabin == "", 0, 1)
table(df_all$Cabin)
0 1
1014 295 “” は0、その他は1のダミー変数に変換しました。
Survivedとの関係をみてみます。
tab <- table(df_all$Cabin, df_all$Survived)
names(dimnames(tab)) <- c("Cabin", "Survived")
tab Survived
Cabin 0 1
0 481 206
1 68 136Cabinが0の人はSurvivedが1の人が多く、Cabinが1の人はSurvivedが1の人が多いです。
カイ二乗検定で、この関係が偶然でなく、CabinとSurvivedが統計的に有意に
関係しているかを確認します。
chisq.test(tab)
Pearson's Chi-squared test with Yates' continuity correction
data: tab
X-squared = 87.941, df = 1, p-value < 2.2e-16p-value < 2.23-16 ですので、CabinとSurvivedは関連あるのは間違いないようです。
グラフで視覚化もしてみます。
barplot(tab,
main = "Survived by Cabin",
legend = TRUE)
左棒グラフがSurvived = 0です。黒い色、Cabin = 0が
圧倒的に多いことがわかります。
Embarked
いよいよ最後の変数、Embarkedまでたどり着きました。
いつものように、NAの有無を確認しましょう。
NAの有無はEDA(Exploratory Data Analysis)の「いろはの”い”」ですよね。
sum(is.na(df_all$Embarked))[1] 0NAは無いようです。
tabled()関数で度数を確認します。
table(df_all$Embarked)
C Q S
2 270 123 914 あれ、なんか、空白のものが2個ありますね。
これは、一番度数の多いSにしてしまいます。
idx <- df_all$Embarked %in% c("C", "Q", "S") # 空白の位置を特定
df_all$Embarked[!idx] # idxが上手く動いているか確認[1] "" ""df_all$Embarked[!idx] <- "S" # ""をSに置換
table(df_all$Embarked) # 結果確認
C Q S
270 123 916 うまく空白のものをSに置換できました。
EmbarkedはC, Q, Sの3種類なので、ファクター型に変換します。
df_all$Embarked <- as.factor(df_all$Embarked)
summary(df_all$Embarked) C Q S
270 123 916 ファクター型にすると、summery()関数でもtable()関数のように度数がわかります。
Survivedとの関係をみてみます。
tab <- table(df_all$Embarked, df_all$Survived)
names(dimnames(tab)) <- c("Embarked", "Survived")
tab Survived
Embarked 0 1
C 75 93
Q 47 30
S 427 219Cの乗客はSurvivedが1の人が多いことがわかります。
prop.table()関数で比率のテーブルにしたほうがわかりやすいです。
round(prop.table(tab, margin = 1), 2) Survived
Embarked 0 1
C 0.45 0.55
Q 0.61 0.39
S 0.66 0.34margin = 1 で横方向の比率の合計値が1になるようにして、
round()関数で小数点以下2桁まえの表示にしました。
Cの乗客は生存率は55%です。
カイ二乗検定してみます。
chisq.test(tab)
Pearson's Chi-squared test
data: tab
X-squared = 25.964, df = 2, p-value = 2.301e-06p-valueは2.301e-06 ですので、EmbarkedとSurvivedは関連ありそうです。
ちなみに、format()関数で、scientific = FALSE にすると、
2.301e-06をe-06を使わないで表示できます。
format(2.301e-06, scientific = FALSE)[1] "0.000002301"これで、全ての変数のEDAが終了しました。
その他のEDA
いままでで、全ての変数のEDAとデータの修正が終了しました。
summary()関数で全体を把握しておきましょう。
summary(df_all) PassengerId Survived Pclass Sex Age Ticket
Min. : 1 Min. :0.0000 1:323 female:466 Min. : 0.17 bt:723
1st Qu.: 328 1st Qu.:0.0000 2:277 male :843 1st Qu.:21.00 gt:586
Median : 655 Median :0.0000 3:709 Median :26.00
Mean : 655 Mean :0.3838 Mean :29.30
3rd Qu.: 982 3rd Qu.:1.0000 3rd Qu.:36.50
Max. :1309 Max. :1.0000 Max. :80.00
NA's :418
Fare Cabin Embarked title agena
economy:337 Min. :0.0000 C:270 Master.: 59 Min. :0.0000
pre_eco:321 1st Qu.:0.0000 Q:123 Miss. :256 1st Qu.:0.0000
biz :328 Median :0.0000 S:916 Mr. :736 Median :0.0000
first :323 Mean :0.2254 Mrs. :191 Mean :0.2009
3rd Qu.:0.0000 others : 67 3rd Qu.:0.0000
Max. :1.0000 Max. :1.0000
family
alone:790
small:437
large: 82
SexとTicketはファクター型ですが、どちらも水準は2つしかありません。
なので、Sexはmaleなら1、femaleなら0のダミー変数に、
Ticketはgtなら1、btなら0のダミー変数にしましょう。
df_all$male <- ifelse(df_all$Sex == "male", 1, 0)
df_all$gt <- ifelse(df_all$Ticket == "gt", 1, 0)
df_all$Sex <- NULL
df_all$Ticket <- NULL
summary(df_all) PassengerId Survived Pclass Age Fare
Min. : 1 Min. :0.0000 1:323 Min. : 0.17 economy:337
1st Qu.: 328 1st Qu.:0.0000 2:277 1st Qu.:21.00 pre_eco:321
Median : 655 Median :0.0000 3:709 Median :26.00 biz :328
Mean : 655 Mean :0.3838 Mean :29.30 first :323
3rd Qu.: 982 3rd Qu.:1.0000 3rd Qu.:36.50
Max. :1309 Max. :1.0000 Max. :80.00
NA's :418
Cabin Embarked title agena family
Min. :0.0000 C:270 Master.: 59 Min. :0.0000 alone:790
1st Qu.:0.0000 Q:123 Miss. :256 1st Qu.:0.0000 small:437
Median :0.0000 S:916 Mr. :736 Median :0.0000 large: 82
Mean :0.2254 Mrs. :191 Mean :0.2009
3rd Qu.:0.0000 others : 67 3rd Qu.:0.0000
Max. :1.0000 Max. :1.0000
male gt
Min. :0.000 Min. :0.0000
1st Qu.:0.000 1st Qu.:0.0000
Median :1.000 Median :0.0000
Mean :0.644 Mean :0.4477
3rd Qu.:1.000 3rd Qu.:1.0000
Max. :1.000 Max. :1.0000
SexとTicketは必要ないので、削除しています。
Survived, Age, Cabin, agena, male, gtの5つは数値データです。 この数値データ同士の相関係数をみてみます。
ok_cols <- c("Survived", "Age", "Cabin", "agena", "male", "gt")
ok_rows <- !is.na(df_all$Survived)
round(cor(df_all[ok_rows , ok_cols]), 2) Survived Age Cabin agena male gt
Survived 1.00 -0.06 0.32 -0.09 -0.54 0.31
Age -0.06 1.00 0.28 -0.07 0.10 0.22
Cabin 0.32 0.28 1.00 -0.14 -0.14 0.43
agena -0.09 -0.07 -0.14 1.00 0.06 -0.06
male -0.54 0.10 -0.14 0.06 1.00 -0.12
gt 0.31 0.22 0.43 -0.06 -0.12 1.00被説明変数のSurvivedと相関が強いのは、
正の相関ではCabin, gt
負の相関ではmaleです。
Ageだけが連続数値データですが
AgeとCabinの相関は0.28, Ageとgtの相関は0.22, Ageとmaleの相関は0.10と
どれもそんなに相関は強くありません。
ですので、Age x Cabin, Age x gt, Age x maleの交差項も作成しておきます。
df_all$AgeCabin <- df_all$Age * df_all$Cabin
df_all$Agegt <- df_all$Age * df_all$gt
df_all$Agemale <- df_all$Age * df_all$male
summary(df_all) PassengerId Survived Pclass Age Fare
Min. : 1 Min. :0.0000 1:323 Min. : 0.17 economy:337
1st Qu.: 328 1st Qu.:0.0000 2:277 1st Qu.:21.00 pre_eco:321
Median : 655 Median :0.0000 3:709 Median :26.00 biz :328
Mean : 655 Mean :0.3838 Mean :29.30 first :323
3rd Qu.: 982 3rd Qu.:1.0000 3rd Qu.:36.50
Max. :1309 Max. :1.0000 Max. :80.00
NA's :418
Cabin Embarked title agena family
Min. :0.0000 C:270 Master.: 59 Min. :0.0000 alone:790
1st Qu.:0.0000 Q:123 Miss. :256 1st Qu.:0.0000 small:437
Median :0.0000 S:916 Mr. :736 Median :0.0000 large: 82
Mean :0.2254 Mrs. :191 Mean :0.2009
3rd Qu.:0.0000 others : 67 3rd Qu.:0.0000
Max. :1.0000 Max. :1.0000
male gt AgeCabin Agegt
Min. :0.000 Min. :0.0000 Min. : 0.000 Min. : 0.00
1st Qu.:0.000 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.: 0.00
Median :1.000 Median :0.0000 Median : 0.000 Median : 0.00
Mean :0.644 Mean :0.4477 Mean : 8.353 Mean :14.84
3rd Qu.:1.000 3rd Qu.:1.0000 3rd Qu.: 0.000 3rd Qu.:30.00
Max. :1.000 Max. :1.0000 Max. :80.000 Max. :80.00
Agemale
Min. : 0.00
1st Qu.: 0.00
Median :22.00
Mean :19.33
3rd Qu.:30.00
Max. :80.00
これでEDAは終わりです。
次はモデル構築になります。
モデル構築
train_X, train_Y, test_Xの作成
glmnetでモデルを構築する関数は、glmnet()関数、cv.glmnet()関数です。
どちらもデータフレームではなくて、マトリックス型のオブジェクトしか
処理できないので、はじめに、df_allというデータフレームから、
train_X, train_Y, test_Xというマトリックス型のオブジェクト作成する
必要があります。
トレーニング用データとテスト用のデータを分ける行のインデックスは
ok_rowsというインデックスを作成してありますので、
Xの列を指定するベクトルを作成します。
model_cols <- c("Pclass", "Age", "Fare", "Cabin", "Embarked", "title",
"agena", "family", "male", "gt", "AgeCabin", "Agegt",
"Agemale")glmnetパッケージにある、make()関数でデータフレームからマトリックス型を
生成できます。
train_X <- makeX(df_all[ok_rows, model_cols])
test_X <- makeX(df_all[!ok_rows, model_cols])
train_Y <- df_all[ok_rows, "Survived"]head()関数ではじめの数行、数個のデータを表示して、確認します。
head(train_X) Pclass1 Pclass2 Pclass3 Age Fareeconomy Farepre_eco Farebiz Farefirst Cabin
1 0 0 1 22 1 0 0 0 0
2 1 0 0 38 0 0 0 1 1
3 0 0 1 26 0 1 0 0 0
4 1 0 0 35 0 0 0 1 1
5 0 0 1 35 0 1 0 0 0
6 0 0 1 26 0 1 0 0 0
EmbarkedC EmbarkedQ EmbarkedS titleMaster. titleMiss. titleMr. titleMrs.
1 0 0 1 0 0 1 0
2 1 0 0 0 0 0 1
3 0 0 1 0 1 0 0
4 0 0 1 0 0 0 1
5 0 0 1 0 0 1 0
6 0 1 0 0 0 1 0
titleothers agena familyalone familysmall familylarge male gt AgeCabin Agegt
1 0 0 0 1 0 1 0 0 0
2 0 0 0 1 0 0 1 38 38
3 0 0 1 0 0 0 0 0 0
4 0 0 0 1 0 0 1 35 35
5 0 0 1 0 0 1 0 0 0
6 0 1 1 0 0 1 0 0 0
Agemale
1 22
2 0
3 0
4 0
5 35
6 26head(test_X) Pclass1 Pclass2 Pclass3 Age Fareeconomy Farepre_eco Farebiz Farefirst
892 0 0 1 34.5 1 0 0 0
893 0 0 1 47.0 1 0 0 0
894 0 1 0 62.0 0 1 0 0
895 0 0 1 27.0 0 1 0 0
896 0 0 1 22.0 0 1 0 0
897 0 0 1 14.0 0 1 0 0
Cabin EmbarkedC EmbarkedQ EmbarkedS titleMaster. titleMiss. titleMr.
892 0 0 1 0 0 0 1
893 0 0 0 1 0 0 0
894 0 0 1 0 0 0 1
895 0 0 0 1 0 0 1
896 0 0 0 1 0 0 0
897 0 0 0 1 0 0 1
titleMrs. titleothers agena familyalone familysmall familylarge male gt
892 0 0 0 1 0 0 1 0
893 1 0 0 0 1 0 0 0
894 0 0 0 1 0 0 1 1
895 0 0 0 1 0 0 1 0
896 1 0 0 0 1 0 0 0
897 0 0 0 1 0 0 1 0
AgeCabin Agegt Agemale
892 0 0 34.5
893 0 0 0.0
894 0 62 62.0
895 0 0 27.0
896 0 0 0.0
897 0 0 14.0head(train_Y)[1] 0 1 1 1 0 0うまくできましたね。
Ridge, LASSO, Elastic-Net
glmnet()関数、cv.glmnet()関数でモデルの学習をするとき、
alphaというパラメータがあります。
alpha = 0 だとRidge,
alpha = 1 だとLASSO,
alpha が0より大きく1より小さいとElastic-Netになります。
数式を使わずに雰囲気で説明すると、
Ridgeは全ての変数を残し、係数を縮小する。
LASSOは不要な変数の係数を0にする。
Elastic-NetはRidgeとLASSOのいいとこどり。
ということです。
百聞は一見に如かず、といいますから、ちょっとやってみます。
Ridge(alpha = 0)
rid_mod <- glmnet(train_X, train_Y, family = "binomial",
alpha = 0)
plot(rid_mod, xvar = "lambda", label = TRUE)
右から左へと-Log(λ)が小さくなるにつれて、それぞれの変数の係数が0に
近づいているのがわかります。
13がプラスの係数では一番大きい値からスタートして、
21がマイナスの係数では一番(絶対値が)大きい値からスタートしています。
どんな変数か確認します。
colnames(train_X)[c(13, 21)][1] "titleMaster." "familylarge" titleMaster.とfamilylargeです。titleがMaster.の人は生存しやすく、
familyがlargeの人は生存しにくい、ということが読み取れます。
LASSO(alpha = 1)
las_mod <- glmnet(train_X, train_Y, family = "binomial",
alpha = 1)
plot(las_mod, xvar = "lambda", label = TRUE)
-Log(λ)が小さくなるにつれて、変数の係数が0になり、モデルから脱落している
ことがわかります。22が一番最後まで生き残っています。
colnames(train_X)[22][1] "male"22番目の変数はmaleでした。やはり、男性か女性かが一番重要なようです。
Elasitic-Net(alpha = 0.5) とりあえず、alpha = 0.5でやってみます。
ela_mod <- glmnet(train_X, train_Y, family = "binomial",
alpha = 0.5)
plot(ela_mod, xvar = "lambda", label = TRUE)
Elastic-Netでも-Log(λ)が小さくなるにつれて、生き残る係数が減っていくことが
わかります。
予測の実験(lambda.minとlambda.1seの比較)
glmnetパッケージでモデルを学習して、予測をする、という一連の流れを
確認しましょう。
まず、train_X、train_Yをモデル学習用、テスト用に分けます。
set.seed(1000)
train_idx <- sample(1:nrow(train_X), nrow(train_X) * 0.7,
replace = FALSE)これで、train_idxの行の観測値はトレーニング用のデータ、それ以外は
テスト用のデータにわかれます。
そうしたら、cv.glmnet()関数でモデルを学習します。
とりあえず、alpha = 0.5 にします。
set.seed(1)
fit <- cv.glmnet(train_X[train_idx, ], train_Y[train_idx],
alpha = 0.5,
family = "binomial")fitというモデルが学習できました。
plot()関数でfitをグラフにしてみます。
plot(fit)
モデルの係数を確認します。
cbind(
coef(fit, s = "lambda.min"),
coef(fit, s = "lambda.1se")
)27 x 2 sparse Matrix of class "dgCMatrix"
lambda.min lambda.1se
(Intercept) 0.70719960 0.52656761
Pclass1 0.47091782 0.25199065
Pclass2 . .
Pclass3 -0.37590315 -0.28650832
Age . .
Fareeconomy -0.10083510 .
Farepre_eco . .
Farebiz . .
Farefirst 0.06982751 .
Cabin 0.35230895 0.50553269
EmbarkedC . .
EmbarkedQ 0.55073495 .
EmbarkedS . .
titleMaster. 1.82097365 0.81591013
titleMiss. . .
titleMr. -0.76434703 -0.75918553
titleMrs. 0.23646668 0.09520891
titleothers -0.44546642 .
agena -0.12004937 .
familyalone . .
familysmall . 0.04734862
familylarge -1.87519618 -0.94110635
male -1.14992767 -0.89855394
gt 0.57980031 0.37290835
AgeCabin 0.01700343 0.00680552
Agegt . .
Agemale -0.03368960 -0.02389263s = “lambda.min”のほうは、クロスバリデーションで最も良い性能（最小のCV誤差）
を示したモデルに対応する λ です。
s = “lambda.1se”のほうは、lambda.min と比較してCV誤差が
1標準誤差以内に収まり、かつより大きなλ（＝より単純なモデル）
を選択するための値です。
一般に、lambda.1se で得られるモデルは正則化が強く、過学習しにくいため、
未知データに対する予測性能が安定しやすい傾向があります。
このfitを使って、予測します。
# s = lambda.minのモデルの予測
pred_min <- predict(fit, train_X[-train_idx, ], type = "class",
s = "lambda.min")
# s = lambda.1seのモデルの予測
pred_1se <- predict(fit, train_X[-train_idx, ], type = "class",
s = "lambda.1se")table()関数で答え合わせをします。
table(pred_min, train_Y[-train_idx])
pred_min 0 1
0 135 36
1 15 82table(pred_1se, train_Y[-train_idx])
pred_1se 0 1
0 133 36
1 17 82このケースでは、s = lambda.minのモデルのほうが135+82=217個正解で
s = lambda.1seのモデルが133+82=215個正解ですので、lambda.minのほうが
未知のデータに対しても予測結果はよかったですね。
でも、この中のset.seed()関数の中の数字を変えてみると、違う結果になります。
set.seed()関数は乱数のタネを決める関数です。
set.seed(6) # <- ここを6に変更しました。
train_idx <- sample(1:nrow(train_X), nrow(train_X) * 0.7,
replace = FALSE)
set.seed(4) # <- ここを4に変更しました。
fit <- cv.glmnet(train_X[train_idx, ], train_Y[train_idx],
alpha = 0.5,
family = "binomial")
# s = lambda.minのモデルの予測
pred_min <- predict(fit, train_X[-train_idx, ], type = "class",
s = "lambda.min")
# s = lambda.1seのモデルの予測
pred_1se <- predict(fit, train_X[-train_idx, ], type = "class",
s = "lambda.1se")
# s = lambda.minの予測結果
table(pred_min, train_Y[-train_idx])
pred_min 0 1
0 152 29
1 17 70# s = lambda.1seの予測結果
table(pred_1se, train_Y[-train_idx])
pred_1se 0 1
0 154 27
1 15 72train_idxを作るときのset.seed()を6に、
cv.glmnet()を走らせるときのset.seed()を4に変えてみたところ、
s = lambda.minの正解は、152+70=222個。
s = lambda.1seの正解は、154+72=226個となり、lambda.1seのほうが良いです。
モデルは、同じ性能ならば単純であればあるほど良いモデルですので、
s = lambda.1seのほうが未知のデータに対しては、良いモデルであると言えます。
alphaの選択
cv.glmnet()関数は、lambdaの水準はCV(クロスバリデーション)で調べてくれますが、
alphaの水準は、自分で決めないといけません。
alphaの水準を入力 >>> cv.glmnet()関数でモデル作成 >>> 予測結果
こんな感じで関数を作ってみます。
alpha_check <- function(x, seed) { # xはalpha, seedはset.seed
# 学習用のインデックスを作成
set.seed(seed)
train_idx <- sample(1:nrow(train_X), nrow(train_X) * 0.7,
replace = FALSE)
# モデル学習
set.seed(seed)
fit <- cv.glmnet(train_X[train_idx, ], train_Y[train_idx],
alpha = x,
family = "binomial")
# s = lambda.1seのモデルの予測
pred_1se <- predict(fit, train_X[-train_idx, ], type = "class",
s = "lambda.1se")
# s = lambda.minの予測結果(正解率)
mean(pred_1se == train_Y[-train_idx])
}alpha_check()関数というものを作りました。
alphaの水準とset.seed()の値を入力すると、cv.glmnet()関数でモデルを作り
s = “lambda.1se” での予測の正解率を出力する関数です。
エラーなく動くか試してみます。
alpha_check(0.5, 1)[1] 0.8246269alpha_check(0.5, 2)[1] 0.8134328alpha_check(0.5, 3)[1] 0.8656716うまく動いていますね。
次は、alpha = 0.1 で、set.seed()を1~100まで試して、
その平均値を計算してみましょう
results_0.1 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.1[i] <- alpha_check(0.1, i)
}
# 正解率の平均値
mean(results_0.1)[1] 0.83130683.1%がalpha=0.1で100試したときの平均正解率です。
alpha = 0.2 ではどうでしょうか？
results_0.2 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.2[i] <- alpha_check(0.2, i)
}
# 正解率の平均値
mean(results_0.2)[1] 0.829514983.0%がalpha=0.2で100回試したときの平均正解率です。
alpha = 0.3 ではどうでしょうか？
results_0.3 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.3[i] <- alpha_check(0.3, i)
}
# 正解率の平均値
mean(results_0.3)[1] 0.828619482.9%でした。
alpha = 0.4 ではどうでしょうか？
results_0.4 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.4[i] <- alpha_check(0.4, i)
}
# 正解率の平均値
mean(results_0.4)[1] 0.8280597正解率は82.8%でした。
alpha = 0.5 ではどうでしょうか？
results_0.5 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.5[i] <- alpha_check(0.5, i)
}
# 正解率の平均値
mean(results_0.5)[1] 0.828097正解率は82.8%でした。
alpha = 0.6 ではどうでしょうか？
results_0.6 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.6[i] <- alpha_check(0.6, i)
}
# 正解率の平均値
mean(results_0.6)[1] 0.8280224正解率は82.8%でした。
alpha = 0.7ではどうでしょうか？
results_0.7 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.7[i] <- alpha_check(0.7, i)
}
# 正解率の平均値
mean(results_0.7)[1] 0.8272015正解率は82.7%でした。
alpha = 0.8ではどうでしょうか？
results_0.8 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.8[i] <- alpha_check(0.8, i)
}
# 正解率の平均値
mean(results_0.8)[1] 0.8272015正解率は82.7%でした。
alpha = 0.9ではどうでしょうか？
results_0.9 <- numeric(0) # 結果を入れるベクトル
# for loopで100回計算
for (i in 1:100) {
results_0.9[i] <- alpha_check(0.9, i)
}
# 正解率の平均値
mean(results_0.9)[1] 0.8262687正解率は82.6%でした。
alpha = 0.1 ~ 0.9まで試したところ、alpha = 0.1が一番正解率が高かったです。
lambdaの選択
cv.glmnet()関数では、通常、約100個のlambdaを試します。
今回は、300個のlambdaを試すようにしてみます。
# lambdaの300個のシークエンスを作成
lambda_seq <- exp(seq(-4, 4, length.out = 300))
set.seed(1)
fit <- cv.glmnet(train_X, train_Y,
family = "binomial",
alpha = 0.1,
lambda = lambda_seq,
nfolds = 10
)plot()関数でグラフにしてみます。
plot(fit)
-Log(λ)が -1よりも小さいところは、選択された変数が0になっていますね。
少し範囲を変えて、0から5の範囲でもう一度やってみます。
# lambdaの300個のシークエンスを作成
lambda_seq <- exp(seq(-5, 0, length.out = 300))
set.seed(1)
fit <- cv.glmnet(train_X, train_Y,
family = "binomial",
alpha = 0.1,
lambda = lambda_seq,
nfolds = 10
)
plot(fit)
モデルの係数を確認します。
cbind(
coef(fit, s = "lambda.1se"),
coef(fit, s = "lambda.min")
)27 x 2 sparse Matrix of class "dgCMatrix"
lambda.1se lambda.min
(Intercept) 0.526590835 0.918469472
Pclass1 0.342299340 0.482861436
Pclass2 . .
Pclass3 -0.398277481 -0.487163822
Age -0.004183394 -0.011106036
Fareeconomy -0.101161027 -0.255962115
Farepre_eco . -0.124329754
Farebiz . 0.096258985
Farefirst 0.112332513 0.288922389
Cabin 0.423721490 0.552693847
EmbarkedC 0.028737394 .
EmbarkedQ 0.003075255 0.191153051
EmbarkedS -0.148921485 -0.195027781
titleMaster. 0.822192007 1.634072238
titleMiss. 0.378168715 0.339783040
titleMr. -0.672372460 -0.811591741
titleMrs. 0.491161847 0.680165743
titleothers -0.258136326 -0.578104210
agena -0.042845860 -0.224342346
familyalone . 0.293335274
familysmall 0.095551331 .
familylarge -1.217302805 -2.231859167
male -0.706278918 -1.000646939
gt 0.407945591 0.551208462
AgeCabin 0.003963810 0.004090083
Agegt 0.001937733 0.001581456
Agemale -0.019313711 -0.023996678alpha = 0.1 とRidge回帰に近いモデルですので、多くの変数が残っています。
s = “lambda.1se” の係数をグラフにしてみます。
# グラフのデータ
variables <- coef(fit, s = "lambda.1se") # 係数のマトリックス
var_values <- as.vector(variables)[-1] # 係数の数値(Interceptは除外)
var_names <- rownames(variables)[-1] # 係数の名前
top_values <- var_values[abs(var_values) > 0.2] # abs(0.2)より大きい係数のみ
top_names <- var_names[abs(var_values) > 0.2] # abs(0.2)より大きい変数名
# グラフ作成
par(mar = c(4, 12, 4, 2)) # 左側の余白を広げる（横向き用）
barplot(
top_values,
names.arg = top_names,
las = 1, # 横向きラベル
horiz = TRUE, # ★ これで 90 度回転（横棒グラフ）
main = "モデルの係数"
)
par(mar = c(5.1, 4.1, 4.1, 2.1)) # 初期設定に戻すfamily-largeやmale, title-Mr.などは生存にはマイナスの影響、
title-Mrs.やtitle-Miss., title-Master.などは生存にプラスの影響だとわかります。
テストデータで予測
cv.glmnet()関数でモデルができました。
いよいよ、テストデータで予測です。
predict()関数を使います。
pred <- predict(fit, newx = test_X, type = "class",
s = "lambda.1se")
head(pred) lambda.1se
892 "0"
893 "1"
894 "0"
895 "0"
896 "1"
897 "0" これをCSVファイルに出力します。
data.frame()関数でデータフレームを作成して、
write.csv()関数でCSVファイルとして出力します。
results <- data.frame(
PassengerId = rownames(pred),
Survived = pred[, 1]
)
head(results) PassengerId Survived
892 892 0
893 893 1
894 894 0
895 895 0
896 896 1
897 897 0write.csv(results, "submit.csv", row.names = FALSE)Kaggleに提出してみたら、スコアは。。。
base R と glmnetで挑んだスコアとしては、まあまあの結果かなと思います。
これでおしまいです。
読んでいただき、ありがとうございました。
