XGBoostとデータ前処理(1) 備忘録

# 概要
Jason先生がきれいに、かつ実用的にまとめられていたので要所を備忘録として残す。本当はこれを紹介する程度にするつもりだったけど、データ前処理にハマったのでタイトルに変更した。
後日、グリッドサーチやBayesian Optimizationを投稿する。

https://machinelearningmastery.com/xgboost-python-mini-course/

- 実施期間: 2021年9月
- 環境：Colab
- パケージ：XGBoost

## Gradient Boosting Decision Tree (GBDT)
Trainingにおいて、勾配降下法(GD)で目的関数の値(Loss)が小さくなるよう、数値最適化で弱いdecision treesを直列につないでいくアルゴリズム。DNNにおいてEpochが進むように、弱いdecision treeが追加されていく。どれくらい追加するかなどはハイパパラメータとなる。
要素としては下記となる。

- 最適化する目的関数としての損失関数：Classificationではcross entropy、RegressionではMSEなどが使用される
- 予測する弱いdecision tree：直列に複数連なったdecision treeで予測する
- 加算モデル：Lossが小さくなるように逐次decision treeが追加されていく。

以前のdecision treesの残留Lossが小さくなるように次のdecision treeが追加される。Lossを小さくするために以前のdecision treesが変更されることはなく、積み重なっていくだけ。

## XGBoost (eXtreme Gradient Boosting)
GBDTの実装の一種類で、特に高速化のため下記の機能が実装されている。

- Trainingにおいて、全CPUコアを使用する並列処理
- クラスタ化されたPCを使用した分散処理
- datasetが大きくメモリに乗らない場合はOut-of-Core処理（外部記憶から適宜読み込み）
- データ構造のキャッシュ最適化

ただ、XGBoostよりもLoghtgbmはさらに高速化がなされている。

XGBoostはColabは最初からインストール済みだが、Colabでなければ下記でインストールできる。

```terminal
# 標準Pythonなら
pip install xgboost

# Anacondaなら
conda install -c conda-forge xgboost
```
## Dataset
今回使用するデータは米国のピマ・インディアンの血を引く21歳以上の女性について、糖尿病の陽性陰性を調査したも
下記から入手可能

https://www.kaggle.com/uciml/pima-indians-diabetes-database?select=diabetes.csv

#### 各列の説明
   1. Pregnancies: 妊娠の回数
   2. Glucose: 経口ブドウ糖負荷試験における2時間後の血漿グルコース濃度
   3. BloodPressure: 拡張期血圧(mmHg)
   4. SkinThickness: 上腕三頭筋皮弁の厚さ(mm)
   5. Insulin: 2時間値の血清インスリン（μU/ml
   6. BMI: 肥満度指数（体重（kg）／（身長（m））^2
   7. DiabetesPedigreeFunction: 糖尿病の血統機能
   8. Age: 年齢(歳)
   9. Outcome: クラス変数(0:陰性(500個)、1:陽性(268個))

合計768個のデータが十分かわからないが、positive/negativeがimbalancedではないので良かった。

### Datasetのロード
```python
import pandas as pd
import numpy as np
from numpy import loadtxt
from xgboost import XGBClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

df = pd.read_csv('diabetes.csv')
df.info()
```

![Screenshot from 2021-09-20 14-51-35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/0777eb08-e152-3e2b-36f8-576c8ed46f7c.png)

## Baseline
Nullはないようなので、このままモデルに入れるとどうなるか見てみよう。
XGBoostモデル用に成形する。

```python
# split data into X and y
X = np.array(df.iloc[:,0:8])
Y = np.array(df.iloc[:,8])

# split data into train and test sets
seed = 7
test_size = 0.33
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=test_size, random_state=seed)
```

### Training
整形したtrainデータでまずaccを求め、無調整のXGBoostのbaselineを確認する。

```python
# fit model on training data
model = XGBClassifier()
model.fit(X_train, y_train)

# make predictions for test data
y_pred = model.predict(X_test)
predictions = [round(value) for value in y_pred]

# evaluate predictions
accuracy = accuracy_score(y_test, predictions)
print("Accuracy: %.2f%%" % (accuracy * 100.0))
```

**Accuracy: 75.20%**
良くはない。
どのように良くないかわからないので、test datasetを使ってTraining中のacc変動を確認する。metricは**error**とする。


```python
eval_set = [(X_test, y_test)]
model.fit(X_train, y_train, eval_metric="error", eval_set=eval_set, verbose=True)
```

