Partial Dependence Plots(PDP)

높은 성능의 Modeling을 위해 RandomForest, Boosting과 같은 Ensemble model을 주로 사용하게 된다.
이런 복잡도가 높은 Model은 Linear model에 비해 해석하기 어렵다.
- Complex model : 이해하기 어렵지만 성능이 좋다.
- Simple model : 이해하기 쉽지만 성능이 아쉽다.
예로 RandomForest, Boosting의 경우 쉽게 Feature importance 값을 얻을 수 있는데, 이를 통해서 알 수 있는 것은 어떤 특성들이 Model의 성능에 중요하다, 많이 쓰인다 정도다.
특성의 값에 따라 Target 값이 증가/감소하느냐와 같은 어떻게 영향을 미치는지에 대한 정보는 알 수 없다.
PDP(Partal dependence plot, 부분의존도그림)를 사용하면 특성들이 타겟에 어떻게 영향을 주는지 파악할 수 있다.

import pandas as pd

# Kaggle 데이터세트에서 10% 샘플링된 데이터
## Source: https://www.kaggle.com/wordsforthewise/lending-club
## 10% of expired loans (loan_status: ['Fully Paid' and 'Charged Off'])
## grades A-D
## term ' 36 months'

df = pd.read_csv('https://ds-lecture-data.s3.ap-northeast-2.amazonaws.com/lending_club/lending_club_sampled.csv')
df['issue_d'] = pd.to_datetime(df['issue_d'], infer_datetime_format=True)

# issue_d로 정렬
df = df.set_index('issue_d').sort_index()

df['interest_rate'] = df['int_rate'].astype(float)
df['monthly_debts'] = df['annual_inc'] / 12 * df['dti'] / 100

# 152 특성 중 6특성만 사용
columns = ['annual_inc', # 연수입
           'fico_range_high', # 신용점수 
           'funded_amnt', # 대출
           'title', # 대출 목적
           'monthly_debts', # 월간 부채
           'interest_rate'] # 이자율

df = df[columns]
df = df.dropna()

# 마지막 10,000 대출은 테스트셋
# 테스트셋 전 10,000 대출이 검증셋
# 나머지는 학습셋
test = df[-10000:]
val = df[-20000:-10000]
train = df[:-20000]


df.columns
'''
Index(['annual_inc', 'fico_range_high', 'funded_amnt', 'title',
       'monthly_debts', 'interest_rate'],
      dtype='object')
'''


test.shape, val.shape, train.shape
'''
((10000, 6), (10000, 6), (76408, 6))
'''


# 타겟은 이자율
target = 'interest_rate' 
features = df.columns.drop('interest_rate')

X_train = train[features]
y_train = train[target]
X_val = val[features]
y_val = val[target]
X_test = test[features]
y_test = test[target]


# 타겟이 약간 right skewed 되어 있으나 큰 문제는 아님
%matplotlib inline
import seaborn as sns
sns.displot(y_train, kde=True);

# 선형회귀 학습
from category_encoders import TargetEncoder
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

linear = make_pipeline(
    TargetEncoder(),  
    LinearRegression()
)

linear.fit(X_train, y_train)
print('R^2', linear.score(X_val, y_val))
'''
R^2 0.17585064958162422
'''


# 선형회귀 결과 해석
## 회귀계수
coefficients = linear.named_steps['linearregression'].coef_
pd.Series(coefficients, features)
'''
annual_inc        -0.000005 # 연수입
fico_range_high   -0.052805 # 신용점수
funded_amnt        0.000021
title              1.007214
monthly_debts      0.000019
dtype: float64
'''

-0.000005 * 10000
'''
-0.05
연수입의 이자율에 대한 영향 : $10k 더 벌 수록 0.05 이율이 줄어든다.
'''

-0.052805 * 100
'''
-5.2805
신용점수의 이자율에 대한 영향 : 100 point 오를 때마다 5% 가까이 이율이 줄어든다.
'''


# Gradiant Boosting
from category_encoders import OrdinalEncoder
from sklearn.metrics import r2_score
from xgboost import XGBRegressor

encoder = OrdinalEncoder()
X_train_encoded = encoder.fit_transform(X_train) # 학습데이터
X_val_encoded = encoder.transform(X_val) # 검증데이터

boosting = XGBRegressor(
    n_estimators=1000,
    objective='reg:squarederror', # default
    learning_rate=0.2,
    n_jobs=-1
)

eval_set = [(X_train_encoded, y_train), 
            (X_val_encoded, y_val)]

boosting.fit(X_train_encoded, y_train, 
          eval_set=eval_set,
          early_stopping_rounds=50
         )
