種々の回帰手法によるダブルクロスバリデーションの実装

ダブルクロスバリデーションによる予測値と各サブモデルの重要度を計算

前回の記事に引き続き、ダブルクロスバリデーションで学習された各サブモデルの特徴量重要度を確認してみるコードです。前回はPLS回帰で実装しましたが、今回は複数の回帰分析手法を切り替えて学習するコードにしてみました。

実装

ライブラリのインポート

! pip install dcekit

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import VarianceThreshold
from sklearn.pipeline import Pipeline
from sklearn.model_selection import KFold, GridSearchCV, cross_validate
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
from sklearn.linear_model import Lasso, Ridge
from sklearn.cross_decomposition import PLSRegression
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from dcekit.variable_selection import search_highly_correlated_variables

データのインポート・前処理

今回はよりMIの実務っぽいデータセットとするために、データセットから80サンプルのみ取り出して分析します。
データ前処理として、フラグメント系の記述子の削除と互いに高い相関関係にある変数の削除をします。

url = "https://datachemeng.com/wp-content/uploads/2017/07/logSdataset1290.csv"

df_in = pd.read_csv(url, index_col=0).sample(80, random_state=123)

x = df_in.iloc[:, 1:]
y = df_in.iloc[:, 0]

keep_index = list(map(lambda x: not x, x.duplicated()))
x = x.loc[keep_index, :]
y = y[keep_index]

selected_features = [i for i in x.columns if "fr_" in i]

x = x.drop("Ipc", axis=1)
x = x.drop(selected_features, axis=1)
x = x.loc[:, x.columns[x.var() != 0]]
x = x.drop(x.columns[search_highly_correlated_variables(x, 0.95)], axis=1)

回帰

まず、使いたい回帰アルゴリズムを指定します。

method = "XGBoost"

クロスバリデーションの設定と、各回帰アルゴリズムのパラメータ探索範囲を指定します。

inner_cv = KFold(n_splits=5, shuffle=True, random_state=0)
outer_cv = KFold(n_splits=5, shuffle=True, random_state=1)

pipe = [("selector", VarianceThreshold())]

if method == "LASSO":
    pipe = pipe + [("scaler", StandardScaler()), ("regressor", Lasso())]
    params = {"regressor__alpha": 10 ** np.arange(-3., 4., 1.),
              "regressor__random_state": [42],
              }

elif method == "Ridge":
    pipe = pipe + [("scaler", StandardScaler()), ("regressor", Ridge())]
    params = {"regressor__alpha": 10 ** np.arange(-3., 4., 1.),
              "regressor__random_state": [42],
              }

elif method == "PLS":
    pipe = pipe + [("scaler", StandardScaler()), ("regressor", PLSRegression())]
    params = {"regressor__n_components": range(1, min(50, x.shape[1]))}

elif method == "RandomForest":
    pipe = pipe + [("regressor", RandomForestRegressor())]
    params = {"regressor__max_features": np.arange(0.1, 1.1, 0.1),
              "regressor__n_jobs": [-1],
              "regressor__random_state": [42],
              }

elif method == "XGBoost":
    pipe = pipe + [("regressor", XGBRegressor())]
    params = {"regressor__objective": ["reg:squarederror"],
              "regressor__learning_rate": [.05],
              "regressor__n_estimators": [500, 1000, 5000],
              "regressor__max_depth": [3, 5, 7],
              "regressor__colsample_bylevel": [.4, .7],
              "regressor__reg:lambda": [.1],
              "regressor__reg:alpha": [.1],
              "regressor__importance_type": ["total_gain"],
              "regressor__verbosity": [1],
              "regressor__n_jobs": [-1],
              "regressor__random_state": [42],
              }

else:
    print(f"{method} is an undefined method.")

勾配ブースティング木では、通常木の数は early_stopping で決めますが、スモールデータの解析では、 eval_set に回せるほどデータがなく、回したとしても非常に少ないデータ数なので、変な木の数で止まる可能性があります。私はデータ数が少ないときは early_stopping を諦めて、Grid Searchで木の数も決めています。(このあたりは完全に独自判断なので、別の考え方がございましたら是非コメントにてご教示お願い致します🙇‍♂️)

さて、pipelineと内側のCVの設定 (HPOのGrid Search) をして、ダブルクロスバリデーションを実行します。

pipe = Pipeline(pipe)
gscv = GridSearchCV(pipe,
                    params,
                    scoring="neg_mean_squared_error",
                    cv=inner_cv,
                    verbose=-1
                    )

res_dcv = cross_validate(gscv,
                         x,
                         y,
                         scoring=[gscv.scoring],
                         cv=outer_cv,
                         verbose=-1,
                         return_train_score=True,
                         return_estimator=True,
                         )

結果の確認

まず、ダブルクロスバリデーションの予測値を計算して、汎化性能を確認してみましょう。