[0]	validation_0-error:0.259843
[1]	validation_0-error:0.26378
[2]	validation_0-error:0.26378
：
[83]	validation_0-error:0.212598
[84]	validation_0-error:0.208661
[85]	validation_0-error:0.204724
[86]	validation_0-error:0.212598
：
decision treeが追加されるたびにtest datasetを使ったerrorが低下していっているが、上昇に転じる箇所が確認できる。本番のfit()ではearly stoppingでtreeの追加を止めることとする。metricは**logloss**とする。loglossが底を打ってから様子を見るため10回までfitは継続させる。

```python
# fit model on training data
model = XGBClassifier()
eval_set = [(X_test, y_test)]
model.fit(X_train, y_train, early_stopping_rounds=10, eval_metric="logloss", eval_set=eval_set, verbose=True)

# make predictions for test data
y_pred = model.predict(X_test)
predictions = [round(value) for value in y_pred]

# evaluate predictions
accuracy = accuracy_score(y_test, predictions)
print("Accuracy: %.2f%%" % (accuracy * 100.0))
```

[0]	validation_0-logloss:0.660186
Will train until validation_0-logloss hasn't improved in 10 rounds.
[1]	validation_0-logloss:0.634854
[2]	validation_0-logloss:0.61224
：
[54]	validation_0-logloss:0.463379
[55]	validation_0-logloss:0.462555
[56]	validation_0-logloss:0.463003
[57]	validation_0-logloss:0.462065
[58]	validation_0-logloss:0.463195
[59]	validation_0-logloss:0.464264
[60]	validation_0-logloss:0.464393
[61]	validation_0-logloss:0.464292
[62]	validation_0-logloss:0.4644
[63]	validation_0-logloss:0.464663
[64]	validation_0-logloss:0.464624
[65]	validation_0-logloss:0.464717
[66]	validation_0-logloss:0.465197
[67]	validation_0-logloss:0.465055
Stopping. Best iteration:
[57]	validation_0-logloss:0.462065
**Accuracy: 77.17%**

Early stoppingでは改善しなかった。pre-processしなかったので当然データに問題がある。

### Feature importance
XGBoostにはどの説明変数を参照していたか確認する機能があるので、現時点のそれを確認する。

```python
from xgboost import plot_importance
from matplotlib import pyplot

_, ax = plt.subplots(figsize=(12, 4))
plot_importance(model, ax=ax,
                importance_type='weight',
                show_values=True)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/0cda7d12-a373-75ca-14fa-1a7a1efe9f9e.png)

説明変数(f0 ~ f7)が見づらいが下記の順である。
['Pregnancies', 'Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI', 'DiabetesPedigreeFunction', 'Age', 'Outcome']

簡単のため、early stoppingでtrainしてfeature importanceを描画するまでを関数にして、以降これを使う。

```python
def xgb_train(my_df):
    num = len(my_df.columns)-1
    X = np.array(my_df.iloc[:,0:num])
    Y = np.array(my_df.iloc[:,num])

    # split data into train and test sets
    seed = 27
    test_size = 0.33
    X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=test_size, random_state=seed)

    eval_set = [(X_test, y_test)]

    # fit model on training data
    model = XGBClassifier()
    eval_set = [(X_test, y_test)]
    model.fit(X_train, y_train, early_stopping_rounds=10, eval_metric="logloss", eval_set=eval_set, verbose=True)

    # make predictions for test data
    y_pred = model.predict(X_test)
    predictions = [round(value) for value in y_pred]

    # evaluate predictions
    accuracy = accuracy_score(y_test, predictions)
    print("Accuracy: %.2f%%" % (accuracy * 100.0))

    _, ax = plt.subplots(figsize=(12, 4))
    plot_importance(model, ax=ax,
                    importance_type='weight',
                    show_values=True)
    plt.show()
```


## EDA
初めて使うデータなので説明変数の分布を確認する。

```python
import matplotlib.pyplot as plt
import seaborn as sns
sns.pairplot(df.iloc[:,0:8])
plt.show()

df.iloc[:,0:8].describe()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/b544af9f-4cf6-85b8-689d-b2939d11f237.png)

![Screenshot from 2021-09-26 09-49-46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/60ec181d-ca8d-d658-d7d5-3e94e97eed0a.png)
ありえない説明変数にゼロが散見されるので調べてみる。