'''
[0]	validation_0-rmse:9.44760	validation_1-rmse:10.18662
Multiple eval metrics have been passed: 'validation_1-rmse' will be used for early stopping.

Will train until validation_1-rmse hasn't improved in 50 rounds.
[1]	validation_0-rmse:7.74920	validation_1-rmse:8.52176
[2]	validation_0-rmse:6.42741	validation_1-rmse:7.23395
[3]	validation_0-rmse:5.41177	validation_1-rmse:6.24461
.
.
.
[100]	validation_0-rmse:2.53256	validation_1-rmse:3.27319
[101]	validation_0-rmse:2.53215	validation_1-rmse:3.27314
[102]	validation_0-rmse:2.53171	validation_1-rmse:3.27313
Stopping. Best iteration:
[52]	validation_0-rmse:2.59294	validation_1-rmse:3.26819

XGBRegressor(base_score=0.5, booster=None, colsample_bylevel=1,
             colsample_bynode=1, colsample_bytree=1, gamma=0, gpu_id=-1,
             importance_type='gain', interaction_constraints=None,
             learning_rate=0.2, max_delta_step=0, max_depth=6,
             min_child_weight=1, missing=nan, monotone_constraints=None,
             n_estimators=1000, n_jobs=-1, num_parallel_tree=1, random_state=0,
             reg_alpha=0, reg_lambda=1, scale_pos_weight=1, subsample=1,
             tree_method=None, validate_parameters=False, verbosity=None)
'''


y_pred = boosting.predict(X_val_encoded)
print('R^2', r2_score(y_val, y_pred))
'''
R^2 0.22947518106177134
'''

한 특성 사용

선형모델은 회귀계수를 이용해 변수와 타겟 사이의 관계를 해석할 수 있지만 트리모델은 할 수 없다.
대신 PDP를 사용해 개별 특성과 타겟 간의 관계를 볼 수 있다.

# dpi(dots per inch) 수치를 조정해 이미지 화질 조정
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 144


# 선형회귀에서의 annual_inc
from pdpbox.pdp import pdp_isolate, pdp_plot

feature = 'annual_inc'


isolated = pdp_isolate(
    model=linear, 
    dataset=X_val, 
    model_features=X_val.columns, 
    feature=feature,
    grid_type='percentile', # default='percentile', or 'equal'
    num_grid_points=10 # default=10
)
pdp_plot(isolated, feature_name=feature);

# Gradiant Boosting에서의 annual_inc
isolated = pdp_isolate(
    model=boosting, 
    dataset=X_val_encoded, 
    model_features=X_val_encoded.columns, 
    feature=feature
)

pdp_plot(isolated, feature_name=feature);

# 일부분 확대
pdp_plot(isolated, feature_name=feature)
plt.xlim((20000,150000));

10개의 ICE(Individual Conditional Expectation) curves

한 ICE 곡선은 하나의 관측치에 대해 관심 특성을 변화시킴에 따른 타겟값 변화 곡선이고 이 ICE의 평균이 PDP이다.

pdp_plot(isolated
         , feature_name=feature
         , plot_lines=True # ICE plots
         , frac_to_plot=0.001 # or 10 (# 10000 val set * 0.001)
         , plot_pts_dist=True) 

plt.xlim(20000,150000);

X_val_encoded['annual_inc'].value_counts()
'''
60000.00    391
50000.00    388
65000.00    290
70000.00    282
40000.00    271
           ... 
52850.00      1
39998.40      1
29498.00      1
65101.71      1
15360.00      1
Name: annual_inc, Length: 1386, dtype: int64
'''

ICE curves -> PDP를 표현하는 GIF(Christoph Molnar)
한 특성에 대해 PDP를 그릴 경우 데이터셋 사이즈에 grid points를 곱한 수만큼 예측해야한다.

isolated = pdp_isolate(
    model=boosting, 
    dataset=X_val_encoded, 
    model_features=X_val.columns, 
    feature=feature,
    # grid point를 크게 주면 겹치는 점이 생겨 Number of unique grid points는 grid point 보다 작을 수 있다.
    num_grid_points=100, # grid 포인트를 더 줄 수 있다. default = 10
)


isolated = pdp_isolate(
    model=boosting, 
    dataset=X_val_encoded, 
    model_features=X_val.columns, 
    feature=feature,
    # grid point를 크게 주면 겹치는 점이 생겨 Number of unique grid points는 grid point 보다 작을 수 있습니다.
    num_grid_points=100, # grid 포인트를 더 줄 수 있습니다. default = 10
)
'''
예측수:  1000000
'''


pdp_plot(isolated
         , feature_name=feature
         , plot_lines=True
         , frac_to_plot=0.01 # ICE curves는 100개
         , plot_pts_dist=True )

plt.xlim(20000,150000);

두 특성 사용

