0. はじめに

機械学習を行う際、いきなりモデルを構築するのではなく、大まかにデータの中身を確認することが多いです。

pandas-profilingを用いた探索的データ解析(EDA)なども非常に有用ですが、特徴量が多すぎると確認が面倒です。

私は普段、以下の流れで特徴量抽出を行い、散布図を作成して中身を確認しております。

フィルター法（同じ値の割合、標準偏差、相関係数）
ラッパー法（RFEによる特徴量選択）
散布図作成
簡単な回帰分析+SHAP値の確認
簡単な逆解析

今回は上記の解析で用いているPythonコードをご紹介します。

1. 参考文献

2. インポート

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import category_encoders as ce
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import RFE

3. フィルター法

まずは設定した基準に基づいて大まかに特徴量を選択します。

① 分散がゼロのカラムを除去

drop_std_zero.py

def drop_std_zero(X):
    """
    分散が0のカラムを除去

    Parameters
    ----------
    X: DataFrame, 説明変数
    
    Returns
    -------
    X_drop_zero_std: DataFarme, 処理結果
    """
    ## 削除するカラムの名前を格納したリスト
    std_zero_variable_cols = X.columns[X.std() < 1e-10]
    
    ## 削除
    X_drop_zero_std = X.drop(std_zero_variable_cols, axis=1)
    return X_drop_zero_std

② 同じ値を持つ割合が閾値よりも大きいものを除外

drop_same_value.py

def drop_same_value(X, threshold=0.9):
    """
    同じ値を持つ割合が閾値よりも大きいものを除去

    Parameters
    ----------
    X: DataFrame, 説明変数
    threshold: float, 閾値

    Returns
    -------
    X_drop_same: DataFrame, 処理結果
    """
    ## 削除するカラムの名前を格納するリスト
    X_same_drop_cols = []

    ## カラムを順に読み込んで処理
    for col in range(X.shape[1]):
        same_value_num = X.iloc[:, col].value_counts()
        if same_value_num.iloc[0] / X.shape[0] >= threshold:
            X_same_drop_cols.append(col)
        X_drop_same = X.drop(X.columns[X_same_drop_cols], axis=1)
    return X_drop_same

③ 相関係数が閾値よりも大きいものを除外

drop_high_corr.py

