Gradient Boosting Treeを使ってみる - R{gbm} (Part.2)
前回に続き、いくつかの変数変換を試してみる。
highschoolstudent.hatenablog.com
Partial Dependence Plotで気づいたことを元に、以下の説明変数に対して変換を実行。
"age" - 数値変数だとノイズに敏感に見えるため、カテゴリー化してみる。
"fnlwgt" - 外れ値にキャップをかけてみる。
"capital.gain" - 外れ値にキャップをかけてみる。
"capital.loss" - 元の変数の意味が分からなくて解釈に困るが、カテゴリー化して様子を見てみる。
"hours.per.week" - 外れ値にキャップをかけてみる。
"native.country" - カテゴリカル変数の多水準問題に2つの方法（両方とも数値化する方法）で取り組んでみる。
キャップと言うのは、連続変数に上限下限を設けること。
例えば、"身長"という変数があって180cmで上限のキャップをかけた場合、それ以上高い観測値はすべて180cmに変換する。
Tree系のモデルは説明変数の外れ値にロバストであり、キャップをかけることによってパフォーマンスが大幅に改善したりするとは考えられない気がするが、それでも念のために外れ値の処理はしておいたほうが良いとのこと（同僚より）なので、実行しておく。
カテゴリー化の理由だが、Plotを見た感じで学習データにオーバーフィッティングしている感じがしたので、カテゴリー化することによりスムーズにしようというのが一つで、もう一つが解釈しやすくなる（Plotで見やすくなる）のではないかという期待から。
いずれにせよそこまで科学的な根拠があると言うわけではないが、良くやられているので、とりあえずやってみる、といった要素が強い。
キャップ関数を定義。変数（データフレームの列）を引数とし、ベクトルを返す。
# Capをかける関数、カラムと最小/最大値を指定
fn_cap <- function(num_col, low=min(num_col), high=max(num_col)){
num_col[num_col<=low] <- low
num_col[num_col>=high] <- high
return(num_col)
}
カテゴリー化の関数を定義。変数（データフレームの列）を引数とし、ベクトルを返す。
分け方だが、上記理由を考慮した完全な主観。
# age(numeric)をカテゴリー化する関数
fn_age_category <- function(age_col){
age_col[age_col<20] <- "-19"
age_col[age_col>=20 & age_col<25] <- "20-24"
age_col[age_col>=25 & age_col<30] <- "25-29"
age_col[age_col>=30 & age_col<35] <- "30-34"
age_col[age_col>=35 & age_col<40] <- "35-39"
age_col[age_col>=40 & age_col<45] <- "40-44"
age_col[age_col>=45 & age_col<50] <- "45-49"
age_col[age_col>=50 & age_col<55] <- "50-54"
age_col[age_col>=55 & age_col<60] <- "55-59"
age_col[age_col>=60 & age_col<65] <- "60-65"
age_col[age_col>=65] <- "65-"
return(age_col)
}
# capital.lossをカテゴリー化する関数
fn_capLoss_category <- function(closs_col){
closs_col[closs_col==0] <- "0"
closs_col[closs_col>=0 & closs_col<1000] <- "<1000"
closs_col[closs_col>=1000 & closs_col<1250] <- "<1250"
closs_col[closs_col>=1250 & closs_col<1500] <- "<1500"
closs_col[closs_col>=1500 & closs_col<1750] <- "<1750"
closs_col[closs_col>=1750 & closs_col<2000] <- "<2000"
closs_col[closs_col>=2000] <- ">=2000"
return(closs_col)
}
データ加工は、adult2_2から出発。
後で実行する多水準カテゴリカル変数の変換のときに必要であるからだが、学習データ（学習＋検証）とテストデータを示すフラグ変数を作成しておく。
前回同様、24420までが学習＋検証（0とする）、24421から最後の32561までがテスト（1とする）。
test_col <- c(rep(0,24420),rep(1,(nrow(adult2_2)-24420)))
adult2_3 <- data.frame(adult2_2, Test=test_col)
head(adult2_3)
age workclass fnlwgt education education.num marital.status
1 34 Private 191856 Some-college 10 Married
2 40 Private 205047 HS-grad 9 Married
3 35 Private 230054 7th-8th 4 Married
4 27 Private 279608 5th-6th 3 Married
5 34 Private 226872 Bachelors 13 Divorced
6 61 State-gov 151459 10th 6 Never-married
occupation relationship race sex capital.gain capital.loss
1 Exec-managerial Wife White Female 7298 0
2 Craft-repair Not-in-family White Male 0 0
3 Machine-op-inspct Husband White Male 0 0
4 Farming-fishing Husband White Male 0 0
5 Sales Not-in-family White Female 0 0
6 Other-service Not-in-family Black Female 0 0
hours.per.week native.country income_class income Test
1 40 United-States >50K 1 0
2 40 United-States >50K 1 0
3 40 Mexico <=50K 0 0
4 40 Mexico <=50K 0 0
5 40 United-States <=50K 0 0
6 38 United-States <=50K 0 0
最終列に0/1データの"Test"カラムが作成される。
ageのカテゴリー化。
age2 <- fn_age_category(adult2_3[,"age"])
fnlwgtに500000でキャップ。
fnlwgt2 <- fn_cap(adult2_3[,"fnlwgt"], high=500000)
capital.gainに10000でキャップ。
capital.gain2 <- fn_cap(adult2_3[,"capital.gain"], high=10000)
capital.lossのカテゴリー化。
capital.loss2 <- fn_capLoss_category(adult2_3[,"capital.loss"])
hours.per.weekに80でキャップ。
hours.par.week2 <- fn_cap(adult2_3$hours.per.week, high=80)
データフレームにまとめる。
adult2_4 <- data.frame(adult2_3, age2, fnlwgt2, capital.gain2, capital.loss2, hours.par.week2)
では、カテゴリカル変数の多水準問題（High Cardinality Problem)への対処を考えてみる。
今回は逆に、数値化することによって、オーバーフィッティングへの対処を試みる。
一つ目の方法が、単純にランキング化する方法。Simple Rankingと呼ぶことにする。
学習データにおいて（これがTest列を作成しておいた理由）水準別に目的変数を集計し、値が高いものから順にランキングを付ける。
"native.country"列は42水準。この42水準別に目的変数である"income"が１となる割合を計算する。高いものから順に並べ、1から42のランクを付ける。
関数群を定義しておく。
############################
###### Simple Ranking ######
# 水準別ランキングをdfとして返す
fn_ranking_table <- function(d_f, y_name, x_name, x_new_name){
# 水準別平均値のテーブル
tbl <- tapply(d_f[,y_name], d_f[,x_name], mean)
# ランキング
rnk <- rank(-tbl, na.last="keep", ties.method="average")
# dfとして返す
tbl_df <- data.frame(names(rnk), rnk)
names(tbl_df) <- c(x_name, x_new_name)
return( tbl_df )
}
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "Rankingによる新しい変数名", "テスト変数名")
# テスト変数 - 0:train&validation 1:test
fn_SimpleRanking2 <- function(d_f, y_name, x_name, x_new_name, test_name=NULL){
#テストデータ(水準1)を除く
if(is.null(test_name)){ # テストデータを指定しない場合
d_f_train <- d_f
}
else if(sum(d_f[,test_name]==1)>0) { # テストデータ用の水準(1)を指定した場合
d_f_train <- d_f[d_f[,test_name]!=1,]
}
else { stop("--- Test data set??? ---") } # 不明な変数を指定した場合
#Rankingテーブルの作成
tbl_df <- fn_ranking_table(d_f_train, y_name, x_name, x_new_name)
return( merge(d_f, tbl_df, by=x_name, all.x=TRUE) )
}
############################
############################
fn_SimpleRanking2(データフレーム, "目的変数(numeric)名", "説明変数(カテゴリー)名", "Rankingによる新しい変数名", "テスト変数名(0:学習, 1:テスト)")として実行する。
temp_df <- fn_SimpleRanking2(adult2_4, y_name="income", x_name="native.country", x_new_name="native.country2_SR", test_name="Test")
では、もう一つの方法を実行。
学習データにおいて、説明変数のある水準を除いたときの目的変数の平均値をその水準の値とする。
例えば、"native.country"変数でJapan水準を除いたときの"Income"変数の平均値をJapanの値とする。
さらに、オーバーフィッティングを押さえるため、適当なノイズを加える。
Mean of Response Variable as Explanetory Variable by Leave One Out、略してLOOと呼ぶことにする。
この方法は、同僚に教えてもらった。
###########################
########### LOO ###########
# 説明変数の水準を指定すると、その水準を除いた目的変数の平均値を返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", 説明変数の水準)
fn_mean_oneLevel <- function(d_f, y_name, x_name, x_level){
return( mean(d_f[d_f[,x_name]!=x_level,y_name]) )
}
# fn_mean_oneLevel関数を、指定した説明変数の水準分イテレーション実施して、dfを返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名")
fn_mean_oneLevel_table <- function(d_f, y_name, x_name, x_new_name){
Num <- c()
for(lvl in levels(d_f[,x_name])){
Num <- c(Num, fn_mean_oneLevel(d_f, y_name, x_name, lvl))
}
loo_df <- data.frame(levels(d_f[,x_name]), Num)
names(loo_df) <- c(x_name, x_new_name)
return( loo_df )
}
# 元のdfにLOOを付け、返す関数
# (df, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名", ノイズの大きさ, "テスト変数名")
# テスト変数 - 0:train&validation 1:test
# 学習データのみにある水準が含まれ、テストデータにその水準が含まれない場合の処理はしていない
fn_LOO2 <- function(d_f, y_name, x_name, x_new_name, noise_ratio=0.1, test_name=NULL){
#テストデータ(水準1)を除く
if(is.null(test_name)){ # テストデータを指定しない場合
d_f_train <- d_f
}
else if(sum(d_f[,test_name]==1)>0) { # テストデータ用の水準(1)を指定した場合
d_f_train <- d_f[d_f[,test_name]!=1,]
}
else { stop("--- Test data set??? ---") } # 不明な変数を指定した場合
#LOOテーブルの作成
loo_table <- fn_mean_oneLevel_table(d_f_train, y_name, x_name, x_new_name)
# merge
d_f_merge <- merge(d_f, loo_table, by=x_name, all.x=TRUE, sort=FALSE)
# ノイズのSD
sdiv <- sd(d_f_merge[,x_new_name])
# sdivをベースに正規分布からノイズを作成、noise_ratioで大きさを調整
noise <- rnorm(nrow(d_f), mean=0, sd=sdiv*noise_ratio)
# ノイズを加える
d_f_merge[,x_new_name] <- d_f_merge[,x_new_name] + noise
return( d_f_merge )
}
###########################
###########################
fn_LOO2(データフレーム, "目的変数(numeric)名", "説明変数(カテゴリー)名", "LOOによる新しい変数名", ノイズの大きさ, "テスト変数名(0:学習, 1:テスト)") として実行。
'ノイズの大きさ'とは、0から1の間の値を取り、0だと全くノイズを与えないという意味。
LOOによって各水準の値を計算した後、それらのばらつき（SD）を基準として、そのSDの何割かをノイズとして加えておく。
rnorm(n, mean=0, SD*'ノイズの大きさ')といったノイズを発生させ、LOOで計算した変数に加えておく。
adult2_5 <- fn_LOO2(temp_df, y_name="income", x_name="native.country", x_new_name="native.country2_LOO", noise_ratio=0.1, test_name="Test")
ノイズは0.1にしてみた。
変数が増えてややこしくなってきたので、順序を入れ替えて整理しておく。
adult3 <- adult2_5[,c("Test", "age", "age2", "workclass", "fnlwgt", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain", "capital.gain2", "capital.loss", "capital.loss2", "hours.per.week", "hours.par.week2", "native.country", "native.country2_SR",
"native.country2_LOO", "income_class","income")]
str(adult3)
'data.frame': 32561 obs. of 24 variables: $ Test : int 0 1 0 0 0 1 0 0 0 0 ... $ age : int 26 41 35 20 27 35 38 37 37 36 ... $ age2 : Factor w/ 11 levels "-19","20-24",..: 3 6 5 2 3 5 5 5 5 5 ... $ workclass : Factor w/ 7 levels "?","Federal-gov",..: 4 4 2 1 4 4 4 3 5 4 ... $ fnlwgt : int 56929 207685 35309 200061 269354 187415 220237 244803 188774 31438 ... $ fnlwgt2 : num 56929 207685 35309 200061 269354 ... $ education : Factor w/ 16 levels "10th","11th",..: 10 12 10 16 10 5 9 12 2 12 ... $ education.num : int 13 9 13 10 13 3 11 9 7 9 ... $ marital.status : Factor w/ 5 levels "Divorced","Married",..: 3 2 3 3 2 2 2 2 2 1 ... $ occupation : Factor w/ 13 levels "?","Adm-clerical",..: 4 13 12 1 11 7 3 3 11 13 ... $ relationship : Factor w/ 6 levels "Husband","Not-in-family",..: 2 1 2 4 1 1 1 1 2 5 ... $ race : Factor w/ 4 levels "Asian-Pac-Islander",..: 2 4 1 4 4 1 4 4 4 4 ... $ sex : Factor w/ 2 levels "Female","Male": 2 2 2 1 2 2 2 2 2 2 ... $ capital.gain : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.gain2 : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.loss : int 0 0 0 0 0 0 0 0 0 0 ... $ capital.loss2 : Factor w/ 6 levels "<1000","<1250",..: 1 1 1 1 1 1 1 1 1 1 ... $ hours.per.week : int 50 45 28 40 25 50 40 40 60 43 ... $ hours.par.week2 : int 50 45 28 40 25 50 40 40 60 43 ... $ native.country : Factor w/ 42 levels "?","Cambodia",..: 1 1 1 1 1 1 1 1 1 1 ... $ native.country2_SR : num 19 19 19 19 19 19 19 19 19 19 ... $ native.country2_LOO: num 0.241 0.24 0.241 0.242 0.241 ... $ income_class : Factor w/ 2 levels "<=50K",">50K": 1 1 1 1 1 1 2 1 1 1 ... $ income : int 0 0 0 0 0 0 1 0 0 0 ...
では、再びモデルを作成してみる。
学習(+検証)、テストデータの分割。
train3 <- adult3[adult3[,"Test"]==0,]
test3 <- adult3[adult3[,"Test"]==1,]
"native.country2_SR"と"native.country2_LOO"で、別々にモデルを作ってみる。
（元々の変数重要度が、最初に作ったモデルでも4.893%しかないので、大きく違ってくるとは考えにくい）
まず、Simple Ranking,"native.country2_SR"。
キャップをかけたりカテゴリー化した変数ももちろん用いる。
x_col_names <- c("age2", "workclass", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain2", "capital.loss2", "hours.par.week2", "native.country2_SR")
x_cols <- train3[, x_col_names]
y_col <- train3[,"income"]
train3データセットのうち、最初の16264行を学習、残りを検証。これは前回と同じ。
for_train <- 16264
前回のPatial Dependence Plotの結果を見て、単調増加(monotone increasing)に見えた変数には、その制約を課してみる。
"education.num"と"capital.gain"をmonotone increasingとする。ベクトルに1を入れる。
単調減少(monotone decreasing)の場合は-1を入れる。
"native.country2_SR"は、学習(+検証)データにおいて目的変数の値が高いものから順に並べれれているので、monotone decreasingの制約をかける。
単調増加/減少にすることにより、オーバーフィッティングの回避が行われるのだと思われる。
monotone <- c(0,0,0,0,1,0,0,0,0,0,1,0,0,-1)
0はこの制約をかけないという意味。
引数var.monotoneに、このベクトルを渡す。
mod2_1 <- gbm.fit(
x = x_cols,
y = y_col,
distribution = "bernoulli",
var.monotone = monotone,
bag.fraction = 0.5,
n.trees = 5000,
interaction.depth = 7,
n.minobsinnode = 100,
shrinkage = 0.005,
nTrain = for_train
)
変数重要度。
summary(mod2_1)
var rel.inf relationship relationship 30.1044827 capital.gain2 capital.gain2 15.2426028 occupation occupation 14.1707533 education education 11.5529362 age2 age2 8.6704758 education.num education.num 4.4322248 capital.loss2 capital.loss2 3.9274653 hours.par.week2 hours.par.week2 3.6641404 fnlwgt2 fnlwgt2 3.1119212 workclass workclass 2.3148916 marital.status marital.status 1.4242388 native.country2_SR native.country2_SR 0.8392975 sex sex 0.3368926 race race 0.2076771
native.countryをランク化したことによって、重要度が落ちた。他はあまり変化がない。
最適なTreeの数。
gbm.perf(mod2_1, plot.it=TRUE)
3132
では、予測を行う。
学習データ。
x_cols_train <- x_cols[1:for_train,]
y_col_train <- y_col[1:for_train]
検証データ。
x_cols_varidation <- x_cols[(for_train+1):nrow(x_cols),]
y_col_varidation <- y_col[(for_train+1):nrow(x_cols)]
学習データに対する予測値。
pred_mod2_1_train <- predict(mod2_1, newdata=x_cols_train, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
検証データに対する予測値。
pred_mod2_1_varidation <- predict(mod2_1, newdata=x_cols_varidation, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
テストデータに対する予測値。
test_x <- test3[,x_col_names]
test_y <- test3[,"income"]
pred_mod2_1_test <- predict(mod2_1, newdata=test_x, n.trees=gbm.perf(mod2_1, plot.it=FALSE), type="response")
学習データでのAUC。
fn_AUC(pred=pred_mod2_1_train, obs=y_col_train)
0.9332073
検証データでのAUC。
fn_AUC(pred=pred_mod2_1_varidation, obs=y_col_varidation)
0.9198591
テストデータでのAUC。
fn_AUC(pred=pred_mod2_1_test, obs=test_y)
0.9192887
fn_AUCに関しては、前回参照。
個人的な経験だが、今回示したようなオーバーフィッティングの回避を試みると、学習とテストデータのパフォーマンス(AUCなど)の幅が狭くなり、検証とテストデータのパフォーマンスもかなり近くなるといった傾向がある。
前回の結果と見比べてみてほしい。
一応、LOOも試しておく。
x_col_names <- c("age2", "workclass", "fnlwgt2", "education", "education.num", "marital.status", "occupation", "relationship", "race", "sex", "capital.gain2", "capital.loss2", "hours.par.week2", "native.country2_LOO")
x_cols <- train3[, x_col_names]
y_col <- train3[,"income"]
for_train <- 16264
monotone <- c(0,0,0,0,1,0,0,0,0,0,1,0,0,1)
"native.country2_LOO"はmonotone increasingとしておく。
同様にgbmを実行し、結果をmod2_2に格納。
変数重要度。
summary(mod2_2)
var rel.inf relationship relationship 30.0720371 capital.gain2 capital.gain2 15.2900749 occupation occupation 14.2201974 education education 11.1203147 age2 age2 8.7793286 education.num education.num 5.1393013 capital.loss2 capital.loss2 3.9263051 hours.par.week2 hours.par.week2 3.6922352 fnlwgt2 fnlwgt2 3.0737037 workclass workclass 2.3722413 marital.status marital.status 1.3842674 sex sex 0.3466526 native.country2_LOO native.country2_LOO 0.3382367 race race 0.2451039
最適なTreeの数。
gbm.perf(mod2_2, plot.it=TRUE)
3328
予測を行い、AUCを計算する。
学習データでのAUC。
0.9331306
検証データでのAUC。
0.9196043
テストデータでのAUC。
0.9192822
以上、今回はSimple RankingやLOO関数のメモを個人的理由とする。
