機械学習の特徴量スケーリングの簡単な観察

この記事は [NTTコムウェア Advent Calendar 2021](https://qiita.com/advent-calendar/2021/nttcomware) 19日目の記事です。

エンジニアリング系の記事が多い中、この記事では機械学習に関する話をしてみようと思います。

機械学習をやってみよう！と学び始めると、ライブラリを使えるようになるまでは比較的簡単なのですが、そこから精度を上げたり正確な（正しい）分析に近づけようとすると、途端に経験知が必要になる部分が増えてくる印象です。
色々ノウハウが書かれている記事もありとても参考になりますが、「百聞は一見に如かず」ということで実際に特徴量スケーリングについて観察してみました。

## この記事に書かれていること
1. 機械学習で扱う特徴量スケーリングのうち、標準化、正規化、対数変換についての簡単な紹介
2.  サンプルデータをもとに、特徴量スケーリングが予測精度に影響するアルゴリズム・影響しないアルゴリズムを観察し、サンプルデータに対して効果的だったスケーリング手法の考察


## 1. 前書き：機械学習で出てくる特徴量のスケーリング
機械学習では、数値変数を扱いやすくするために何らかの変換処理をかけることがあります。この記事ではそんな変換処理のうち、スケーリングと呼ばれるような変換をテーマにします。

細かい説明は省略し、まずはなんとなくのイメージを掴むことを試みます。
例えば、機械学習で家賃を予測することを考えてみます。このとき、入力データ（特徴量、説明変数）として「部屋数」と「駅からの距離 (m)」を使い[^1]、それぞれの値が

* 部屋数：だいたい1部屋 ~ 5部屋くらい
* 駅からの距離：だいたい100m ~ 2000mくらい

というような範囲をとるとします。このようなデータをそのまま線形回帰にかけると、

* 家賃 [万円] = 部屋数 × _1.50_ + 駅からの距離 × (_－0.003_)</u> + 8.00

というような予測式が得られます[^2]。

ここで問題になるのが、斜体にしている係数部分です。それぞれ学習によって推定された数値なのですが、1.50と－0.003では数値の大きさがかなり異なります。このように係数の数値の大きさに差があると、

* 過学習を防ぐために係数の大きさにペナルティを課すとき（正則化）、ペナルティの効きやすさが変わってしまう[^3]
* 学習がうまく進まないケースがある

といった問題が出てくることがあります。（このあたりは書籍『Kaggleで勝つデータ分析の技術』に詳しくまとまっています [[参考文献1]](#%E5%8F%82%E8%80%83%E6%96%87%E7%8C%AE)）

このような係数の大きさの違いは、もともとの特徴量の数値の大きさの違いによるものなので、係数を平等に評価するために

* 特徴量の各数値列を平均0、分散1に変換する（標準化）
* 特徴量の各数値列を最小値0、最大値1に変換する（正規化）

ような変換がよく用いられます[^4]。また、裾が長い分布のデータ[^5]に対しては対数変換も用いられることがあります。

このあたりの数値変換を実際に試してみるのが今回のテーマです。


## 2. 実際にやってみた
タイトルの通りですが、いくつかのアルゴリズムに対して特徴量スケーリングの影響を観察してみます。実はツリー系アルゴリズムは特徴量スケーリングの影響をほぼ受けないことが知られていますが[^6]、それを実際に確認してみるのもサブテーマです。

### 2-1. データ準備
#### 2-1-1. 使用データ
今回はscikit-learnライブラリから入手できる[カリフォルニア住宅価格データセット](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.fetch_california_housing.html) を使用します。
数値変数8個、20,640レコードのデータセットですが、緯度経度は扱いにくいので除外し、残りの6個の特徴量 (MedInc, HouseAge, AveRooms, AveBedms, Population, AveOccup) から住宅価格を予測します。

```python
# データ読み込み
from sklearn.datasets import fetch_california_housing
import pandas as pd

df = pd.DataFrame(housing.data, columns = housing.feature_names)
target = pd.DataFrame(housing.target, columns = ['housing_price'])
```


#### 2-1-2. 特徴量スケーリング方法
* 特徴量の各数値列を平均0、分散1に変換する（標準化）
    * $x' = \displaystyle \frac{x - \mu}{\sigma}$ 　　　　　($\mu$：平均,  $~\sigma$：標準偏差)
* 特徴量の各数値列を最小値0、最大値1に変換する（正規化）
    * $x' = \displaystyle \frac{x - x_{\min}}{x_\max - x_\min}$　　 ($x_\min$：最小値,  $~x_\max$：最大値)
* 特徴量を対数変換してから標準化

の3パターンを試しました。いずれもscikit-learnのライブラリを使うことで簡単に実施できます。

```python
from sklearn import preprocessing
import numpy as np

# 標準化変換
ss = preprocessing.StandardScaler()
df_std = pd.DataFrame(ss.fit_transform(df), 
                      index = df.index, columns = df.columns)

# 正規化
mm = preprocessing.MinMaxScaler()
df_mm = pd.DataFrame(mm.fit_transform(df), 
                     index = df.index, columns = df.columns)

# 対数変換＋標準化
df_log = df.apply(np.log)
ss = preprocessing.StandardScaler()
df_log_std = pd.DataFrame(ss.fit_transform(df_log), 
                          index = df_log.index, columns = df_log.columns)
```

"MedInc" 列を例に変換前後の様子を見てみます。変換前の分布はこのような形をしています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/91b52355-e6ec-3c3b-ab8d-4941f1aec04b.png)
それぞれの変換をすると、分布は以下のように変わります（左から標準化、正規化、対数変換＋正規化）。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/2506527e-a1a2-45ba-4057-6116ccf1ed64.png)

