特徴量選択のまとめ

Kaggle Advent Calendar その2の23日目の記事です。

私はkaggleを始めたばかりでテーブルデータのコンペはTitanicしかやったことがないため、特徴量をどのように選べばいいのかよくわからなかったのでまとめます。

# 特徴量選択手法のまとめ

特徴量選択とは、機械学習のモデルを使用する際に有効な特徴量の組み合わせを探索するプロセスのことを表しています。

特徴量選択を行うことによりいくつかのメリットが得られます。

- 変数を少なくすることで解釈性を上げる
- 計算コストを下げて、学習時間を短縮する
- 過適合を避けて汎用性を向上させる
- 高次元データによって、パフォーマンスが下がることを防ぐ。

## 特徴量選択の種類

特徴量選択の手法は大別して3つ存在します。

- Filter Method
- Wrapper Method
- Emedded Method

### Filter Method

**Filter Method**とは、機械学習モデルを使用せずにデータセットのみで完結する手法であり、データの性質に依存します。そのためたいていはどの機械学習モデルに対しても有効であり、なおかつ処理が高速です。

まずは変数の1つ1つを独立に処理していく手法を紹介します。手順としては以下になります。

1. 評価指標に則って、特徴量を1つ1つランク付けする
2. 上位のランクの特徴量を選択して使用する

あくまで変数を独立に計算していくので、変数間の関係などは考えません。

評価指標として、変数ごとに統計テストを行っていきます。

- カイ2乗スコア
- フィッシャースコア
- anova
- 変数の分散

また、変数間同士の関係も計算する手法も存在し、例として以下のようなものになります。

- 無駄な(説明変数に影響を与えない)特徴量を削除する
- 同じ傾向を示す($y1 = x1 + a$のような)特徴量を削除する
- 相関が非常に高い特徴量を削除する

この**Filter Method**は、特徴量を削減するときに真っ先に行う処理です。

### Wrapper Method

Filter Methodでは、特徴量の性質を1つ1つ見ていき、それが予測にどの程度寄与しているのか等は考慮していませんでした。しかし実際には、最適な特徴量の組み合わせは使用する機械学習モデルに依存します。<br>

**Wrapper Method**とは、機械学習モデルを使用して特徴量の組み合わせを評価することです。こうすることでFilter Methodではわからなかった変数間の関係を探し出し、それぞれのモデルに最適な特徴量の組み合わせを探し出すことができます。

Wrapper Methodを行うには以下を繰り返して実行していきます。

1. 特徴量の組み合わせを選択する
2. [1]で選択した特徴量を使用してモデルを学習させる
3. 性能を評価する

しかしどのように特徴量の組み合わせを探せばいいのか、またいづ探索を切り上げればいいのでしょうか。

特徴量の組み合わせを見つけ出す方法は大きく3つ存在しています。

1. Forward Feature Selection
    - イテレーションごとに特徴量を1つずつ**追加**していく手法
2. Backward feature Elimination
    - イテレーションごとに特徴量を1つずつ**削除**していく方法
3. Exhaustive Feature Search : 
    - すべての組み合わせを試す

この組み合わせの探索方法からわかるように、Wrapper Methodは、Filter Methodと比較して、計算コストが非常に高いです。しかしFilter Methodよりも精度の高い特徴量選択を行うことができます。

また組み合わせの探索を切り上げるための評価指標は、タスクに合わせて選択する必要があります。

### Embedded Method

**Embedded Method**とは、特徴量選択をモデルの学習時に行う手法です。Filter Methodでは計算することのできなかった変換の関係も、この手法で計算することができます。また学習時に探索するため、Wrapper Methodよりも計算コストはかなり低いです。

手順は以下のようになります。

1. モデルを学習させる
2. 特徴量の重要度を算出する
3. 重要でない特徴量を削除する

よく使用されるモデルとして、**Lasso**モデルや、**木モデル**での`feature importance`、**回帰係数**を含む線形モデルなどが存在します。

## 実際に解析で使用してみる

