LoginSignup
1
1

More than 1 year has passed since last update.

@TSYK11361

機械学習モデルを用いたEDA

Posted at

はじめに

「仕事で始める機械学習」で、線形回帰係数や、shap値などをEDAに用いる手法がありました. 実用的に見えたので、肩慣らしのため手持ちのデータで試してみます.

機械学習を用いたEDA

使用するデータの説明　

今回使用したデータはcausalmlのsynthetic_data()から生成した合成データです. 以下がその概要です.

  1. 各観測値（n個）について、p個の独立変数を持つ行列Xが生成されます。各独立変数は$(0, 1)$の一様分布から生成されます。

  2. 期待値bは、以下の数式に基づいて計算されます：
    $$
    b = \sin(\pi \cdot X_{0} \cdot X_{1}) + 2 \cdot (X_{2} - 0.5)^{2} + X_{3} + 0.5 \cdot X_{4}
    $$

  3. 傾向スコアeは、以下の数式に基づいて計算されます：
    $$
    e = \max(\eta, \min(\sin(\pi \cdot X_{0} \cdot X_{1}), 1 - \eta))
    $$
    ここで、$\eta$は小さな正の値で、eが0または1に極端に近づくことを防ぎます。

  4. 傾向スコアeは、調整項adjによって修正されます：
    $$
    e = \text{expit}(\text{logit}(e) - \text{adj})
    $$

  5. 個々の観測値について、treatmentフラグwがeに基づいて生成されます。wは0または1の値を持ち、以下のように生成されます：
    $$
    w \sim \text{Binomial}(1, e)
    $$

  6. 目的変数yは、以下のように生成されます：
    $$
    y = b + (w - 0.5) \cdot \tau + \sigma \cdot \epsilon
    $$
    ここで、$\tau$は個々の観測値に対する処置効果であり、$\sigma$は誤差項の標準偏差です。$\epsilon$は標準正規分布から生成された誤差です。

データの準備

データセットは、25の特徴量と1つの目的変数から構成されており、これを用いてさまざまな分析手法を適用します。まず、データを可視化し、基本的な統計情報を調べます。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from causalml.dataset.regression import synthetic_data
plt.style.use('fivethirtyeight')

# データの準備
n_features = 25
n_samples = 10000
y, X, w, tau, b, e = synthetic_data(mode=1, n=n_samples, p=n_features, sigma=0.5)

# 特徴量の名前付けとデータフレーム化
feature_names = [f'feature_{i}_informative' if i <= 4 else f'feature_{i}_irrelevant' for i in range(X.shape[1])]
source_df = pd.DataFrame()
source_df[feature_names] = X
source_df['treatment'] = w
source_df['y'] = y 

source_df.dtypes 

"""
feature_0_informative    float64
feature_1_informative    float64
feature_2_informative    float64
feature_3_informative    float64
feature_4_informative    float64
feature_5_irrelevant     float64
feature_6_irrelevant     float64
feature_7_irrelevant     float64
feature_8_irrelevant     float64
feature_9_irrelevant     float64
feature_10_irrelevant    float64
feature_11_irrelevant    float64
feature_12_irrelevant    float64
feature_13_irrelevant    float64
feature_14_irrelevant    float64
feature_15_irrelevant    float64
feature_16_irrelevant    float64
feature_17_irrelevant    float64
feature_18_irrelevant    float64
feature_19_irrelevant    float64
feature_20_irrelevant    float64
feature_21_irrelevant    float64
feature_22_irrelevant    float64
feature_23_irrelevant    float64
feature_24_irrelevant    float64
treatment                  int64
y                        float64
dtype: object
"""

データの整形

データセットからカテゴリカルな特徴量と数値的な特徴量を分け、不要な特徴量を削除させた上で、合体させます. 今回のケースではカテゴリカルな特徴量がなかったので、実質的にはほぼ何も手を加えていません.

# データの整形
single_value_column = source_df.nunique()
single_value_column = source_df.nunique() == 1

categorical_df = source_df.select_dtypes(include=['object'])
if len(categorical_df.columns.values) == 0:
    categorical_df = categorical_df
else:
    categorical_df = pd.get_dummies(categorical_df, drop_first=True)