(참고: PDPBox version <= 0.20 과 몇몇 matplotlib 버전에서 pdp_interact_plot에서plot_type='contour' 설정시 에러가 발생할 수 있다. TypeError: clabel() got an unexpected keyword argument 'contour_label_fontsize' https://github.com/SauceCat/PDPbox/issues/40) —

from pdpbox.pdp import pdp_interact, pdp_interact_plot

features = ['annual_inc', 'fico_range_high']

interaction = pdp_interact(
    model=boosting, 
    dataset=X_val_encoded,
    model_features=X_val.columns, 
    features=features
)

pdp_interact_plot(interaction, plot_type='grid', 
                  feature_names=features);

# Plotly로 3D 구현
features
'''
['annual_inc', 'fico_range_high']
'''


# 2D PDP dataframe
interaction.pdp

type(interaction.pdp)
'''
pandas.core.frame.DataFrame
'''

# 위에서 만든 2D PDP를 테이블로 변환(using Pandas, df.pivot_table)하여 사용

pdp = interaction.pdp.pivot_table(
    values='preds', # interaction['preds']
    columns=features[0], 
    index=features[1]
)[::-1] # 인덱스를 역순으로 만드는 slicing


pdp

# 양단에 극단적인 annual_inc를 drop
pdp = pdp.drop(columns=[1764.0, 1500000.0])


import plotly.graph_objs as go

surface = go.Surface(
    x=pdp.columns, 
    y=pdp.index, 
    z=pdp.values
)


layout = go.Layout(
    scene=dict(
        xaxis=dict(title=features[0]), 
        yaxis=dict(title=features[1]), 
        zaxis=dict(title=target)
    )
)

fig = go.Figure(surface, layout)
fig.show()

PDP Categiry 특성 사용

카테고리 특성을 학습할 때 Ordina Encoder, Target Encoder 등의 인코더를 사용한다
Encoding을 하면 학습 후 PDP를 그릴 때 Encoding된 값이 나오게 되어 카테고리 특성의 실제값을 확인하기 어렵다.

PDP에 Encoding되기 전 Category 값을 보여주기 위한 방법

# Titanic Dataset
import seaborn as sns
from sklearn.ensemble import RandomForestClassifier
from sklearn.pipeline import make_pipeline

df = sns.load_dataset('titanic')
df['age'] = df['age'].fillna(df['age'].median())
df = df.drop(columns='deck') # NaN 77%
df = df.dropna()

target = 'survived'
features = df.columns.drop(['survived', 'alive'])

X = df[features]
y = df[target]


pipe = make_pipeline(
    OrdinalEncoder(), 
    RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)
)
pipe.fit(X, y);


encoder = pipe.named_steps['ordinalencoder']
X_encoded = encoder.fit_transform(X)
rf = pipe.named_steps['randomforestclassifier']


import matplotlib.pyplot as plt
from pdpbox import pdp
feature = 'sex'
pdp_dist = pdp.pdp_isolate(model=rf, dataset=X_encoded, model_features=features, feature=feature)
pdp.pdp_plot(pdp_dist, feature); # 인코딩된 sex 값 확인

# encoder 맵핑을 확인, {male:1, female:2}
encoder.mapping
'''
[{'col': 'sex',
  'mapping': male      1
  female    2
  NaN      -2
  dtype: int64,
  'data_type': dtype('O')},
 {'col': 'embarked',
  'mapping': S      1
  C      2
  Q      3
  NaN   -2
  dtype: int64,
  'data_type': dtype('O')},
 {'col': 'class',
  'mapping': Third     1
  First     2
  Second    3
  NaN      -2
  dtype: int64,
  'data_type': CategoricalDtype(categories=['First', 'Second', 'Third'], ordered=False)},
 {'col': 'who',
  'mapping': man      1
  woman    2
  child    3
  NaN     -2
  dtype: int64,
  'data_type': dtype('O')},
 {'col': 'embark_town',
  'mapping': Southampton    1
  Cherbourg      2
  Queenstown     3
  NaN           -2
  dtype: int64,
  'data_type': dtype('O')}]
'''


pdp.pdp_plot(pdp_dist, feature)

# xticks labels 설정을 인코딩된 코드리스트와, 카테고리 값 리스트를 수동으로
plt.xticks([1, 2], ['male', 'female',]);

# 이번에는 PDP 카테고리값 맵핑 자동으로

feature = 'sex'
for item in encoder.mapping:
    if item['col'] == feature:
        feature_mapping = item['mapping'] # Series
        
feature_mapping = feature_mapping[feature_mapping.index.dropna()]
category_names = feature_mapping.index.tolist()
category_codes = feature_mapping.values.tolist()

pdp.pdp_plot(pdp_dist, feature)

# xticks labels 설정을 위한 리스트를 직접 넣지 않아도 됨
plt.xticks(category_codes, category_names);

