読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

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となる割合を計算する。高いものから順に並べ、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関数のメモを個人的理由とする。