はじめに

特徴量選択をシステマチックにできないかなと思っていろいろ調べていたら下記の論文を見つけました。材料の特性予測についての論文で、どういう手順で特徴量選択をしたのかが詳しく書かれています。この論文について、DeepL、ChatGPT3.5、Claude3 Sonnetを駆使して実装していこうと思います。
誤読を含む可能性があります。間違いあればご指摘をお願いします。
https://pubs.aip.org/aip/jcp/article/159/19/194106/2921572/Gradient-boosted-and-statistical-feature-selection
(本論文では、この後、特徴量選択で使用したGBDTと異なるGBDT系のモデルまたはニューラルネットを使用してモデルの学習、ハイパーパラメータの調整を実施しています。今回は、特徴量選択の部分でどうしているのかを見たかったのでそこまでは追わないことにします。)

論文の概要

筆者の上げている問題点は下記のとおりです。

特徴の関連性や相互作用が見落とされてモデルの解釈性が落ちる
モデルの汎化性能を得るために正則化技術が必要になる

これらを解決するために勾配ブースティング(GBDT)と統計的な分析を活用しています。また、多重共線性の削除のため、相関の確認と階層クラスター分析も実施しています。

(Gradient boosted and statistical feature selection workflow for materials property predictions https://pubs.aip.org/aip/jcp/article/159/19/194106/2921572/Gradient-boosted-and-statistical-feature-selection )

実装

実際の論文では材料のデータセットを使用していますが、今回はsklearnからとれる糖尿病の進行具合のデータセットを使いたいと思います。
評価関数は、説明変数の数が大きく動く可能性があるので、自由度調整済み決定係数を用います。

from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
columns = ['age', 'sex', 'bmi', 'bq', 's1', 's2', 's3', 's4', 's5', 's6']
df = pd.DataFrame(columns=columns, data=diabetes.data)
df['target'] = diabetes.target

STEP1：Gradient boosted feature selection

GBDT系のモデルを使って特徴量選択を行います。
手順としては、

あるGBDT系のモデルで学習してfeature_importanceを確認する。
別のGBDT系のモデルで、1.で得たfeature_importanceの重要度が高い順に特徴量を増やしていって精度の改善がなくなったところまでの特徴量を必要な特徴量とする。

です。GBDT系のモデルというだけで特にモデルの指定はなかったように思うので、1.XGBoost、2.LightGBMで実施します。また、1.での学習はパラメータはデフォルトで実施します。
今回のデータセットでは最後まで精度の改善が見られますが、6個でうちきります。
（重要度確認の前に相関の確認・特徴量の削減をしてもいいような気がしています。）

import lightgbm as lgb
import xgboost as xgb

def select_features(train, target, model):
    X = train.drop(target, axis=1)
    y = train[target]
    model.fit(X, y)
    pred = model.predict(X)
    print('Adj-R2:', adjusted_r2(y, pred, X.columns))
    print('R2:', r2_score(y, pred))
    df = pd.DataFrame(columns=['importance'], index=X.columns, data=model.feature_importances_)
    df = df.sort_values('importance', ascending=False)
    return df

lgb_model = lgb.LGBMRegressor(random_state=42)
xgb_model = xgb.XGBRegressor(random_state=42)
train_ = df.copy()
selected_features = select_features(train_, 'target', xgb_model)

def select_n_features(df, train, target):
    y = train[target]
    scores = []
    for i in range(df.shape[0]):
        #print(list(xgbDf.index[:i+1]))
        X = train[list(df.index[:i+1])]
        lgb_model.fit(X, y)
        pred = lgb_model.predict(X)
        #print(adjusted_r2(y,pred,X.columns))
        scores.append(adjusted_r2(y, pred,X.columns))
    return scores

scores = select_n_features(selected_features, train_, 'target')

"""
精度が改善されなくなる場所がどこが確認する
"""
plt.figure(figsize=(10,5))
plt.plot(np.arange(0, len(scores), 1), scores)
plt.xlim(0, len(scores))
plt.grid()
plt.title('target')

# ベースの特徴量の決定
n_features = {'target':6}
features = list(selected_features.index[:n_features['target']])

STEP2：Feature analyses

STEP1(GBFS:Gradient Boosting Feature Selection)で得られた特徴量について統計的な分析を実施しています。このSTEPに関してはおそらくベーシックな分析を実施しています。下記の記事が非常に参考になります。

STEP3：Features engineering

特徴量エンジニアリングと多重共線性の回避をしていきます。特徴量エンジニアリングはSTEP2までで得られている特徴量にから、2変数ずつ取り出して2変数の比率を計算して新たな特徴量とします。比率を計算するので、追加される特徴量の個数は順列で計算できます。今回、STEP2で6つの特徴量を選んでいるので、追加される特徴量の数は$_6P_2 = 30$になります。
（コード中の特徴量生成前の出力は、特徴量を6つに絞る前のものなので11になっています。）

def create_features(train, features, target):
    target_col = train[target]
    train = train[features]
    for col in features:
        for col2 in features:
            if col != col2:
                #print(col, col2)
                col_name = col + '_' + col2
                train[col_name] = train[col] / train[col2]
        #train2.columns
    train[target] = target_col
    return train

print('特徴量生成前：', train_.shape)
train_ = create_features(train_, features, 'target')
print('特徴量生成後：',train_.shape)

"""
出力
特徴量生成前： (442, 11)
特徴量生成後： (442, 43)
"""

STEP4：Multicollinearity reduction

多重共線性への対応は、相関行列と階層クラスタリングの2つを実施しています。

相関行列

STEP3までで得た特徴量の相関行列を確認して相関が0.85以上の特徴量を削除します。

def drop_features(train, target):
    # 相関係数の計算
    train_ = train.drop([target], axis=1)
    corr_matrix_ = train_.corr().abs()
    corr_matrix = train.corr().abs()

    # 相関係数が0.85以上の変数ペアを抽出
    high_corr_vars = np.where(np.triu(corr_matrix_, k=1) > 0.85)
    high_corr_pairs = [(train_.columns[x], train_.columns[y]) for x, y in zip(*high_corr_vars)]

    # 目的変数との相関係数が小さい方の変数を削除
    for pair in high_corr_pairs:
        var1_corr = corr_matrix.loc[target, pair[0]]
        var2_corr = corr_matrix.loc[target, pair[1]]

        try: # 既に消した変数が入ってきたとき用
            if var1_corr < var2_corr:
                train = train.drop(pair[0], axis=1)
            else:
                train = train.drop(pair[1], axis=1)
        except:
            pass
    return train

train_2 = drop_features(train_, 'target')
print(train_2.shape)
"""
出力
(442, 33)
"""

階層クラスタリング

各特徴量をそれ自身のクラスタとみなして、クラスタ分析を実施します。特徴量間の非類似度の指標はスピアマンの順位相関を使い、距離関数はウォード法を使っています。本論文内では、距離が1以下のクラスタ内の特徴量は十分に相関があるとみなし、クラスタ内から1つだけ特徴量を選択します。どれを取り出すかについては言及がなかったため、今回の実装ではクラスタ内の先頭ものを取り出しています。

from scipy.cluster.hierarchy import linkage, fcluster, dendrogram
from scipy.spatial.distance import squareform
from scipy.stats import spearmanr

def remove_collinear_features(train, target, threshold=1.0):
    X = train.drop(target, axis=1)
    y = train[target]
    cols = X.columns
    # 特徴量間の非類似性距離行列を計算
    X_ = X.apply(lambda x: (x - x.mean()) / x.std(), axis=0)  # 標準化
    distances = np.zeros((X_.shape[1], X_.shape[1]))
    for i in range(X_.shape[1]):
        for j in range(i+1, X_.shape[1]):
            corr, _ = spearmanr(X_.iloc[:, i], X_.iloc[:, j])
            distances[i, j] = distances[j, i] = 1 - abs(corr)
    np.fill_diagonal(distances, 0)  # 対角成分をゼロに設定
    distances = squareform(distances)
    
    # Ward の最小分散基準で階層的クラスター分析
    clusters = linkage(distances, method='ward')
    cluster_labels = fcluster(clusters, threshold, criterion='distance')
    # クラスター内で1つの特徴量のみ残す
    unique_cluster_labels = np.unique(cluster_labels)
    unique_features = []
    for label in unique_cluster_labels:
        features = X.columns[cluster_labels == label]
        print(f"同じクラスタの特徴量は{features}です。")
        if len(features) > 1:
            print(f"選ばれたのは{features[0]}でした。")
            unique_features.append(features[0])
        else:
            print(f"選ばれたのは{features}でした。")
            unique_features.extend(features)
    
    df = X[unique_features]
    df[target] = y

    return df, clusters, cols

train_3, clusters, columns = remove_collinear_features(train_2, 'target',threshold=1.0)
print(train_3.shape)

"""
出力
同じクラスタの特徴量はIndex(['s5', 'bmi', 's6', 's1', 'bq', 's3'], dtype='object')です。
選ばれたのはs5でした。
同じクラスタの特徴量はIndex(['s5_s6', 'bmi_s6', 's6_s5', 's6_bmi', 's6_s3', 's3_s6'], dtype='object')です。
選ばれたのはs5_s6でした。
同じクラスタの特徴量はIndex(['s5_s3', 's1_s3', 'bq_s3'], dtype='object')です。
選ばれたのはs5_s3でした。
同じクラスタの特徴量はIndex(['s5_bmi', 'bmi_s5', 'bmi_s3', 's3_bmi'], dtype='object')です。
選ばれたのはs5_bmiでした。
同じクラスタの特徴量はIndex(['s5_bq', 'bmi_bq', 'bq_s5', 'bq_bmi', 'bq_s6'], dtype='object')です。
選ばれたのはs5_bqでした。
同じクラスタの特徴量はIndex(['s1_s5', 's3_s5'], dtype='object')です。
選ばれたのはs1_s5でした。
同じクラスタの特徴量はIndex(['s6_s1', 's1_s6'], dtype='object')です。
選ばれたのはs6_s1でした。
同じクラスタの特徴量はIndex(['bmi_s1', 's1_bmi'], dtype='object')です。
選ばれたのはbmi_s1でした。
同じクラスタの特徴量はIndex(['s1_bq', 'bq_s1'], dtype='object')です。
選ばれたのはs1_bqでした。
(442, 10)
"""

STEP5：The final step of the feature-selection process by RFE

REFCVを使って特徴量をさらに絞ります。特徴量の重要度を決める外部推定期はGBDTを使用して、10分割交差検証のMAEを評価の基準としています。
REFCVについては下記の@FukuharaYoheiさんの記事が詳しかったです。（ほかの特徴量選択についてもまとめられており勉強になりました。）

from sklearn.feature_selection import RFECV
from sklearn.ensemble import GradientBoostingRegressor

def select_features_by_REF(train, target, n_features):
    X = train.drop(target, axis=1)
    y = train[target]
    cv = KFold(10)
    rfe = RFECV(estimator=GradientBoostingRegressor(random_state=0), min_features_to_select=n_features, step=0.5, cv=cv, scoring='neg_mean_absolute_error')
    rfe.fit(X, y)
    train_selected = pd.DataFrame(columns=X.columns.values[rfe.support_], data=train[X.columns.values[rfe.support_]])#data=rfe.transform(X))
    train_selected[target] = y
    return train_selected

train_4 = select_features_by_REF(train_3, 'target', n_features=5)
print(train_4.shape)

"""
出力
(442, 6)
"""

まとめ

最終的に['s5', 's5_s6', 's5_bmi', 's1_s5', 's1_bq']の5つの特徴量が残りました。特徴量の数の変遷としては、10 → 6 → 37 → 32 → 10 → 6 となりました。特徴量選択に行き詰ったときや初手の部分で試してみるのはいいのかなと思いました。

参考記事

下記の論文、記事を参考にさせていただきました。非常に勉強になりました。

Gradient boosted and statistical feature selection workflow for materials property predictions https://pubs.aip.org/aip/jcp/article/159/19/194106/2921572/Gradient-boosted-and-statistical-feature-selection

全コード

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import r2_score
from sklearn.model_selection import cross_val_score, KFold, StratifiedKFold
import warnings
warnings.simplefilter('ignore')

target_list = ['target1', 'target2']

def adjusted_r2(y, pred, columns):
    diff = (y - pred)**2
    diff_ = (y - y.mean())**2
    r = sum(diff) / sum(diff_)
    n_cols = len(columns)
    n_sample = len(y)
    adjusted_r2_score = 1 - ((n_sample-1)/(n_sample-n_cols-1)) * r
    return adjusted_r2_score

from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
columns = ['age', 'sex', 'bmi', 'bq', 's1', 's2', 's3', 's4', 's5', 's6']
df = pd.DataFrame(columns=columns, data=diabetes.data)
df['target'] = diabetes.target
df.head(2)

import lightgbm as lgb
import xgboost as xgb

def select_features(train, target, model):
    X = train.drop(target, axis=1)
    y = train[target]
    model.fit(X, y)
    pred = model.predict(X)
    print('Adj-R2:', adjusted_r2(y, pred, X.columns))
    print('R2:', r2_score(y, pred))
    df = pd.DataFrame(columns=['importance'], index=X.columns, data=model.feature_importances_)
    df = df.sort_values('importance', ascending=False)
    return df

lgb_model = lgb.LGBMRegressor(random_state=42)
xgb_model = xgb.XGBRegressor(random_state=42)
train_ = df.copy()
selected_features = select_features(train_, 'target', xgb_model)

def select_n_features(df, train, target):
    y = train[target]
    scores = []
    for i in range(df.shape[0]):
        #print(list(xgbDf.index[:i+1]))
        X = train[list(df.index[:i+1])]
        lgb_model.fit(X, y)
        pred = lgb_model.predict(X)
        #print(adjusted_r2(y,pred,X.columns))
        scores.append(adjusted_r2(y, pred,X.columns))
    return scores

scores = select_n_features(selected_features, train_, 'target')

plt.figure(figsize=(10,5))
plt.plot(np.arange(0, len(scores), 1), scores)
plt.xlim(0, len(scores))
plt.grid()
plt.xlabel('Number of features')
plt.ylabel('Adj-R2')
plt.title('target')

n_features = {'target':6}
features = list(selected_features.index[:n_features['target']])

def create_features(train, features, target):
    target_col = train[target]
    train = train[features]
    for col in features:
        for col2 in features:
            if col != col2:
                #print(col, col2)
                col_name = col + '_' + col2
                train[col_name] = train[col] / train[col2]
        #train2.columns
    train[target] = target_col
    return train

print('特徴量生成前：', train_.shape)
train_ = create_features(train_, features, 'target')
print('特徴量生成後：',train_.shape)

def drop_features(train, target):
    # 相関係数の計算
    train_ = train.drop([target], axis=1)
    corr_matrix_ = train_.corr().abs()
    corr_matrix = train.corr().abs()

    # 相関係数が0.85以上の変数ペアを抽出
    high_corr_vars = np.where(np.triu(corr_matrix_, k=1) > 0.85)
    high_corr_pairs = [(train_.columns[x], train_.columns[y]) for x, y in zip(*high_corr_vars)]

    # 目的変数との相関係数が小さい方の変数を削除
    for pair in high_corr_pairs:
        var1_corr = corr_matrix.loc[target, pair[0]]
        var2_corr = corr_matrix.loc[target, pair[1]]

        try: # 既に消した変数が入ってきたとき用
            if var1_corr < var2_corr:
                train = train.drop(pair[0], axis=1)
            else:
                train = train.drop(pair[1], axis=1)
        except:
            pass
    return train

train_2 = drop_features(train_, 'target')
print(train_2.shape)

from scipy.cluster.hierarchy import linkage, fcluster, dendrogram
from scipy.spatial.distance import squareform
from scipy.stats import spearmanr

def remove_collinear_features(train, target, threshold=1.0):
    X = train.drop(target, axis=1)
    y = train[target]
    cols = X.columns
    # 特徴量間の非類似性距離行列を計算
    X_ = X.apply(lambda x: (x - x.mean()) / x.std(), axis=0)  # 標準化
    distances = np.zeros((X_.shape[1], X_.shape[1]))
    for i in range(X_.shape[1]):
        for j in range(i+1, X_.shape[1]):
            corr, _ = spearmanr(X_.iloc[:, i], X_.iloc[:, j])
            distances[i, j] = distances[j, i] = 1 - abs(corr)
    np.fill_diagonal(distances, 0)  # 対角成分をゼロに設定
    distances = squareform(distances)
    
    # Ward の最小分散基準で階層的クラスター分析
    clusters = linkage(distances, method='ward')
    cluster_labels = fcluster(clusters, threshold, criterion='distance')
    # クラスター内で1つの特徴量のみ残す
    unique_cluster_labels = np.unique(cluster_labels)
    unique_features = []
    for label in unique_cluster_labels:
        features = X.columns[cluster_labels == label]
        print(f"同じクラスタの特徴量は{features}です。")
        if len(features) > 1:
            print(f"選ばれたのは{features[0]}でした。")
            unique_features.append(features[0])
        else:
            print(f"選ばれたのは{features}でした。")
            unique_features.extend(features)
    
    df = X[unique_features]
    df[target] = y

    return df, clusters, cols

train_3, clusters, columns = remove_collinear_features(train_2, 'target',threshold=1.0)
print(train_3.shape)

fig2, ax2 = plt.subplots(figsize=(20,5))
ax2 = dendrogram(clusters, labels=columns)
fig2.show()

from sklearn.feature_selection import RFECV
from sklearn.ensemble import GradientBoostingRegressor

def select_features_by_REF(train, target, n_features):
    X = train.drop(target, axis=1)
    y = train[target]
    cv = KFold(10)
    rfe = RFECV(estimator=GradientBoostingRegressor(random_state=0), min_features_to_select=n_features, step=0.5, cv=cv, scoring='neg_mean_absolute_error')
    rfe.fit(X, y)
    train_selected = pd.DataFrame(columns=X.columns.values[rfe.support_], data=train[X.columns.values[rfe.support_]])#data=rfe.transform(X))
    train_selected[target] = y
    return train_selected

train_4 = select_features_by_REF(train_3, 'target', n_features=5)
print(train_4.shape)

train_4.columns