```python
clmn = list(df.columns)
for i in clmn:
    print(f'{i}: {sum(df[i] == 0)}')
```

![Screenshot from 2021-09-26 08-40-16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/8fcdc498-e997-e1cf-1498-e599de33704b.png)
'Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI'の５つの説明変数はゼロになりえないのでデータの欠損とみなされ、これらのせいで精度が悪かったのは間違いない。
スポットのヒストグラムで確認する。

```python
df_hist = df.loc[:,['Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI']]
df_hist.hist(figsize = (8,8), bins=20)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/d26c2f8d-c0b1-f35a-8dc0-54da8b676975.png)
SkinThicknessとInsulinのゼロが特に目立つ。Feature importanceで両方の出現回数が低かったのはこのゼロのためか。

作戦としては次の２つが思いつく。
**作戦1. 分布を参考にmeanやaverageで全ゼロを置き換える**
**作戦2. ゼロを含むデータを容赦なく削除する**

## 作戦１
### ゼロの補完
ヒストグラムより分布が左右対象なものは現分布の平均値で置き換え、裾野が片方に伸びたものは現分布の中央値で置き換えることとする。
つまり、'Glucose', 'Insulin'は中央値、'BloodPressure', 'SkinThickness', 'BMI'は平均値で置き換える。

```python
Glucose_median = int(np.median(df[df['Glucose'] != 0]['Glucose']))
Insulin_median = int(np.median(df[df['Insulin'] != 0]['Insulin']))
BloodPressure_ave = int(np.average(df[df['BloodPressure'] != 0]['BloodPressure']))
SkinThickness_ave = int(np.average(df[df['SkinThickness'] != 0]['SkinThickness']))
BMI_ave = np.average(df[df['BMI'] != 0]['BMI'])
print(f'{Glucose_median}, {Insulin_median}, {BloodPressure_ave}, {SkinThickness_ave}, {BMI_ave}')
```

117, 125, 72, 29, 32.457463672391015

```python
df['Glucose'] = df['Glucose'].replace(0, Glucose_median)
df['Insulin'] = df['Insulin'].replace(0, Insulin_median)
df['BloodPressure'] = df['BloodPressure'].replace(0, BloodPressure_ave)
df['SkinThickness'] = df['SkinThickness'].replace(0, SkinThickness_ave)
df['BMI'] = df['BMI'].replace(0.0, BMI_ave)

clmn = list(df.columns)
for i in clmn:
    print(f'{i}: {sum(df[i] == 0)}')