def drop_high_corr(X, threshold=0.95):
    """
    ピアソン相関係数が閾値よりも大きいカラムのうち片方を除外

    Parameters
    ----------
    X: DataFrame, 説明変数
    threhold: float, 閾値

    Returns
    -------
    X_drop_high_corr: DataFrame, 処理結果
    """
    ## ピアソン相関係数を取得
    X_corr = X.corr()

    ## 削除するインデックスを格納するリスト
    X_corr_drop_cols = []

    ## 閾値を超える要素がなくなるまで削除
    for i in range(X_corr.shape[0]):
        # 相関係数の最大値
        corr_max = max(list(X_corr.max())

        if corr_max >= threshold:
            # 最大値を示す変数のリストを取得
            max_idx_list = X_corr[X_corr==corr_max].dropna(how='all').index
            # 選ばれた2つのうち、相関係数が小さいほうを選択
            corr_sum_1 = X_corr.loc[:, max_idx_list[0]].sum()
            corr_sum_2 = X_corr.loc[:, max_idx_list[1]].sum()

            if corr_sum_1 >= corr_sum_2:
                X_corr_drop_cols.append(max_idx_list[0])
            else:
                X_corr_drop_cols.append(max_idx_list[1])

            # 削除した特徴量の相関係数を0に変更(直近で追加した分を処理)
            X_corr.loc[:, X_corr_drop_cols[-1]] = 0
            X_corr.loc[X_corr_drop_cols[-1], :] = 0
        else:
            break
    X_drop_high_corr = X.drop(X_corr_drop_cols, axis=1)
    return X_drop_high_corr

最後に、上記全てをまとめて実行する関数を用意して使いまわします。

drop_auto.py

def drop_auto(X, same_threshold=0.9, corr_threshold=0.95):
    """
    drop_std_zero, drop_same_value, drop_high_corrをまとめて実行

    Parameters
    ----------
    X: DataFrame, 説明変数
    same_threshold: float, 閾値(drop_same_value)
    corr_threshold: float, 閾値(drop_high_corr)

    Returns
    -------
    X_filtered: DataFrame, 処理結果
    """
    ## 実行
    X_filtered = drop_std_zero(X)
    X_filtered = drop_same_value(X_filterd, same_threshold)
    X_filtered = drop_high_corr(X_filtered, corr_threshold)
    
    ## 削除したカラム名を表示
    print(f'{list(set(X.columns)-set(X_filtered.columns))}を削除しました')
    return X_filtered

4. ラッパー法

次に、擬似的なモデル作成結果に基づき重要度の高いと思われる特徴量を選択します。
方法論はいくつかございますが、ここではRecursive Feature Elimination(RFE)を用います。
詳細はこちらのサイトなどをご参照ください。

select_feature_by_rfe_rfr.py

def select_feature_by_rfe_rfr(X, y, num=5):
    """
    RFE, RandomForestRegressorを用いて特徴量選択

    Parameters
    ----------
    X: DataFrame, 説明変数
    y: ndarray or Series or DataFrame, 目的変数(1次元)
    num: int, 選択する特徴量の数

    Returns
    -------
    X_rfe_selected: DataFrame, 処理結果
    """
    ## モデルの準備
    rfr = RandomForestRegressor(random_state=0)
    selector = RFE(rfr, n_feature_to_select=5).fit(X, y)

    ## 選択
    X_rfe_selected = selector.transform(X)
    selected_cols = [X.columns[i] for i in range(len(X.columns)) if selector.get_support()[i]]

    ## DataFrameの作成
    X_rfe_selected = pd.DataFrame(X_rfe_selected, columns=selected_cols, index=X.index)
    return X_rfe_selected

5. 散布図作成

seabornのpairplotを使用しても良いのですが、カラム数が多いとどうしても可視性が低下するため、普段は自分なりの関数を用意してザックリ確認しています。

以下の関数では、Xとy(1次元)を指定すると、Xの列を2つずつ選択してyとの相関を可視化していきます。

draw_scatter_per_two_cols.py

def draw_scatter_per_two_cols(X, y, fontsize=8):
    """
    Xとyの相関を散布図で可視化

    Parameters
    ----------
    X: DataFrame, 説明変数
    y: Series, 目的変数(1次元)
    fontsize: float, フォントサイズ

    Returns
    -------
    None.
    """
    ## 文字サイズの指定
    plt.rcParams['font.size'] = fontsize

    ## カラムを2つずつ選択して散布図を作成
    for i in range(len(X.columns)//2 + 1):
        if len(X.columns) > 2*i+1:
            fig = plt.figure(figsize=(10, 4))
            plt.subplots_adjust(wspace=0.4) # 図の間のスペースを調整

            ax1 = fig.add_subplot(121)
            sns.scatterplot(data=X, x=X.columns[2*i], y=y, ax=ax1) # 散布図の作成

            ax2 = fig.add_subplot(122)
            sns.scatterplot(data=X, x=X.columns[2*i+1], y=y, ax=ax2) # 散布図の作成

            plt.show() # 描画

        elif len(X.columns) > 2 * i :
            fig = plt.figure(figsize=(10, 4))
            plt.subplots_adjust(wspace=0.4) # 図の間のスペースを調整

            ax1 = fig.add_subplot(121)
            sns.scatterplot(data=X, x=X.columns[2*i], y=y, ax=ax1) # 散布図の作成

            plt.show() # 描画

6. 簡単な回帰分析 + SHAP

次に大まかな回帰分析を行い、SHAP値でモデル精度への寄与度を確認します。

① 作図用の関数

regression_train_test_scatter.py

def regression_train_test_scatter(ax, y_train, y_test,
                                y_test_pred, y_test_pred):
    """
    実測値と予測値をプロット

    Parameters
    ----------
    ax: matplotlib.axes._axes.Axes
    y_test, y_train: 目的変数
    y_test_pred, y_train_pred: 予測結果

    Returns
    -------
    None.
    """
    ## グラフの最大値、最小値
    y_max = max(y_test.max(), y_train.max())
    y_min = min(y_test.min(), y_train.min())

    ## 対角線
    ax.plot([y_min-0.05 * (y_max-y_min),
            y_max+0.05 * (y_max-y_min)],
            [y_min-0.05 * (y_max-y_min),
            y_max+0.05 * (y_max-y_min)],
            alpha=0.5,
            color='black',
            ls='--')
    
    ## プロット
    ax.scatter(y_train, y_train_pred, alpha=0.5, label='train')
    ax.scatter(y_test, y_test_pred, alpha=0.5, label='test')
    ax.legend()
    ax.set_xlabel('obs')
    ax.set_ylabel('pred')

② train_test_split → RFR → 結果確認 → SHAP

simple_rfr_shap.py

def simple_rfr_shap(X, y, test_size=0.25, random_state=0):
    """
    train_test_split → RFR → 結果確認 → SHAP

    Parameters
    ----------
    X: DataFrame, 説明変数
    y: Series, 目的変数
    test_size: float, テストサイズ
    random_state: int, ランダムシード

    Returns
    -------
    None.
    """
    ## trainとtestに分割
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=test_size,
                                                        random_state=random_state)

    ## モデル作成
    model = RandomForestRegressor(random_state=random_state).fit(X_train, y_train)

    ## 予測
    y_train_pred = model.predict(X_train)
    y_test_pred = model.predict(X_test)

    ## 結果の可視化
    fig = plt.figure(figsize=(5, 5))
    ax = fig.add_subplot(111)
    regression_train_test_scatter(ax,
                                y_train, y_test,
                                y_train_pred, y_test_pred)
    plt.show()

    ## 予測精度の表示
    print(f'train R2: {r2_score(y_train, y_train_pred)}')
    print(f'test R2: {r2_score(y_test, y_test_pred)}')

    ## SHAP
    explainer = shap.TreeExplainer(model=model)
    shap_values = explainer.shap_values(X=X)
    shap.summary_plot(shap_values, X)

7. 簡単な逆解析

最後に、以下の関数を用いて逆解析用のデータを生成し、作成したモデルを用いて逆解析を行います。

make_simple_inverse_analysis_data.py

def make_simple_inverse_analysis_data(X, col, num):
  """
  指定したカラム(col)と目的変数との相関を逆解析で確認するためのデータを生成
  逆解析で使用するデータセットは以下のような構成とする
  - col : linspace(max, min, num)
  - col以外 : mode

  Parameters
  ----------
  X : DataFrame, 説明変数
  col : str, 逆解析を行うカラム名
  num : int, 逆解析に用いるデータ点数

  Returns
  -------
  inverse_X : DataFrame, 逆解析用のデータ
  """
  ## 逆解析で着目するカラムについて処理
  focus_df = pd.DataFrame(np.linspace(X[col].min(), X[col].max(), num), columns=[col])
  focus_df.index = [0]*len(focus_df)

  ## 上記以外のカラムについて処理
  mode_without_focus_df = pd.DataFrame(X.drop(col, axis=1).mode().iloc[0]).T

  ## 結合
  inverse_X = pd.merge(focus_df, mode_without_focus_df, how='outer', left_index=True, right_index=True)
  inverse_X = inverse_X.reset_index(drop=True)

  ## 並びを元に戻す
  inverse_X = inverse_X.reindex(columns=X.columns)
  return inverse_X

8. さいごに

非常に大雑把な方法ではありますが、変数の数が多すぎてデータの全体を眺めるのが面倒な際には重宝しています。
もっと良い方法があるよ、という方は是非教えてください！