特徴量削減手法について、大体何かしらを漏らしてしまうことが多いので、一覧をメモ。

基本的に

1. Filter法で大まかに削る
2. Wrapper法もしくはEmbedded法によって削減する
   • 線形のモデルを用いるべきか非線形を用いるべきか

と言った流れで行うことが多い。

もっと網羅的に手法については以下が詳しい。

• 特徴量選択のまとめ - Qiita

Pythonでの実装

まずサンプルデータセットを準備。

import os
import sys
from urllib import request

import numpy as np
import pandas as pd
import sklearn
import dcekit
import scipy

バージョン

f'Python: {sys.version}'

'Python: 3.8.13 | packaged by conda-forge | (default, Mar 25 2022, 06:05:16) \n[Clang 12.0.1 ]'

for package in (np, pd, sklearn, scipy):
    print(f'{package.__name__}: {package.__version__}')

numpy: 1.23.1
pandas: 1.4.3
sklearn: 1.1.2
scipy: 1.9.0

• dcekit==2.9.1
• boruta==0.3

金子研究室のデータセットを使用。

# 配布ページ: https://datachemeng.com/pythonassignment/
url = 'https://datachemeng.com/wp-content/uploads/2017/07/logSdataset1290.csv'

dirpath_cache = os.path.abspath('./_cache')
if not os.path.isdir(dirpath_cache):
    os.mkdir(dirpath_cache)

fpath_csv_cache = os.path.join(dirpath_cache, os.path.basename(url))
if not os.path.isfile(fpath_csv_cache):
    with request.urlopen(url) as response:
        content = response.read().decode('utf-8-sig')
    with open(fpath_csv_cache, 'w', encoding='utf-8-sig') as f:
        f.write(content)

df_data = pd.read_csv(fpath_csv_cache, index_col=0)
df_data.head()

col_tareget = 'logS'
X = df_data.drop(col_tareget, axis=1)
y = df_data[col_tareget]

Filter法

分散が0の特徴量を削除

これは正規化（標準化）するときの妨げにもなるのでまず最初にやることが多い。

from sklearn.feature_selection import VarianceThreshold

vselector = VarianceThreshold(threshold=0.0).fit(X)
X.shape[1], np.sum(vselector.get_support())

(196, 176)

分散が0の特徴量は今後のためにも削除しておく。

X = X.iloc[:, vselector.get_support()]

相関係数が特定の値より大きいペアのうち片方を削除

金子先生のdcekitにも搭載。dcekitでは、相関係数が高いペアを削除する関数が実装されている。

• dcekit - Github

dcekitでは削除する変数の決め方として、最も絶対値の大きい相関係数をとるペアを抽出し、そのうち、他の変数との相関係数の絶対値の和が大きい方を削除する、といったアルゴリズムになっている。

from dcekit.variable_selection.basic_variable_selection import search_highly_correlated_variables
from contextlib import redirect_stdout

# printを表示しないため。
with redirect_stdout(open(os.devnull, 'w')):
    indices = search_highly_correlated_variables(X, threshold_of_r=0.9)
flag_not_highly_correlated = np.array([i not in set(indices) for i in range(X.shape[1])])
X.shape[1], np.sum(flag_not_highly_correlated)

(176, 139)

Filter法の結果を反映。

X = X.iloc[:, flag_not_highly_correlated]

以降標準化しておいた方が都合が良いことがあるので、先にここでしておく。

X_scaled = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
y_scaled = (y - y.mean()) / y.std(ddof=1)

Wrapper法

Boruta

scikit-learnに準拠した特徴量削減器と似たようなインターフェースで用いることができる。一般にとても効果的で気に入って使っている。

真面目にBorutaのハイパーパラメータを調節することは過学習を引き起こしかねないため、自分はあまりやらない。percを100, 90, 80と試す程度。

実装のときの注意点としては、fitする際、DataFrameを受け付けてくれないこと。

なので、DataFrameを用いたいときは、以下のコードの通り.valuesでnp.ndarrayを取り出すのが良さそう。

from sklearn.ensemble import RandomForestRegressor
from boruta import BorutaPy

boruta_selector = BorutaPy(RandomForestRegressor(), perc=100, n_estimators='auto', verbose=0, random_state=334)
boruta_selector.fit(X.values, y.values)
X.shape[1], np.sum(boruta_selector.support_)

(139, 23)

Embedded法

PLSによるVIP値

VIP (Variable Importance in Projection) という、PLSのモデルから算出できる特徴量の重要どのようなものを基準に削減する。

正直分類がWrapperとどっちかわからなかったのだけれど、fitした結果その中の指標から特徴量を削減することからこちらに分類されると思い、こちらに示した。

これについて詳しく書いているサイトが見つからなかったので念の為式とともにメモしておく。

部分最小二乗回帰（PLS）モデルは、次の式に示すように、元のデータ行列$X$を、多数の潜在指標の負荷量（$L$、ローディング）とスコア（$T$、主成分）の2つの直交行列と残差行列$E$に分解する。

\[X = TL^T + E = \sum_{j=1}^a t_jl_j^T + E\]

ここで、$a$は成分数を表す。

スコア行列$T$は次の式に示すように回帰行列$b$を介して応答ベクトル$D_e$に関係する。$F$は$D_e$の残差ベクトルである。