標準化と正規化は線形変換なので分布の形は変わらず、値の範囲のみが変わっていることがわかります。
一方、対数変換は非線形変換なので分布の形も変わり、右裾が長い分布から左右対称に近い形に変わっています。分布に偏りがないほうがよいと言われるので、対数変換をしたケースの効果も確認してみましょう。


### 2-2. アルゴリズムごとの評価をするための準備
#### 2-2-1. 比較アルゴリズム
以下のアルゴリズムを比較します。LightGBM以外はscikit-learnを使用しています。LightGBMはツリー系アルゴリズムなのでスケーリングの影響がないのですが、その確認のために比較対象に入れました。

* 線形回帰（正則化項なし） : `LinearRegression()`
* Elastic Net（正則化項あり線形回帰） : `ElasticNet()`
* ニューラルネットワーク : `MLPRegressor()`
* LightGBM :  `LGBMRegressor()`


#### 2-2-2. 評価方法
それぞれについて20分割クロスバリデーションでRMSEによる精度評価を行い[^7]、得られた精度20個の分布を確認します。まずは「結果に影響するか」のみに着目するため、ハイパーパラメータチューニングは省略します。

```python
from sklearn.model_selection import KFold
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import ElasticNet
from sklearn.neural_network import MLPRegressor
import lightgbm as lgb
from sklearn.metrics import mean_squared_error

### クロスバリデーションで精度検証する関数を定義 ###
def predict_cv(df, target, model):
    kf = KFold(n_splits = 20, shuffle = True, random_state = 42)
    rmse_list = []
    
    for i, (tr_idx, va_idx) in enumerate(kf.split(df)):
        X_train, X_test = df.iloc[tr_idx], df.iloc[va_idx]
        y_train, y_test = target.iloc[tr_idx], target.iloc[va_idx]
        
        model.fit(X_train.values, y_train.values.ravel())
        pred = model.predict(X_test.values)
        
        # 精度(RMSE)
        rmse_list.append(np.sqrt(mean_squared_error(y_test, pred)))
    
    return rmse_list


### 線形回帰 ###
# 無変換
rmse_lr = predict_cv(df, target, LinearRegression())
# 標準化
rmse_lr_std = predict_cv(df_std, target, LinearRegression())
# 正規化
rmse_lr_mm = predict_cv(df_mm, target, LinearRegression())
# 対数変換
rmse_lr_log_std = predict_cv(df_log_std, target, LinearRegression())


### Elastic Net ###
rmse_en = predict_cv(df, target, ElasticNet(random_state = 0))
rmse_en_std = predict_cv(df_std, target, ElasticNet(random_state = 0))
rmse_en_mm = predict_cv(df_mm, target, ElasticNet(random_state = 0))
rmse_en_log_std = predict_cv(df_log_std, target, ElasticNet(random_state = 0))


### ニューラルネットワーク (MLPRegressor) ###
rmse_nn = predict_cv(df, target, 
                     MLPRegressor(hidden_layer_sizes = (30, 30, ), random_state = 1, max_iter = 10000, early_stopping = True))
rmse_nn_std = predict_cv(df_std, target,
                         MLPRegressor(hidden_layer_sizes = (30, 30, ), random_state = 1, max_iter = 10000, early_stopping = True))
rmse_nn_mm = predict_cv(df_mm, target, 
                        MLPRegressor(hidden_layer_sizes = (30, 30, ), random_state = 1, max_iter = 10000, early_stopping = True))
rmse_nn_log_std = predict_cv(df_log_std, target,
                             MLPRegressor(hidden_layer_sizes = (30, 30, ), random_state = 1, max_iter = 10000, early_stopping = True))


### LightGBM ###
rmse_lgb = predict_cv(df, target, lgb.LGBMRegressor(seed = 42))
rmse_lgb_std = predict_cv(df_std, target, lgb.LGBMRegressor(seed = 42))
rmse_lgb_mm = predict_cv(df_mm, target, lgb.LGBMRegressor(seed = 42))
rmse_lgb_log_std = predict_cv(df_log_std, target, lgb.LGBMRegressor(seed = 42))
```