numerical_df = source_df.select_dtypes(include=['int64', 'float64'])
numerical_df.drop(labels=['y'], axis=1, inplace=True) 

def has_column(df): 
    length = len(df.columns.values)
    return length >= 1

list_df = [
    categorical_df, 
    numerical_df
]
list_has_column = []
for i, df in enumerate(list_df): 
    has_column_true = has_column(df) 
    if has_column_true: 
        list_has_column.append(1)
    elif has_column_true ==False: 
        list_has_column.append(0) 
    else: 
        list_has_column.append(0)

list_valid_df = []
for df, has_column_true in zip(list_df, list_has_column): 
    if has_column_true: 
        list_valid_df.append(df)

if len(list_valid_df) == 0: 
    print("No valid DataFrame")
if len(list_valid_df) == 1: 
    converted_df = list_valid_df[0]
else: 
    converted_df = pd.DataFrame()
    for i, df in enumerate(list_valid_df): 
        converted_df = pd.concat(
            converted_df, 
            df, 
            axis=1
        )

相関関係の可視化

データセット内の特徴量間の相関関係をヒートマップで可視化します。

# 相関関係の確認
plt.figure(figsize=(15, 15))
plt.imshow(converted_df.corr(), interpolation='nearest')
plt.colorbar()
plt.xticks(range(len(converted_df.columns)), converted_df.columns, rotation='vertical')
plt.yticks(range(len(converted_df.columns)), converted_df.columns)
plt.show()

image.png

上図のとおり、treatmentのみ$X_0$と$X_1$と緩やかな相関が見て取れなくもないという結果となりました. 元データの生成過程と照らし合わせて違和感ありません.

線形回帰

線形回帰モデルを用いて目的変数 y と特徴量の関係を調べます。また、特徴量の重要性を示します。

# 線形回帰
import sklearn.preprocessing
import statsmodels.api

scaler = sklearn.preprocessing.MinMaxScaler()
standardization_df = pd.DataFrame(scaler.fit_transform(converted_df), index=converted_df.index, columns=converted_df.columns)
ols_df = standardization_df.copy()
ols_df['const'] = 1
ols_model = statsmodels.api.OLS(y, ols_df)
fit_results = ols_model.fit()
fit_summary = fit_results.summary2()
print(fit_summary) 
"""
                   Results: Ordinary least squares
=====================================================================
Model:                OLS              Adj. R-squared:     0.519     
Dependent Variable:   y                AIC:                16845.6461
Date:                 2023-09-04 06:56 BIC:                17040.3252
No. Observations:     10000            Log-Likelihood:     -8395.8   
Df Model:             26               F-statistic:        415.9     
Df Residuals:         9973             Prob (F-statistic): 0.00      
R-squared:            0.520            Scale:              0.31474   
---------------------------------------------------------------------
                       Coef.  Std.Err.    t    P>|t|   [0.025  0.975]
---------------------------------------------------------------------
feature_0_informative  0.5320   0.0208 25.5264 0.0000  0.4912  0.5729
feature_1_informative  0.5128   0.0210 24.4730 0.0000  0.4718  0.5539
feature_2_informative -0.0197   0.0194 -1.0160 0.3097 -0.0576  0.0183
feature_3_informative  0.9753   0.0193 50.5049 0.0000  0.9374  1.0131
feature_4_informative  0.5099   0.0195 26.2117 0.0000  0.4718  0.5480
feature_5_irrelevant   0.0009   0.0194  0.0480 0.9617 -0.0372  0.0390
feature_6_irrelevant   0.0034   0.0193  0.1780 0.8587 -0.0344  0.0413
feature_7_irrelevant   0.0090   0.0194  0.4625 0.6437 -0.0291  0.0470
feature_8_irrelevant   0.0225   0.0195  1.1512 0.2497 -0.0158  0.0608
feature_9_irrelevant   0.0070   0.0195  0.3601 0.7188 -0.0312  0.0453
feature_10_irrelevant  0.0024   0.0195  0.1211 0.9036 -0.0358  0.0405
feature_11_irrelevant -0.0044   0.0194 -0.2285 0.8192 -0.0425  0.0337
feature_12_irrelevant -0.0069   0.0195 -0.3541 0.7233 -0.0450  0.0313
...
=====================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors
is correctly specified.
"""