\[D_e = Tb^T + F\]

各プロセス$m$の指標の値は、プロセスベクトル$x_m$として表される。プロセスベクトル$x_m$は、次の式で与えられるように、重みベクトル$w_j$を介してスコアベクトルに関連付けることができる。

\[t_j = w_j^T x_m\]

特定の指標の投影における変数重要度（The Variable Importance in Projection, VIP）は、回帰係数$b$、重みベクトル$w_j$、スコアベクトル$t_j$を用いて、次の式で与えられるように計算される。

\[\mathrm{VIP}_K = \sqrt{\frac{\displaystyle \sum_{j=1}^a b^2_j t_j^T t_j \frac{w_{kj}}{\|w_j\|}}{\displaystyle \sum_{j=1}^a b_j^2 t_j^T t_j}}\]

ここで、$w_{kj}$は重みベクトル$w_j$の$k$番目の要素である。

PLS-VIP は、集計指標$D_e$に影響を与える各指標の重要性を特定するために使用される。

指標の相対的な変動が結果に影響するのを避けるため、VIP では正規化された指標を使用する。

VIPスコアが低い指標は$D_e$にほとんど影響を与えず、VIPスコアが高い指標は$D_e$に最も貢献する。VIPスコアの2乗の平均は1になる。

VIPスコアが1より大きい場合、指標の相対的な重要性を判断する基準として一般的に使用される。

上記は以下を参考にした。

• Mukherjee, R., Sengupta, D., & Sikdar, S. K. (2015). Selection of Sustainable Processes Using Sustainability Footprint Method: A Case Study of Methanol Production from Carbon Dioxide. In Computer Aided Chemical Engineering (Vol. 36, pp. 311-329). Elsevier. DOI: 10.1016/B978-0-444-63472-6.00012-4

これを実装する。

ライブラリのインポート

from sklearn.cross_decomposition import PLSRegression

関数の実装。

def calc_vip(estimator:PLSRegression)->np.ndarray:
    """Calculate VIP (The Variable Importance in Projection)

    References
    - Mukherjee, R., Sengupta, D., & Sikdar, S. K. (2015). Selection of Sustainable Processes Using Sustainability Footprint Method: A Case Study of Methanol Production from Carbon Dioxide. In Computer Aided Chemical Engineering (Vol. 36, pp. 311-329). Elsevier. DOI: [10.1016/B978-0-444-63472-6.00012-4](https://doi.org/10.1016/B978-0-444-63472-6.00012-4)
    - [Variable Importance in Projection](https://www.sciencedirect.com/topics/engineering/variable-importance-in-projection)
    - [PLSRegression VIP score calculation #7050 - Github](https://github.com/scikit-learn/scikit-learn/issues/7050#issuecomment-345208503)

    Parameters
    ----------
    estimator : PLSRegression
        fitted pls regression model

    Returns
    -------
    np.ndarray
        VIP of each variable (n_features,)
    """
    # score: 潜在変数、主成分
    t = estimator.x_scores_     # (n_samples, n_components)
    w = estimator.x_weights_    # (n_features, n_components)
    # loading: 因子負荷量
    q = estimator.y_loadings_   # (n_targets, n_components)
    p, h = w.shape
    s = np.diag(t.T @ t @ q.T @ q).reshape(h, -1)
    return np.sqrt(p * (np.square(w / np.linalg.norm(w, axis=0)) @ s).ravel() / np.sum(s))

pls = PLSRegression().fit(X_scaled, y_scaled)

閾値を1に設定して、

flag_high_vip = calc_vip(pls) > 1.0
X.shape[1], np.sum(flag_high_vip)

(139, 28)

LASSOの回帰係数

L1正規化に伴い、回帰係数が0になりやすい性質を利用して回帰係数が0になってしまったものを削除する。

今回はLassoCV（交差検証を勝手にして、ハイパーパラメータをよしなにしてくれる）を用いる。

from sklearn.linear_model import LassoCV
from sklearn.feature_selection import SelectFromModel

lasso = LassoCV(random_state=334, n_jobs=-1, cv=5)
# SelectFromModelでは、thresholdと等しい値でも残すことになる。
# したがって、thresholdを0に指定するとすべての特徴量が残ってしまうことになる。
# そこで、scikit-learnのデフォルトと同じ1e-5を使用する。
lasso_selector = SelectFromModel(lasso, threshold=1e-5).fit(X_scaled, y_scaled)
X.shape[1], np.sum(lasso_selector.get_support())

(139, 116)

ElasticNetの回帰係数

Lassoと同様。少し制限が緩くなるイメージ。

from sklearn.linear_model import ElasticNetCV

en = ElasticNetCV(random_state=334, n_jobs=-1, cv=5)
en_selector = SelectFromModel(en, threshold=1e-5).fit(X_scaled, y_scaled)
X.shape[1], np.sum(en_selector.get_support())

(139, 118)

参考

• Mukherjee, R., Sengupta, D., & Sikdar, S. K. (2015). Selection of Sustainable Processes Using Sustainability Footprint Method: A Case Study of Methanol Production from Carbon Dioxide. In Computer Aided Chemical Engineering (Vol. 36, pp. 311-329). Elsevier. DOI: 10.1016/B978-0-444-63472-6.00012-4