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

# はじめに 
「仕事で始める機械学習」で、線形回帰係数や、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つの目的変数から構成されており、これを用いてさまざまな分析手法を適用します。まず、データを可視化し、基本的な統計情報を調べます。

```python
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
"""
```

## データの整形

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

```python
# データの整形
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
        ) 
```

## 相関関係の可視化

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

```python
# 相関関係の確認
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/85196c58-be21-bdb6-1a3b-49a60372c148.png)

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

## 線形回帰

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

```python
# 線形回帰
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値でソートをかけます. 

```python 
# 線形回帰続き　
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` と特徴量の関係を視覚化します。

```python
# 決定木
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/e3298000-c301-91f9-ddcf-ae076333fee5.png)

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

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


## ランダムフォレスト

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

```python
# ランダムフォレスト 特徴量重要度
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）を用いて、各特徴量が予測にどれだけ寄与しているかを可視化します。

```python
# ランダムフォレスト 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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/de766356-b983-3f3d-3150-ddc112ae0676.png)

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