使用するデータは2018年に開催された[Santandar](https://www.kaggle.com/c/santander-value-prediction-challenge)のデータを使用します。

```python
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.feature_selection import VarianceThreshold

data = pd.read_csv('train.csv')
print(data.shape)
-> (4459, 4993)
```

データの訓練数と特徴の数は確認できたので、今度はNULL値が含まれているかどうかを確認します。

```python
null_cols = [col for col in data.columns if data[col].isnull().sum() > 0]
print(null_cols)
-> []
```

NULL値が含まれていないことが確認できました。次は解析に使用するために、データを学習データと検証データに分けていきます。

```python
X_train, X_test, y_train, y_test = train_test_split(
    data.drop(['ID', 'target'], axis=1),
    data['target'],
    test_size=0.2,
    random_state=2018
)
```

## Filter Methodでの解析例

### 特徴量の値のみ選択する

ここでは特徴量の分散や同じ特徴量が存在するのかを確認していきます。

なお処理は**個別に行っているわけではなく、上から順番に**処理を行っています。

#### 分散が0の場合

```python
# 分散が0（すべて同じ値）のデータは削除します
sel = VarianceThreshold(threshold=0)
sel.fit(X_train)

# get_supportで保持するデータのみをTrue値、そうでないものはFalse値を返します
print(sum(sel.get_support()))
# -> 4692
```

これで元の4993個存在していた特徴量が、4692個とおよそ300個ほど特徴量を削減することができました。<br>
あとはこの操作を訓練データとテストデータの両方に適用させます。

```python
# numpy に変換したい場合はこちら
X_train = sel.transform(X_train)
X_test = sel.transform(X_test)

# pandasのまま保持して次の処理を行い場合はこちら
X_train = X_train.loc[:, sel.get_support()]
X_test = X_test.loc[:, sel.get_support()]
```

#### 分散がほぼ0の場合

ここでは特徴量の値の大部分が同じデータのみの場合に削除するようにしていきます。当然このやり方は必ず正しいというわけではなく、元のデータの特徴を観察して削除するかどうか決めてください。

```python
sel = VarianceThreshold(threshold=0.1)
sel.fit(X_train)

print(sum(sel.get_support()))
# -> 4692
```

```python
# numpy に変換したい場合はこちら
X_train = sel.transform(X_train)
X_test = sel.transform(X_test)

# pandasのまま保持して次の処理を行い場合はこちら
X_train = X_train.loc[:, sel.get_support()]
X_test = X_test.loc[:, sel.get_support()]
```

この処理では特徴量の数は削減されなかったので、分散が0.1以下の特徴量は存在していないことがわかります。

#### 特徴量がほかの特徴量と完全に一致している場合

ここでは`Pandas`の`duplicated`メソッドを使用します。

```python
# indexとcolumnsを入れ替える
X_train_T = X_train.T

print(X_train_T.duplicated().sum())
# -> 9

# 同じ特徴量の名前を取得したい場合はこちら
duplicated_features = X_train_T[X_train_T.duplicated()].index.values

print(len(X_train.columns))
# -> 4683
```

この手法で9つの特徴量が削減できました。

### 特徴量の相関係数を使用する

よい特徴量ならば、ターゲットの値と高い相関を持っているという仮定で特徴量を削減していきます。<br>
注意点としては、特徴量にはターゲットといい相関を持っていてほしいですが、同時にほかの特徴量同士はあまり相関を持たないようにすることです。こうすることで冗長な特徴量を削減することができます。

参考：[Correlation-based Feature Selection for Machine Learning](https://www.cs.waikato.ac.nz/~mhall/thesis.pdf)

この手法を使用すればいくつかのメリットが得られます。

- 互いに相関の高い特徴量の片方を削除することで、精度にあまり影響を与えずに特徴量空間の次元を下げることができる。
- 線形モデルの解釈性を上げることができる。

相関はピアソン相関係数を用いて評価します。

```python
threshold = 0.8

feat_corr = set()
corr_matrix = X_train.corr()
for i in range(len(corr_matrix.columns)):
    for j in range(i):
        if abs(corr_matrix.iloc[i, j]) > threshold:
            feat_name = corr_matrix.columns[i]
            feat_corr.add(feat_name)

print(len(set(feat_corr)))
# -> 658

X_train.drop(labels=feat_corr, axis='columns', inplace=True)
X_test.drop(labels=feat_corr, axis='columns', inplace=True)

print(len(X_train.columns))
# -> 4025
```

### 統計的評価指標を用いる

流れとしては以下のようになります。

1. 特徴量をある評価指標で測る
2. 高いランクの特徴量だけを残す

あくまで1つの特徴量とターゲットを使用して性能を評価するので、特徴量間の関係などは基本的にはわかりません。

```python
# 使用するデータは直前までのFilter Methodを適用させたものが対象です

print(X_train.shape)
# -> (3567, 4025)
print(X_test.shape)
# -> (892, 4025)
```

#### Mutual Infomartion

相関係数を使用して特徴量を削減した手法と似ていますが、より一般的な手法になります。

ある特徴量Xと別の特徴量Yの間の同時分布$P(X, Y)$と個々の分布の積$P(X)P(Y)$がでれだけ似ているのかを計算します。もし互いに独立であれば`Mutual Infomartion`は0になります。

$$
\sum_{i,y} P(xi, yj) \times log \frac{p(xi, yj)}{P(xi)P(yj)}
$$

```python
# 使用するモジュール
from sklearn.feature_selection import mutual_info_regression
from sklearn.feature_selection import SelectKBest, SelectPercentile

MI = mutual_info_regression(X_train, y_train)
MI = pd.Series(MI)
MI.index = X_train.columns
MI.sort_values(ascending=False).plot(kind='bar', figsize=(20,10))

# KBest : 抽出する特徴量の"数"を指定
kbest_sel_ = SelectKBest(mutual_info_regression, k=10)
print(len(kbest_sel_.get_support()))

# Percentile : 抽出する特徴量の割合を指定
percentile_sel_ = SelectPercentile(mutual_info_regression, percentile=10)
print(len(percentile_sel_.get_support()))
```

![mutual_info_plot.PNG](https://qiita-image-store.s3.amazonaws.com/0/194085/9f1a3a33-cf57-c328-fe5c-a723255b529b.png)


#### カイ2乗、フィッシャー係数

これは特徴量がカテゴリ変数でありターゲットが2値の場合によく使用される手法です。

例えば2値の分類タスクであるカテゴリの特徴量を評価する際、カテゴリごとの2値の出現回数などをカテゴリごとに比較したり、全体の出現回数などと比較したりします。

$$
S_{i} = \frac{\sum_{k=1}^{K} n_{j}(\mu_{ij}-\mu_{i]})}{\sum_{k=1}^{K} n_{j} \rho_{ij}^2}
$$


```python
# Santandarのデータには適用できないので、簡単な紹介のみ

from sklearn.feature_selection import chi2

# fisher score
fscore = chi2(X_train.fillna(0), y_train)
```

#### Univariate

ターゲットが2値である必要があり、よく連続値の特徴量に対して使用されます。2変数に対して行う手法は**ANOVA**などと呼ばれます。

前提として（1）変数とターゲットに線形の関係があること、（2）変数が正規分布している必要があります。

```python
# Santandarのデータには適用できないので、簡単な紹介のみ

from sklearn.feature_selection import f_regression
from sklearn.feature_selection import SelectKBest, SelectPercentile

univariate = f_regression(X_train, y_train)
sel_ = SelectKBest(f_regression, k=1000).fit(X_train, y_train)
```

#### Univariate ROC_AUC / RMSE

2変数間の依存性を評価しますが、この手法では機械学習モデルを使用します。

**Univariate**とは異なり、変数に特に前提条件は必要ありません。

流れとしては以下のようになります。

1. 1変数のみを使用して学習させる
2. すべての変数に対して計算にランキング付けする
3. 高いランクの変数だけを選択する

```python
from sklearn.tree import DecesionTreeRegressor
from sklearn.metrics import mean_squared_error, roc_auc_score

MSE_features = []
for feature in X_train.columns:
    clf = DecisionTreeRegressor()
    clf.fit(X_train[feature].to_frame(), y_train)
    y_pred = clf.predict(X_test[feature].to_frame())
    MSE_features.append(mean_squared_error(y_test, y_pred))

MSE_series = pd.Series(MSE_features)
MSE_series.index = X_train.columns
MSE_series.sort_values(ascending=False).plot(kind='bar', figsize=(20, 10))
```

![univariate_mse.PNG](https://qiita-image-store.s3.amazonaws.com/0/194085/ba1ca579-a474-227d-abab-c107547d437b.png)


平均2乗誤差なので低いほうから選択して選ぶ必要があります。

## Wrapper Methodでの解析例

使用するデータは簡単にFilter Methodを適用させたものを対象とします。

```python
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.feature_selection import VarianceThreshold

from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.metrics import roc_auc_score
from mlxtend.feature_selection import SequentialFeatureSelector as SFS

data = pd.read_csv('train.csv')
X_train, X_test, y_train, y_test = train_test_split(
    data.drop(['ID', 'target'], axis=1),
    data['target'],
    test_size=0.2,
    random_state=2018
)

# 1. delete const, quisi-const, duplicated 
# 2. correlated 0.8
```

評価指標としては、例として分類問題なら`ROC_AUC`を選択し、回帰問題なら`R squared`など選択します。

### Step Forward

1. すべての特徴量から1つずつ独立して学習させる。
2. 最も評価指標が高くなったものを選択する
3. [1]で選択した特徴量に、ほかに特徴量をどれか1つのみ追加して学習させる。
4. 最も評価指標が高くなった組み合わせを選択する。
5. 上記を特徴量を1つずつ加えながら繰り返す。

```python
sfs1 = SFS(RandomForestClassifier(n_jobs=4),
           k_features=10,   # いくつの特徴量を出力させたいのか
           forward=True,
           floating=False,
           verbose=2,
           scoring='roc_auc', # 'r2' なども選択できる。
           cv=3)

sfs1 = sfs1.fit(np.array(X_train.fillna(0)), y_train)

selected_feat = X_train.columns[list(sfs1.k_feature_idx_)]
```

最後に抽出された特徴量を使用して学習を行ってみます。

```python
clf = RandomForestClassifier(n_estimators=200, random_state=2018, max_depth=4)
clf.fit(X_train[selected_feat], y_train)
print(roc_auc_score(y_train, clf.predict_proba(X_train[selected_feat])[:, 1]))
print(roc_auc_score(y_test, clf.predict_proba(X_test[selected_feat])[:, 1]))
```

### Step backward

1. すべての特徴量を使用して評価する。
2. 特徴量を1つずつ削除して評価する。
3. 特徴量をさらに1つ削除して評価する。
4. 上記を繰り返す。

```python
sfs1 = SFS(RandomForestClassifier(n_jobs=4),
           k_features=15,   # いくつの特徴量を出力させたいのか
           forward=False,   # False -> Step Backward
           floating=False,
           verbose=2,
           scoring='roc_auc', # 'r2' なども選択できる。
           cv=3)

sfs1 = sfs1.fit(np.array(X_train.fillna(0)), y_train)

selected_feat = X_train.columns[list(sfs1.k_feature_idx_)]
```

最後に評価を行います。

```python
clf = RandomForestClassifier(n_estimators=200, random_state=2018, max_depth=4)
clf.fit(X_train[selected_feat], y_train)
print(roc_auc_score(y_train, clf.predict_proba(X_train[selected_feat])[:, 1]))
print(roc_auc_score(y_test, clf.predict_proba(X_test[selected_feat])[:, 1]))
```

### Exhaustive Search

- すべての組み合わせを試して最も性能が高くなった組み合わせを選択する。

これは計算量が非常に多くなります。例えば4つの特徴量を考えれば、全部で15通りを計算する必要があり、特徴量が多いデータセットには適用させることは難しいです。

```python
from mlxtend.feature_selection import ExhaustiveFeatureSelector as EFS

efs1 = EFS(RandomForestClassifier(n_jobs=4, random_state=2018,
           min_features=1,
           max_features=4,
           scoring='roc_auc',
           print_progress=True,
           cv=2))

efs1 = efs1.fit(np.array(X_train[X_train.columns[0:4]].fillna(0)), y_train)

clf = RandomForestClassifier(n_estimators=200, random_state=2018, max_depth=4)
clf.fit(X_train[selected_feat], y_train)
print(roc_auc_score(y_train, clf.predict_proba(X_train[selected_feat])[:, 1]))
print(roc_auc_score(y_test, clf.predict_proba(X_test[selected_feat])[:, 1]))
```

## Embedded Method での解析例

使用するデータは簡単なFilter Methodを適用させたものにします。

### Lasso正則化

`L1`正則化には重要でない特徴量に対して、その係数をゼロにしていく性質が存在します。そこでこの性質を利用して特徴量を削減していくことを考えます。

解析の際は線形モデルなので特徴量を正規化させておきます。

```python
from sklearn.linear_model import Lasso
from sklearn.feature_selection import SelectFromModel
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaler.fit(X_train)

sel_ = SelectFromModel(Lasso(alpha=100))
sel_.fit(scaler.transform(X_train, y_train))

print(len(sel_.get_support()))
```

### 木モデル

```python
from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import SelectFromModel

sel_ = SelectFromModel(RandomForestRegressor(n_estimators=200))
sel_.fit(X_train, y_train)

rest_features = X_train.columns[sel_.get_support()]
```

## 参考文献

- [A survey on Feature Selection Methods](http://newiranians.ir/Chandrashekar%20-%20Feature%20Selection%20Methods%20-%202014.pdf)
- [Feature-Selection-for-Machine-Learning](https://github.com/anujdutt9/Feature-Selection-for-Machine-Learning)
- [The 2009 Knowledge Discovery in Data Competition (KDD Cup 2009)](http://www.mtome.com/Publications/CiML/CiML-v3-book.pdf)
- [Feature Selection: A literature Review](https://www.cc.gatech.edu/~hic/CS7616/Papers/Kumar-Minz-2014.pdf)
- [An Introduction to Variable and Feature Selection](http://www.jmlr.org/papers/volume3/guyon03a/guyon03a.pdf)
- [Feature Selection for High-Dimensional Data: A Fast Correlation-Based Filter Solution](https://www.aaai.org/Papers/ICML/2003/ICML03-111.pdf)
- [Data Preprocessing for Supervised Learning
](https://www.researchgate.net/publication/228084519_Data_Preprocessing_for_Supervised_Learning)
- [Correlation-based Feature Selection for Machine Learning](https://www.cs.waikato.ac.nz/~mhall/thesis.pdf)
- [A review of feature selection methods with applications](https://ieeexplore.ieee.org/document/7160458)
- [Feature Selection for Classification: A Review](https://pdfs.semanticscholar.org/310e/a531640728702fce6c743c1dd680a23d2ef4.pdf)

ほかにも**Hybrid Method**や**Recursive**なやり方などいろいろあるけど力尽きました。