# 2D PDP
features = ['sex', 'age']

interaction = pdp_interact(
    model=rf, 
    dataset=X_encoded, 
    model_features=X_encoded.columns, 
    features=features
)

pdp_interact_plot(interaction, plot_type='grid', feature_names=features);

# 2D PDP 를 Seaborn Heatmap으로 그리기 위해 데이터프레임으로 만듭니다
pdp = interaction.pdp.pivot_table(
    values='preds', 
    columns=features[0], 
    index=features[1]
)[::-1]

pdp = pdp.rename(columns=dict(zip(category_codes, category_names)))
plt.figure(figsize=(6,5))
sns.heatmap(pdp, annot=True, fmt='.2f', cmap='viridis')
plt.title('PDP decoded categorical');

SHAP

어떤 ML Model이든 단일 관측치로부터 특성들의 기여도(Feature attribution)를 계산하기 위한 방법이다.
Shapley valu는 원래 게임이론에서 나온 개념이지만 복잡한 ML model의 예측을 설명하기 위한 유용한 방법이다.

# 회귀 모델 예시
import numpy as np
import pandas as pd

# 킹카운티 주택가격
df = pd.read_csv('https://ds-lecture-data.s3.ap-northeast-2.amazonaws.com/kc_house_data/kc_house_data.csv')

# price, longitude, latitude 양 끝단 값 1% 제거
# Remove the most extreme 1% prices,
# the most extreme .1% latitudes, &
# the most extreme .1% longitudes
df = df[(df['price'] >= np.percentile(df['price'], 0.5)) & 
        (df['price'] <= np.percentile(df['price'], 99.5)) & 
        (df['long'] >= np.percentile(df['long'], 0.05)) & 
        (df['long'] <= np.percentile(df['long'], 99.95)) &
        (df['lat'] >= np.percentile(df['lat'], 0.05)) & 
        (df['lat'] < np.percentile(df['lat'], 99.95))]

# split train/test, 2015-03-01 기준으로 나누기
df['date'] = pd.to_datetime(df['date'], infer_datetime_format=True)
cutoff = pd.to_datetime('2015-03-01')
train = df[df['date'] < cutoff]
test  = df[df['date'] >= cutoff]


train.shape, test.shape
'''
((16660, 21), (4691, 21))
'''


train.columns
'''
Index(['id', 'date', 'price', 'bedrooms', 'bathrooms', 'sqft_living',
       'sqft_lot', 'floors', 'waterfront', 'view', 'condition', 'grade',
       'sqft_above', 'sqft_basement', 'yr_built', 'yr_renovated', 'zipcode',
       'lat', 'long', 'sqft_living15', 'sqft_lot15'],
      dtype='object')
'''


features = ['bedrooms', 'bathrooms', 'long', 'lat']
target = 'price'
X_train = train[features]
y_train = train[target]
X_test = test[features]
y_test = test[target]


from scipy.stats import randint, uniform
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import RandomizedSearchCV

param_distributions = { 
    'n_estimators': randint(50, 500), 
    'max_depth': [5, 10, 15, 20, None], 
    'max_features': uniform(0, 1), 
}

search = RandomizedSearchCV(
    RandomForestRegressor(random_state=2), 
    param_distributions=param_distributions, 
    n_iter=5, 
    cv=3, 
    scoring='neg_mean_absolute_error', 
    verbose=10, 
    return_train_score=True, 
    n_jobs=-1, 
    random_state=2
)

search.fit(X_train, y_train);
'''
Fitting 3 folds for each of 5 candidates, totalling 15 fits
'''


print('최적 하이퍼파라미터: ', search.best_params_)
print('CV MAE: ', -search.best_score_)
model = search.best_estimator_
'''
최적 하이퍼파라미터:  {'max_depth': 15, 'max_features': 0.6327377306009369, 'n_estimators': 166}
CV MAE:  101224.42224844794
'''

Shapley values

게임이론에서 같은 팀 선수들(특성들)이 게임 목표(예측) 달성을 위해 각자 자신의 역할(기여)을 한다고 할 때 게임 목표 달성 후 받은 포상을 어떻게 하면 그들의 기여도에 따라 공평하게 나누어 줄 수 있을 것인가? 라는 질문과 연관된다.
Shapley value를 ML의 특성 기여도산정에 활용

특성 갯수가 많아질 수록 Shapley value를 구할 때 필요한 계산량이 기하급수적으로 늘어난다.
SHAP에서는 샘플링을 이용해 근사적으로 값을 구합니다

# Test set 두번째 sample의 Shap value
row = X_test.iloc[[1]]  # 중첩 brackets을 사용하면 결과물이 DataFrame
row

# 실제 집값
y_test.iloc[[1]] # 2번째 데이터를 사용했습니다
'''
9    323000.0
Name: price, dtype: float64
'''