### 2-3. 結果確認
#### 2-3-1. 線形回帰
![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/90bb219e-a4d8-f441-1552-b14319d52a69.png)

* 対数変換をすると予測誤差 (RMSE) の幅が小さくなっていますが、残り3種類の変換は影響がないようです。
* 理論上は、正則化項なしの単純な線形回帰では係数の大きさがどこにも影響しないので、今回の観察でもそれが確認できた形です。

#### 2-3-2. Elastic Net
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/45a1a616-5923-2182-857c-593592f2a169.png)

* それぞれ、変換によって結果が変わっていることがわかります。
* Elastic Netは線形回帰に正則化項（ペナルティ項）がついたものなので、前書きで書いたように係数のスケールがペナルティの効きやすさに影響するためです。
* パラメータチューニングをしていないので、この結果をもって「XXの変換で精度が上がる or 下がる」というような結論はできないことにご注意ください[^8]。「スケーリングが結果に影響するか否か」のみに着目します。

#### 2-3-3. ニューラルネットワーク
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/c32d6d5d-9386-f4dc-c54b-14fa8166124d.png)

* Elastic Net同様に、変換によって結果が変わっていることがわかります。理由もElascit Netと同様です。

#### 2-3-4. LightGBM
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/f5ebd6d9-8d96-c39f-e86d-4f900ff7d892.png)

* 変換有無による影響がほとんどないことがわかります。
* 「ツリー系アルゴリズムは特徴量スケーリングが影響しない」ことが改めて確認できました。


### 2-4. 精度比較
せっかくなので、特徴量スケーリングの影響を受けるElastic Net, ニューラルネットワークについて、ハイパーパラメータチューニングまで含めて比較をしてみました。
まだまだ厳密な比較というには足りない部分もありそうですが[^9]、それぞれの数値スケールに合うように前提条件（ハイパーパラメータ）を揃えるので、多少の参考値にはなるのではないでしょうか。

今回はチューニングにはOptunaを使用し、全データを8:2に分割したうえで試行回数100回でもっとも精度がよかったパラメータを採用しました[^10]。

#### 2-4-1. Elastic Net
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/22fbf29f-e892-e9e2-4171-c3047056ffa4.png)