def calc_dcv_prediction(inner_gscvs, X, cv):
    fold = np.zeros(X.shape[0], dtype=int)
    dcv_prediction = np.zeros(X.shape[0], dtype=float)
    for i, (inner_gscv, indecies) in enumerate(zip(inner_gscvs, cv.split(X)), 1):
        model = inner_gscv.best_estimator_
        test_index = indecies[1]
        fold[test_index] = i
        dcv_prediction[test_index] = model.predict(X.iloc[test_index, :]).reshape(-1)

    dcv_prediction = {
        "fold": fold,
        "predicted": dcv_prediction,
        }
    dcv_prediction = pd.DataFrame(dcv_prediction,
                                  index=X.index
                                  )

    return dcv_prediction


pred_y_dcv = calc_dcv_prediction(res_dcv["estimator"], x, outer_cv)

def yy_plot(observed, predicted, fold=None):
    if fold is not None:
        fold = fold.astype("category")

    fig, ax = plt.subplots()
    sns.scatterplot(x=observed, y=predicted, hue=fold, palette="deep", ax=ax)
    ax.set_xlabel("Observed")
    ax.set_ylabel("Predicted")

    axis_min = min(*ax.set_xlim(), *ax.set_ylim())
    axis_max = max(*ax.set_xlim(), *ax.set_ylim())
    ax.set_xlim(axis_min, axis_max)
    ax.set_ylim(axis_min, axis_max)
    ax.plot([axis_min, axis_max], [axis_min, axis_max], "k-")
    if fold is not None:
        plt.legend(loc="upper left", bbox_to_anchor=(1, 1), title=fold.name)
    ax.set_aspect("equal", adjustable="box")

    return fig, ax


def calc_metrics(observed, predicted):
    metrics = {
        "R2": r2_score(observed, predicted),
        "MAE": mean_absolute_error(observed, predicted),
        "RMSE": mean_squared_error(observed, predicted, squared=False),
    }
    metrics = pd.DataFrame(metrics, index=["score"])

    return metrics


yy_plot(y,
        pred_y_dcv["predicted"],
        pred_y_dcv["fold"]
        )
plt.show()
display(calc_metrics(y, pred_y_dcv["predicted"]))

$R^2$	MAE	RMSE
0.80	0.64	0.81

前回の記事で紹介したように、OOF予測値と実測値を突き合わせて評価指標を計算しています。yy-plotと評価指標から、サンプルサイズの割にはそれなりの汎化性能は担保できていそうです。

続いて、各foldで学習されたXGBoostモデルの特徴量重要度を確認してみましょう。

def get_importance_of_each_fold(inner_gscvs, X):
    importances_ = []
    for i, inner_gscv in enumerate(inner_gscvs):
        best_pipe = inner_gscv.best_estimator_
        best_steps = dict(best_pipe.steps)
        try:
            if hasattr(best_steps["regressor"], "coef_"):
                importance = best_steps["regressor"].coef_.reshape(-1)
            elif hasattr(best_steps["regressor"], "feature_importances_"):
                importance = best_steps["regressor"].feature_importances_.reshape(-1)
        except:
            print("Algorithm that cannot compute feature importance.")
        feature_names = (best_steps["selector"]
                         .feature_names_in_[best_steps["selector"].get_support()]
                         .reshape(-1)
                         )
        importance = pd.DataFrame(importance,
                                  index=feature_names,
                                  columns=[f"fold{i + 1}"]
                                  )
        importances_.append(importance)

    importances_ = pd.concat(importances_, axis=1, join="outer")

    return importances_


def plot_cv_feature_importance(cv_feature_importance, sns_plotting_method, max_show=20):
    fold_names = cv_feature_importance.columns
    cv_feature_importance.index.name = "feature_name"
    cv_feature_importance = cv_feature_importance.reset_index()

    _plotting_data = cv_feature_importance.melt(id_vars=["feature_name"],
                                                value_vars=fold_names,
                                                var_name="fold",
                                                value_name="importance"
                                                ).dropna(axis=0)

    _order = _plotting_data.groupby("feature_name").mean().abs().sort_values("importance", ascending=False).index
    if len(_order) > max_show:
        _order = _order[:max_show]

    plt.figure(figsize=(6.4, max(0.25 * len(_order), 4.8)))
    plot_args = {"data"     : _plotting_data,
                 "x"        : "importance",
                 "y"        : "feature_name",
                 "order"    : _order,
                 "palette"  : "husl",
                 }
    sns_plotting_method(**plot_args)
    plt.show()


cv_coefs = get_importance_of_each_fold(res_dcv["estimator"], x)
plot_cv_feature_importance(cv_coefs, sns.boxenplot, 20)

各foldで特徴量重要度が安定しており、与えられたデータ全体の傾向を捉えた学習ができていそうです。
この後、SHAPなどのXAI手法でより詳細にモデルの挙動を確認したい場合、全データ学習をしてモデルを再学習することになることが多いと思います。
このとき、全データ学習のモデルと各CVのモデル間で特徴量重要度の傾向が似ていれば、全データ学習のモデルの学習の傾向はダブルクロスバリデーション時の学習の傾向と似ていると考えられ、XAIによる解釈もそれなりに妥当なものになっているのではないか？と考えています。