# 모델 예측값
model.predict(row)
'''
array([341878.50142523])
'''


# 모델이 이렇게 예측한 이유를 알기 위하여 SHAP Force Plot을 그린다.
import shap

explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(row)

shap.initjs()
shap.force_plot(
    base_value=explainer.expected_value, 
    shap_values=shap_values,
    features=row
)

# 집 가격 평균값 base value
explainer.expected_value[0]
'''
525264.9249674568
'''


# 이 관측치의 예측값이 왜 341,878.50 이 나오게 되었는지 각 특성(bathrooms, lat, bedrooms)의 영향을 시각화

# 예측함수 정의
def predict(bedrooms, bathrooms, longitude, latitude):

    # 함수 내에서 예측에 사용될 input 생성
    df = pd.DataFrame(
        data=[[bedrooms, bathrooms, longitude, latitude]], 
        columns=['bedrooms', 'bathrooms', 'long', 'lat']
    )

    # 예측
    pred = model.predict(df)[0]

    # Shap value 계산
    explainer = shap.TreeExplainer(model)
    shap_values = explainer.shap_values(df)

    # Shap value, 특성이름, 특성값을 가지는 Series 생성
    feature_names = df.columns
    feature_values = df.values[0]
    shaps = pd.Series(shap_values[0], zip(feature_names, feature_values))

    # 결과 프린트
    result = f'평균가격: ${explainer.expected_value[0]:,.0f} \n'
    result += f'예측가격: ${pred:,.0f}. \n'
    result += shaps.to_string()
    print(result)


    # SHAP Force Plot
    shap.initjs()
    return shap.force_plot(
        base_value=explainer.expected_value, 
        shap_values=shap_values, 
        features=df
    )


# lat 분포
df['lat'].plot.hist(bins=100, figsize=(4, 2));

# 적당한 지역의 방 3개인 집 가격
predict(3, 1, -121.35, 47.55)
'''
평균가격: $525,265 
예측가격: $382,123. 
(bedrooms, 3.0)     -25124.723574
(bathrooms, 1.0)   -142083.969321
(long, -121.35)     -21022.116137
(lat, 47.55)         45088.940133
'''

# 같은 지역에 방 2개 집 가격을 예측해 보면, 지역(lat) 수치가 같음에도 영향은 달라짐
predict(2, 1, -122.35, 47.55)
'''
평균가격: $525,265 
예측가격: $281,714. 
(bedrooms, 2.0)     -45592.182391
(bathrooms, 1.0)   -118088.150603
(long, -122.35)     -62003.221763
(lat, 47.55)        -17867.074728
'''

# 같은 지역 방 1개인 집 가격
predict(1, 1, -122.35, 47.55)
'''
평균가격: $525,265 
예측가격: $277,940. 
(bedrooms, 1.0)     -54303.671169
(bathrooms, 1.0)   -120596.652648
(long, -122.35)     -60103.932089
(lat, 47.55)        -12320.268028
'''

# SHAP plot으로 각 특성이 어떤 값 범위에서 영향을 주는지 확인
# 100개 테스트 샘플에 대해서 각 특성들의 영향 확인
# 샘플 수를 너무 크게 잢으면 계산이 오래걸리니 주의
shap_values = explainer.shap_values(X_test.iloc[:100])
shap.force_plot(explainer.expected_value, shap_values, X_test.iloc[:100])

shap_values = explainer.shap_values(X_test.iloc[:300])
shap.summary_plot(shap_values, X_test.iloc[:300])

shap.summary_plot(shap_values, X_test.iloc[:300], plot_type="violin")

shap.summary_plot(shap_values, X_test.iloc[:300], plot_type="bar")

SHAP value 분류 문제에 적용

# Lending Club 데이터 사용
# 대출 상태가 'charged off'(상각) 인지 'fully paid'(완납) 인지 예측하는 문제
import pandas as pd

# Kaggle 데이터셋에서 10% 샘플링된 데이터입니다.
## Source: https://www.kaggle.com/wordsforthewise/lending-club
## 10% of expired loans (loan_status: ['Fully Paid' and 'Charged Off'])
## grades A-D
## term ' 36 months'
df = pd.read_csv('https://ds-lecture-data.s3.ap-northeast-2.amazonaws.com/lending_club/lending_club_sampled.csv', index_col=0)


# 2-class 타겟 ('Fully Paid' or 'Charged Off')
target = 'loan_status'
X = df.drop(columns=target)
y = df[target]


# 데이터셋 분리
from sklearn.model_selection import train_test_split

X_train_val, X_test, y_train_val, y_test = train_test_split(
    X, y, test_size=10000
    , stratify=y
    , random_state=2)