df_hist = df.loc[:,['Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI']]
img = df_hist.hist(figsize = (8,8), bins=20)
```

Pregnancies: 111
Glucose: 0
BloodPressure: 0
SkinThickness: 0
Insulin: 0
BMI: 0
DiabetesPedigreeFunction: 0
Age: 0
Outcome: 500
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/968160cb-4f13-c889-6a4f-8268c5776ae6.png)
ゼロは取れたが、SkinThicknessとInsulinの分布がいびつで気になる。
Feature importanceがどれだけ改善されるか、これでモデルに入れて精度を確認する。

### Training

```python
xgb_train(df)
```

:
[56]	validation_0-logloss:0.44372
[57]	validation_0-logloss:0.443915
[58]	validation_0-logloss:0.44359
[59]	validation_0-logloss:0.445002
[60]	validation_0-logloss:0.445071
[61]	validation_0-logloss:0.445515
[62]	validation_0-logloss:0.446372
[63]	validation_0-logloss:0.446125
[64]	validation_0-logloss:0.446569
[65]	validation_0-logloss:0.445882
[66]	validation_0-logloss:0.445944
[67]	validation_0-logloss:0.445446
[68]	validation_0-logloss:0.444948
Stopping. Best iteration:
[58]	validation_0-logloss:0.44359

Accuracy: **79.53%**

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/20697da2-3713-2b61-2442-72e081039d29.png)
改善効果：　77.17% -> 79.53%
また、Feature importanceは上図となった。説明変数(f0 ~ f7)が見づらいが下記の順である。
['Pregnancies', 'Glucose', 'BloodPressure', 'SkinThickness', 'Insulin', 'BMI', 'DiabetesPedigreeFunction', 'Age', 'Outcome']


SkinThickness(f3)、Insulin(f4)のtree出現回数（参照回数）はbaselineより多くなり改善と言える。
SkinThickness(f3)が22回と低いが、そもそも出現回数の多いBMI(f5)とSkinThicknessに相関があるような気がするので確認する。


```python
print(np.corrcoef(np.array(df['SkinThickness']), np.array(df['BMI'])))
print(np.corrcoef(np.array(df['Insulin']), np.array(df['BMI'])))
```

[[1.         0.54316229]
 [0.54316229 1.        ]]
[[1.         0.18016977]
 [0.18016977 1.        ]]

BMIとSkinThicknessはやや相関はあると言える。そのためBMIに食われてSkinThicknessの出現回数が減ったと思われる。
(Insulin-BMIの相関は参考値)
説明変数からSkinThicknessを除いてもよいのかもしれない。

### 冗長な説明変数(SkinThickness)の削除

```python
df_copy = df.copy(deep = True)
df_copy = df_copy.drop('SkinThickness', axis=1)
xgb_train(df_copy)
```

[48]	validation_0-logloss:0.444248
[49]	validation_0-logloss:0.443776
[50]	validation_0-logloss:0.442832
[51]	validation_0-logloss:0.443451
[52]	validation_0-logloss:0.443629
[53]	validation_0-logloss:0.44426
[54]	validation_0-logloss:0.44394
[55]	validation_0-logloss:0.444451
[56]	validation_0-logloss:0.44435
[57]	validation_0-logloss:0.446117
[58]	validation_0-logloss:0.444625
[59]	validation_0-logloss:0.444344
[60]	validation_0-logloss:0.443666
Stopping. Best iteration:
[50]	validation_0-logloss:0.442832

Accuracy: **80.71%**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/84d63a00-c06c-443c-096a-0683a871f5d2.png)
説明変数(f0 ~ f6)が見づらいが下記の順である。
['Pregnancies', 'Glucose', 'BloodPressure', 'Insulin', 'BMI', 'DiabetesPedigreeFunction', 'Age', 'Outcome']

さらに良くなり下記である。やはりSkinThicknessは不要だった。
改善効果：　**77.17% -> 79.53% -> 80.71%**
なんとなく、'Insulin'のゼロの置き換え方がイケてないような気がしてならない。今日はココまでとして他に手法がないか後日調べてみよう。

## 作戦２
### ゼロの削除
まず、SkinThicknessが冗長であることがわかったので、これは使わない。
従い、先の４つの説明変数についてゼロが含まれるデータを削除する。


```python
def multiply_list(lst):
    result = 1
    for i in lst:
        result *= i
    return result

df = pd.read_csv('diabetes.csv')
df = df.drop('SkinThickness', axis=1)
arr = np.array(df.loc[:, ['Glucose', 'BloodPressure', 'Insulin', 'BMI']])
df_del = pd.DataFrame(columns=df.columns)

# ゼロでなければ何らかの数値が入っている。
for i in range(len(arr)):
    if (multiply_list(arr[i]) != 0):
        print(list(df.iloc[i]))
        df_del = df_del.append(df.iloc[i])

print(len(df_del))
df_del.describe()

sns.pairplot(df_del, hue = 'Outcome', diag_kind='hist')
plt.show()
```

![Screenshot from 2021-09-26 16-28-22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/487f6d8b-db41-6e75-aaf4-006973a662ca.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/eb39c596-f51e-25ee-25bd-3f481c6e8888.png)

データ数は392個に減ってしまったが、分布はとても自然

```python
xgb_train(df_del)
```

:
[26]	validation_0-logloss:0.512833
[27]	validation_0-logloss:0.511745
[28]	validation_0-logloss:0.510679
[29]	validation_0-logloss:0.513842
[30]	validation_0-logloss:0.514283
[31]	validation_0-logloss:0.516222
[32]	validation_0-logloss:0.518303
[33]	validation_0-logloss:0.51677
[34]	validation_0-logloss:0.514266
[35]	validation_0-logloss:0.516112
[36]	validation_0-logloss:0.515658
[37]	validation_0-logloss:0.518273
[38]	validation_0-logloss:0.52124
Stopping. Best iteration:
[28]	validation_0-logloss:0.510679

Accuracy: **76.15%**

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/e1e447dd-f76c-fb31-44c5-54a0a3c7a205.png)

説明変数(f0 ~ f6)が見づらいが下記の順である。
['Pregnancies', 'Glucose', 'BloodPressure', 'Insulin', 'BMI', 'DiabetesPedigreeFunction', 'Age', 'Outcome']

baselineよりも悪くなってしまった。データ数が少なくなりすぎたせいだろう。

以上