見やすくするため以下のコードで係数のp値でソートをかけます.

# 線形回帰続き　
sorted_ols_coef = fit_summary.tables[1].sort_values(
    by='Coef.', key=lambda t:abs(t), ascending=False
)
sorted_ols_coef = sorted_ols_coef[sorted_ols_coef['P>|t|'] < 0.5] 

table = pd.DataFrame(sorted_ols_coef[:10]).to_markdown()
import pyperclip 
pyperclip.copy(table)
Coef. Std.Err. t P>abs(t) [0.025 0.975]
feature_3_informative 1.00143 0.0194547 51.475 0 0.963295 1.03956
treatment 0.700441 0.0129727 53.9935 0 0.675012 0.72587
feature_0_informative 0.556961 0.0209712 26.5584 3.17086e-150 0.515854 0.598069
feature_4_informative 0.527073 0.0192631 27.3617 5.27696e-159 0.489313 0.564833
feature_1_informative 0.526586 0.0210706 24.9916 9.30429e-134 0.485284 0.567889
const -0.139034 0.0492756 -2.82156 0.00478852 -0.235624 -0.042444
feature_8_irrelevant -0.0382736 0.0194372 -1.9691 0.0489698 -0.0763744 -0.000172862
feature_15_irrelevant 0.0302639 0.0194537 1.55569 0.119814 -0.00786927 0.0683971
feature_2_informative -0.028409 0.0193807 -1.46584 0.142723 -0.066399 0.00958101
feature_12_irrelevant 0.0282152 0.0192578 1.46513 0.142917 -0.009534 0.0659644

決定木

決定木モデルを用いて目的変数 y と特徴量の関係を視覚化します。

# 決定木
import sklearn.tree
dt_model = sklearn.tree.DecisionTreeRegressor(max_depth=3, random_state=42)
dt_model.fit(converted_df, y)

plt.figure(figsize=(50, 10))
sklearn.tree.plot_tree(dt_model, feature_names=converted_df.columns, filled=True)

image.png

dtreevizでの可視化は以下の通りです.

# 決定木 続き　
import dtreeviz 
viz = dtreeviz.model(
    dt_model, 
    converted_df, 
    y, 
    target_name='y', 
    feature_names=converted_df.columns
) 
viz.view()

ランダムフォレスト

ランダムフォレストを用いて特徴量の重要性を評価します。

# ランダムフォレスト 特徴量重要度
import sklearn.ensemble
rf_model = sklearn.ensemble.RandomForestRegressor(n_estimators=300, min_samples_leaf=100, max_depth=5, n_jobs=1, random_state=42)
rf_model.fit(converted_df, y)

def check_coef(column_names, coef_list, intercept=None):
    weights = dict(zip(column_names, coef_list))
    if intercept:
        weights['intercept'] = intercept
    df = pd.DataFrame.from_dict(weights, orient='index')
    df.columns = ['coef']
    df.sort_values(by='coef', key=lambda t:abs(t), inplace=True, ascending=False)
    print(df.head(10))
check_coef(converted_df.columns, rf_model.feature_importances_) 
"""
                           coef
treatment              0.588229
feature_3_informative  0.216718
feature_1_informative  0.075873
feature_0_informative  0.071717
feature_4_informative  0.039205
feature_2_informative  0.007162
feature_8_irrelevant   0.000097
feature_13_irrelevant  0.000093
feature_15_irrelevant  0.000091
feature_19_irrelevant  0.000088
"""

SHAP（SHapley Additive exPlanations）を用いて、各特徴量が予測にどれだけ寄与しているかを可視化します。

# ランダムフォレスト SHAP
import shap
rf_model = sklearn.ensemble.RandomForestRegressor(n_estimators=300, max_depth=5, n_jobs=1, random_state=42)
rf_model.fit(converted_df, y)
shap.initjs()
explainer = shap.TreeExplainer(rf_model)
shap_values = explainer.shap_values(converted_df)
shap.summary_plot(shap_values, converted_df)

image.png

まとめ

  • 今回使用した合成データでは、合成データの生成過程と平仄のとれた形で元データを解釈することが可能なEDAとなった.
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1