X_train, X_val, y_train, y_val = train_test_split(
    X_train_val, y_train_val, test_size=10000
    , stratify=y_train_val
    , random_state=42)

print('X_train shape', X_train.shape)
print('y_train shape', y_train.shape)
print('X_val shape', X_val.shape)
print('y_val shape', y_val.shape)
print('X_test shape', X_test.shape)
print('y_test shape', y_test.shape)
'''
X_train shape (77591, 150)
y_train shape (77591,)
X_val shape (10000, 150)
y_val shape (10000,)
X_test shape (10000, 150)
y_test shape (10000,)
'''


# test ids를 저장하고 SHAP분석시 사용
test_id = X_test['id']


def wrangle(X):
    X = X.copy()

    # to datetime
    X['issue_d'] = pd.to_datetime(X['issue_d'], infer_datetime_format=True)
    
    # 개설 날짜 - 최초 신용 개설
    X['earliest_cr_line'] = pd.to_datetime(X['earliest_cr_line'], infer_datetime_format=True)
    X['earliest_cr_line'] = X['issue_d'] - X['earliest_cr_line']
    X['earliest_cr_line'] = X['earliest_cr_line'].dt.days

    # Engineer issue_d_year
    X['issue_d_year'] = X['issue_d'].dt.year
    
    # Engineer issue_d_year
    X['issue_d_month'] = X['issue_d'].dt.month
            
    # non-digit 문자 치환 -> '', float변환
    X['emp_length'] = X['emp_length'].str.replace(r'\D','').astype(float)
        
    # Get length of free text fields
    X['title'] = X['title'].str.len()
    X['desc'] = X['desc'].str.len()
    X['emp_title'] = X['emp_title'].str.len()
    
    # sub_grade 숫자 치환
    sub_grade_ranks = {'A1': 1.1, 'A2': 1.2, 'A3': 1.3, 'A4': 1.4, 'A5': 1.5, 
                       'B1': 2.1, 'B2': 2.2, 'B3': 2.3, 'B4': 2.4, 'B5': 2.5, 
                       'C1': 3.1, 'C2': 3.2, 'C3': 3.3, 'C4': 3.4, 'C5': 3.5, 
                       'D1': 4.1, 'D2': 4.2, 'D3': 4.3, 'D4': 4.4, 'D5': 4.5}
    X['sub_grade'] = X['sub_grade'].map(sub_grade_ranks)
    
    # 크게 의미 없는 특성 삭제
    X = X.drop(columns='id')        # Always unique
    X = X.drop(columns='url')       # Always unique
    X = X.drop(columns='grade')     # Duplicative of sub_grade
    X = X.drop(columns='zip_code')  # High cardinality
    X = X.drop(columns='issue_d')   # date
    
    # drop null > 70%
    null_frac = X.isnull().mean().sort_values(ascending=False)
    X = X.drop(columns = sorted(list(null_frac[null_frac > 0.7].index)))
    
    # Keep list (https://www.kaggle.com/pileatedperch/predicting-charge-off-from-initial-listing-data)
    # 잠재적인 투자자에게만 제공되는 특성으로만 제한합니다
    keep_list = keep_list = [
        'addr_state', 'annual_inc', 'application_type', 'dti', 'earliest_cr_line', 'emp_length'
        , 'emp_title', 'fico_range_high', 'fico_range_low', 'grade', 'home_ownership', 'id'
        , 'initial_list_status', 'installment', 'int_rate', 'issue_d', 'loan_amnt', 'loan_status'
        , 'mort_acc', 'open_acc', 'pub_rec', 'pub_rec_bankruptcies', 'purpose', 'revol_bal'
        , 'revol_util', 'sub_grade', 'term', 'title', 'total_acc', 'verification_status', 'zip_code']
    drop_list = [col for col in X.columns if col not in keep_list]
    X = X.drop(labels=drop_list, axis=1)
        
    # Reset index
    X = X.reset_index(drop=True)
    
    return X

X_train = wrangle(X_train)
X_val   = wrangle(X_val)
X_test  = wrangle(X_test)

print('X_train shape', X_train.shape)
print('X_val shape', X_val.shape)
print('X_test shape', X_test.shape)
'''
X_train shape (77591, 26)
X_val shape (10000, 26)
X_test shape (10000, 26)
'''


# 클래스의 비율
y_train.value_counts(normalize=True)
'''
Fully Paid     0.847328
Charged Off    0.152672
Name: loan_status, dtype: float64
'''


ratio = 0.15/0.84
ratio
'''
0.17857142857142858
'''


from category_encoders import OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from xgboost import XGBClassifier

processor = make_pipeline(
    OrdinalEncoder(), 
    SimpleImputer(strategy='median')
)

X_train_processed = processor.fit_transform(X_train)
X_val_processed = processor.transform(X_val)

