【機械学習】カイ二乗検定を用いたカテゴリ変数の特徴量選択

# はじめに

この間統計検定２級を取得したのですが、その試験対策の一環で「カイ二乗検定」という統計学的手法について勉強しました。
この「カイ二乗検定」は機械学習の特徴量選択に有用なんじゃないかと思い調べたところ、scikit-learnのSelectKBestで既に実装されていたため、その使い方をメモしたいと思います。


# カイ二乗検定とは

別名「独立性の検定」と呼ばれ、事象Aと事象Bが「独立」であることを検定するためのものです。

事象Aと事象Bが独立と仮定した場合、実測値がどれだけ「ありえない値」となっているかを確認し、「ありえない」＝「事象Aと事象Bは独立ではなく、何かしらの関係がある」ということを検定するもの、というのが私の理解です。

もう少し詳しい説明は以下のサイトが大変参考になります。

> [独立性の検定―最もポピュラーなカイ二乗検定(統計WEB)](https://bellcurve.jp/statistics/blog/14038.html)

# Pythonによる実装例と正解率の比較

## 今回使うデータ

説明変数にカテゴリデータを含む以下のデータを使用します。
このデータは、ある人の年収が50Kを超えているかどうかをその人の年齢・職業・性別などの情報と共に記録されています。

|||
|--------|-------------|
|データ名|Adult Data Set|
|URL|https://archive.ics.uci.edu/ml/datasets/adult|
|データ数| 32561|

上記データを読み込むコードを以下に記録しておきます。
なお、カイ二乗検定はある特徴量の度数が期待値とどれだけ外れているかを検定するためのものなので、数値データ(量的データ)は除外し、カテゴリデータのみを用いて学習します。

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

#データの読み込みと列名指定
#データの読み込みとデータ処理
columns=['age','workclass','fnlwgt','education','education-num',
 'marital-status','occupation','relationship','race','sex','capital-gain',
 'capital-loss','hours-per-week','native-country','Income']
data = pd.read_csv('adult.data.csv', header=None).sample(frac=1).dropna()
data.columns=columns
print(data.shape)
data = data.replace({' <=50K':0,' >50K':1})
data.head()


# 説明変数(カテゴリデータ)と目的変数で分離
data_x = data[['workclass','education','marital-status','occupation','relationship','race','sex','native-country']]
data_y = data['Income']

# カテゴリデータはダミー変数に
data_dm = pd.get_dummies(data_x)

# 学習データとテストデータに分割
X_train, X_test, Y_train, Y_test = train_test_split(data_dm, data_y, train_size=0.7)
```

## 特徴量選択をしなかった場合

まずは特徴量選択をしなかった場合、つまり、全ての特徴量を使ってモデルに学習させた場合の結果を確認してみます。

今回使うモデルは勾配ブースティングで、パラメータチューニングは行わずにデフォルトで学習させます。

```python
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import GradientBoostingClassifier
from imblearn.pipeline import Pipeline 

scaler = StandardScaler()
clf = GradientBoostingClassifier()

# 一連の処理(Pipeline)を定義
# 標準化→識別器訓練
estimator = [
    ('scaler',scaler),
    ('clf',clf)
]

pipe = Pipeline(estimator)
pipe.fit(X_train, Y_train)
print('正解率: ' + str(pipe.score(X_test, Y_test)))
```

>正解率: 0.828743986078

パラメータチューニングも何もしていないので、まぁこんなものかという精度です。
次に、カイ二乗検定を用いた特徴量選択を行なった場合の結果を見てみましょう。

## カイ二乗検定を用いた特徴量選択

カイ二乗検定による特徴量選択は、scikit-learnの「[SelectKBest](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.SelectKBest.html)」を利用すると簡単に実装できます。

このクラスはパラメータscore_funcで指定した評価指標に基づき、特徴量の数を絞ってくれます。使うメソッドはいつものfit・transformです。絞った後の特徴量の数はパラメータkで指定できます。

選択できる評価指標はいくつかありますが、ここではカイ二乗検定を用いるため、「chi2」を指定します。

最初に全ての特徴量で学習させた後、特徴量を一つずつ減らしていき、モデルの正解率が下がり出したタイミングでイテレーションを終了する手法を取っています。（「Backward feature Elimination」という手法らしいです）

```python
from sklearn.feature_selection import SelectKBest, chi2
max_score = 0
for k in range(len(data_dm.columns)):
    print(len(data_dm.columns)-k)
    select = SelectKBest(score_func=chi2, k=len(data_dm.columns)-k) 
    scaler = StandardScaler()
    clf = GradientBoostingClassifier()

    # 一連の処理(Pipeline)を定義
    # 特徴量抽出→標準化→識別器訓練
    estimator = [
        ('select',select),
        ('scaler',scaler),
        ('clf',clf)
    ]

    pipe_select = Pipeline(estimator)
    pipe_select.fit(X_train, Y_train)
    score = pipe_select.score(X_test, Y_test)
    if score < max_score:
        break
    else:
        max_score = score
        pipe_fix = pipe_select

print('正解率: ' + str(pipe_fix.score(X_test, Y_test)))
```
> 正解率: 0.828948715324

少しですが正解率が上がりました。元々のデータには不要な特徴量が含まれていたようで、カイ二乗検定による特徴量選択により、その特徴量を除外できているようです。

ちなみに、除外された特徴量は以下で確認できます。

```python
mask = -pipe_fix.steps[0][1].get_support()
data_dm.iloc[:,mask].columns
```
> Index(['native-country_ Greece', 'native-country_ Holand-Netherlands','native-country_ Thailand'],dtype='object')


## モデルベースの特徴量選択

比較のため、別の特徴量選択の手法も紹介しておきます。
機械学習で用いるモデルには、どの特徴量が精度に大きく寄与しているかを保持するものがあります。例えばLasso回帰, Ridge回帰などです。ここでは、RandomForestによる特徴量選択での実装例とその結果も見てみましょう。

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

select = SelectFromModel(RandomForestClassifier(n_estimators=100, n_jobs=-1))
scaler = StandardScaler()
clf = GradientBoostingClassifier()

# 一連の処理(Pipeline)を定義
# 特徴量抽出→標準化→識別器訓練
estimator = [
    ('select',select),
    ('scaler',scaler),
    ('clf',clf)
]

pipe_rf = Pipeline(estimator)
pipe_rf.fit(X_train, Y_train)
print('正解率: ' + str(pipe_rf.score(X_test, Y_test)))
```
> 正解率: 0.826287235132

全て特徴量を使った時より正解率が落ちています。

特徴量選択で利用しているモデルのパラメータチューニングを行なっていないためにイマイチな精度となっているのではないかと思っています。本題ではないため、本記事では深掘りしません。

# まとめ

本記事では、カイ二乗検定の簡単な説明と、SelectKBestを用いた実装例を紹介しました。
カイ二乗検定による特徴量選択は、特徴量がカテゴリデータでかつ分類問題のみしか有効ではない、というのが今の私の理解です。
特徴量に数値データを含んでいたり、回帰問題の場合、ピアソンの相関係数やANOVAといった評価指標による特徴量選択を用いましょう。
（ちなみに私は、カテゴリデータ・数値データ問わず使えて楽なランダムフォレストによる特徴量選択をよく利用しています。）

# 最後に
記載に何か気になる点・不明点などあれば、コメントください。


# [付録1]SelectKBestにfitする前にダミー変数化している理由

SelectKBestのscore_funcで指定しているchi2のドキュメント([sklearn.feature_selection.chi2](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.chi2.html#sklearn.feature_selection.chi2))を読むと、以下のような記載があります。

> X : {array-like, sparse matrix}, shape = (n_samples, n_features_in)

パラメータとして渡すデータはスパースな配列、つまりほとんどの要素が0の配列と書いています。

ソースを読んでみると以下のようなコードがあり、列ごとの合計値、つまり、その特徴量の度数を計算していることがわかります。この計算を行うために、ダミー変数化する前の文字列を含むデータをパラメータとして渡すとエラーが発生することがわかります。

```python
feature_count = X.sum(axis=0).reshape(1, -1)
```

また、特徴量に数値データ(量的データ)が含まれていても、エラーは発生しませんが謎の計算が行われることになります。
(例えば年齢という特徴量があった場合、全行の年齢の合計を計算し、それを度数として認識してしまいます）

# [付録2]SelectKBestで選択された特徴量の確認

fitした後のSelectKBestインスタンスは、get_support()メソッドでどの特徴量を選択したかのリストを持ちます。
transformメソッドを実行すると、入力データの特徴量のうち、maskがTrueとなっている特徴量を選択して返却する、という動作をしているようです。

以下のソースは、特徴量を10個に絞った場合の実装例です。

```python
from sklearn.feature_selection import SelectKBest, chi2
select = SelectKBest(score_func=chi2, k=10) 
select.fit(X_train, Y_train)

mask = select.get_support()
mask
```
>array([ True,  True,  True,  True,  True,  True, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False,  True, False,  True,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False,  True,
       False, False,  True, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False], dtype=bool)

どの特徴量を選択したかは、例えば以下のようなコードで確認できます。

```python
data_dm.iloc[:,mask].columns
```
> Index(['age', 'fnlwgt', 'education-num', 'capital-gain', 'capital-loss','hours-per-week', 'marital-status_ Married-civ-spouse','marital-status_ Never-married', 'relationship_ Husband','relationship_ Own-child'],dtype='object')