* 正則化の強さ `alpha` とL1正則化項に対するペナルティの割合 `l1_ratio` をチューニング
* 対数変換をかけると箱ひげ図のひげが短くなっていますが、精度差はほとんどなさそうです。
* 仕組み的にも、Elastic Netでは係数の大きさが影響するのは正則化の強さのみなので、正則化の強さをチューニングすれば標準化、正規化の影響はなくなるのは自然な結果です。


#### 2-4-2. ニューラルネットワーク
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2370717/f6f9d7fe-6add-fdc4-f95e-68dcf07c2fcc.png)

* 隠れ層は2層に固定して2層のユニット数も揃える条件のもとで、ユニット数 (`hidden_layer_sizes`) と正則化の強さ `alpha` をチューニング
* 特徴量スケーリングの効果が表れています。標準化と正規化を見比べると、標準化のほうが少し精度がよい傾向です。対数変換をしたほうが少し箱の位置が低いですが、明確な差があるとまでは言えないですかね。
* アルゴリズムの仕組み的にも、正則化や学習の進みやすさを考えるとスケールを揃えた方がよいので、スケーリングによって精度が良くなるのは自然な結果です。


## 3. まとめ
機械学習における特徴量スケーリングの簡単な紹介と観察をしてみました。

* ツリー系ではなく正則化が効くアルゴリズム（Elascit Netやニューラルネットワーク）では、特徴量スケーリングの有無によって結果が変わる
* 今回のサンプルデータセットでは、ニューラルネットワークにおいて特徴量スケーリングによる精度向上の効果が見られた

あたりが観察で確認できたことです。今後機械学習のアルゴリズムを使うときの参考になれば幸いです。


[^1]: 実際の家賃はもっと複雑な要因によって決まります。良さそうな物件はだいたい家賃も比例していて、うまく価格付けされているものです。
[^2]: 適当に作った数式なので、この予測式はあてにしないでください。
[^3]: 家賃の例で正則化ペナルティを0.5としてみると、部屋数係数に対しては0.75のペナルティが課されるのに対して、距離係数に対しては0.0015のペナルティしか課されなくなります。ペナルティが大きい部屋数係数が必要以上に削減されてしまうことになります。
[^4]: このあたりの用語の使い分けは完全には統一されていないようです。今回は『Kaggleで勝つデータ分析の技術』などで採用されている、比較的多数派と思われる呼び方に合わせました。
[^5]: 日本人の年収の分布のような、ごく一部のデータがかなり極端な値をとるようなデータ
[^6]: 雑な言い方をすると、ツリー系アルゴリズムは「数値がしきい値以上か、しきい値未満か」というIF文の固まりなので、自身以外の特徴量の数値の大きさには影響を受けないです。
[^7]: データをランダムに20グループに分割し、そのうち19グループで学習、残りの1グループを予測して予測精度を確認しています。
[^8]: デフォルトパラメータをそのまま使っているので、変換によってたまたまデフォルトパラメータに適した値（スケール）になっている可能性があります。逆もしかり。
[^9]: そもそも1データセットのみの観察なので、たまたまこのデータセットではこのような結果になった可能性も大いにあります。
[^10]: 厳密にはハイパーパラメータチューニングに使ったデータで精度評価をするのは良くないのですが、簡略化のためこの方法を取りました。


## 参考文献
1. 門脇大輔, 阪田隆司, 保坂桂佑, 平松雄司, 『Kaggleで勝つデータ分析の技術』, 技術評論社, 2019.
2. [Feature Scalingはなぜ必要？ - Qiita](https://qiita.com/ttskng/items/2a33c1ca925e4501e609)
3. [【特徴量スケーリング】いつも紛らわしい「標準化」と「正規化」の違いを理解する - Qiita](https://qiita.com/oki_kosuke/items/02ec7bee92c85cf10ac2)
4. [機械学習でなぜ正規化が必要なのか - Qiita](https://qiita.com/yShig/items/dbeb98598abcc98e1a57)