eval_set = [(X_train_processed, y_train), 
            (X_val_processed, y_val)]

# XGBoost 분류기를 학습
# 클래스 비율을 맞추기 위해 scale_pos_weight= #C harged Off / # Fully Paid
model = XGBClassifier(n_estimators=1000, verbosity=0, n_jobs=-1, scale_pos_weight=ratio)
model.fit(X_train_processed, y_train, eval_set=eval_set, eval_metric='auc', 
          early_stopping_rounds=10)
'''
[0]	validation_0-auc:0.68092	validation_1-auc:0.66370
Multiple eval metrics have been passed: 'validation_1-auc' will be used for early stopping.

Will train until validation_1-auc hasn't improved in 10 rounds.
[1]	validation_0-auc:0.68988	validation_1-auc:0.67363
[2]	validation_0-auc:0.69617	validation_1-auc:0.67690
[3]	validation_0-auc:0.70074	validation_1-auc:0.67876
.
.
.
[35]	validation_0-auc:0.77165	validation_1-auc:0.68821
[36]	validation_0-auc:0.77249	validation_1-auc:0.68820
[37]	validation_0-auc:0.77412	validation_1-auc:0.68758
Stopping. Best iteration:
[27]	validation_0-auc:0.75906	validation_1-auc:0.68906

XGBClassifier(base_score=0.5, booster=None, colsample_bylevel=1,
              colsample_bynode=1, colsample_bytree=1, gamma=0, gpu_id=-1,
              importance_type='gain', interaction_constraints=None,
              learning_rate=0.300000012, max_delta_step=0, max_depth=6,
              min_child_weight=1, missing=nan, monotone_constraints=None,
              n_estimators=1000, n_jobs=-1, num_parallel_tree=1, random_state=0,
              reg_alpha=0, reg_lambda=1, scale_pos_weight=0.17857142857142858,
              subsample=1, tree_method=None, validate_parameters=False,
              verbosity=0)
'''


from sklearn.metrics import roc_auc_score
X_test_processed = processor.transform(X_test)
X_val_processed = processor.transform(X_val)
class_index = 1
y_pred_proba = model.predict_proba(X_test_processed)[:, class_index]
print(f'Test AUC for class "{model.classes_[class_index]}":')
print(roc_auc_score(y_test, y_pred_proba)) # 범위는 0-1, 수치는 높을 수록 좋습니다
'''
Test AUC for class "Fully Paid":
0.6839267781606986
'''


# Confution matrix를 확인해 봅시다
from sklearn.metrics import classification_report
y_test_pred = model.predict(X_test_processed)
print(classification_report(y_test, y_test_pred))
'''
              precision    recall  f1-score   support

 Charged Off       0.24      0.66      0.35      1527
  Fully Paid       0.91      0.61      0.73      8473

    accuracy                           0.62     10000
   macro avg       0.57      0.64      0.54     10000
weighted avg       0.81      0.62      0.67     10000
'''


# 예측값 실제값 비교
df_p = pd.DataFrame({
    'id': test_id, 
    'pred_proba': y_pred_proba, # 예측확률 
    'status_group': y_test # 실제값
})

df_p = df_p.merge(
     df[['id','issue_d','sub_grade','total_pymnt','funded_amnt', 'term','int_rate']],
     how='left'
)


df_p.head()

fully_paid = df_p['status_group'] == 'Fully Paid'
charged_off = ~fully_paid
right = (fully_paid) == (df_p['pred_proba'] > 0.50)
wrong = ~right


# 대출은 Fully Paid, 예측이 맞는 경우
df_p[fully_paid & right].sample(n=10, random_state=1).sort_values(by='pred_proba')

# 테스트셋에서 인덱스 1 샘플의 예측
# 우선 모든 특성 수치를 본다
row = X_test.iloc[[3160]]
row

# SHAP 그래프로 예측 설명
## UnicodeDecoderError 발생시 xgboost 1.1-> 1.0 다운그레이드 (conda install -c conda-forge xgboost=1.0)
import xgboost
import shap

explainer = shap.TreeExplainer(model)
row_processed = processor.transform(row)
shap_values = explainer.shap_values(row_processed)

shap.initjs()
shap.force_plot(
    base_value=explainer.expected_value, 
    shap_values=shap_values, 
    features=row, 
    link='logit' # SHAP value를 확률로 변환
)

# 예측을 SHAP 그래프를 통해 설명하는 함수
feature_names = row.columns
feature_values = row.values[0]
shaps = pd.Series(shap_values[0], zip(feature_names, feature_values))


pros = shaps.sort_values(ascending=False)[:3].index
cons = shaps.sort_values(ascending=True)[:3].index


print('fully paid 예측에 대한 Positive 요인 Top 3 입니다:')
for i, pro in enumerate(pros, start=1):
    feature_name, feature_value = pro
    print(f'{i}. {feature_name} : {feature_value}')

print('\n')
print('Negative 요인 Top 3 입니다:')
for i, con in enumerate(cons, start=1):
    feature_name, feature_value = con
    print(f'{i}. {feature_name} : {feature_value}')
'''
fully paid 예측에 대한 Positive 요인 Top 3 입니다:
1. sub_grade : 1.1
2. int_rate : 5.42
3. fico_range_low : 745.0


Negative 요인 Top 3 입니다:
1. open_acc : 16.0
2. purpose : major_purchase
3. addr_state : PA
'''


def explain(row_number):
    positive_class = 'Fully Paid'
    positive_class_index = 1

    # row 값을 변환합니다
    row = X_test.iloc[[row_number]]
    row_processed = processor.transform(row)

    # 예측하고 예측확률을 얻습니다 
    pred = model.predict(row_processed)[0]
    pred_proba = model.predict_proba(row_processed)[0, positive_class_index]
    pred_proba *= 100
    if pred != positive_class:
        pred_proba = 100 - pred_proba

    # 예측결과와 확률값을 얻습니다
    print(f'이 대출에 대한 예측결과는 {pred} 으로, 확률은 {pred_proba:.0f}% 입니다.')
    
    # SHAP를 추가합니다
    shap_values = explainer.shap_values(row_processed)

    # Fully Paid에 대한 top 3 pros, cons를 얻습니다
    feature_names = row.columns
    feature_values = row.values[0]
    shaps = pd.Series(shap_values[0], zip(feature_names, feature_values))
    pros = shaps.sort_values(ascending=False)[:3].index
    cons = shaps.sort_values(ascending=True)[:3].index

    # 예측에 가장 영향을 준 top3
    print('\n')
    print('Positive 영향을 가장 많이 주는 3가지 요인 입니다:')
    
    evidence = pros if pred == positive_class else cons
    for i, info in enumerate(evidence, start=1):
        feature_name, feature_value = info
        print(f'{i}. {feature_name} : {feature_value}')

    # 예측에 가장 반대적인 영향을 준 요인 top1
    print('\n')
    print('Negative 영향을 가장 많이 주는 3가지 요인 입니다:')
    
    evidence = cons if pred == positive_class else pros
    for i, info in enumerate(evidence, start=1):
        feature_name, feature_value = info
        print(f'{i}. {feature_name} : {feature_value}')

    # SHAP
    shap.initjs()
    return shap.force_plot(
        base_value=explainer.expected_value, 
        shap_values=shap_values, 
        features=row, 
        link='logit'
    )


explain(3160)
'''
이 대출에 대한 예측결과는 Fully Paid 으로, 확률은 90% 입니다.


Positive 영향을 가장 많이 주는 3가지 요인 입니다:
1. sub_grade : 1.1
2. int_rate : 5.42
3. fico_range_low : 745.0


Negative 영향을 가장 많이 주는 3가지 요인 입니다:
1. open_acc : 16.0
2. purpose : major_purchase
3. addr_state : PA
'''

# 대출 결과는 Fully Paid, 예측이 잘못된 경우
df_p[fully_paid & wrong].sample(n=10, random_state=1).sort_values(by='pred_proba')

explain(3559)
'''
이 대출에 대한 예측결과는 Charged Off 으로, 확률은 83% 입니다.


Positive 영향을 가장 많이 주는 3가지 요인 입니다:
1. sub_grade : 4.4
2. installment : 1300.55
3. int_rate : 19.99


Negative 영향을 가장 많이 주는 3가지 요인 입니다:
1. home_ownership : MORTGAGE
2. emp_length : 10.0
3. revol_bal : 21434.0
'''

# 대출 결과는 Charged Off, 예측이 잘못된 경우
df_p[charged_off & wrong].sample(n=10, random_state=1).sort_values(by='pred_proba')

explain(7896)
'''
이 대출에 대한 예측결과는 Fully Paid 으로, 확률은 91% 입니다.


Positive 영향을 가장 많이 주는 3가지 요인 입니다:
1. sub_grade : 1.3
2. int_rate : 6.89
3. dti : 5.72


Negative 영향을 가장 많이 주는 3가지 요인 입니다:
1. annual_inc : 31000.0
2. emp_length : 3.0
3. earliest_cr_line : 8705
'''

Feature Importances, PDP, SHAP의 특징 구분

서로 관련이 있는 모든 특성들에 대한 전역적인(Global) 설명
- Feature Importances
- Drop-Column Importances
- Permutaton Importances
타겟과 관련이 있는 개별 특성들에 대한 전역적인 설명
- Partial Dependence plots
- 개별 관측치에 대한 지역적인(local) 설명
Shapley Values