[sklearn][SHAP][モデル可視化]マーケティングデータを分析してみた

概要

データ分析の一連を解説してみます。

はじめに

一通り、データセットがある状態から、
前処理、モデル作成、結果の解釈まで一通りよく行われるデータ分析の流れについて解説していきたいと思います。

使用するデータセット

https://www.kaggle.com/prabhakar01/telemarketing-case
これはテレマーケティングのケーススタディです。

ケースの目的は、
利用可能なフィールドのすべての可能な組み合わせを使用して販売パターンを特定し、全体的な売上（収益）を促進する学習を実装することです。
機会を利用して、予測モデリングのためのビジネス開発ケースを準備します。ビジネスに実装できるさまざまな利用可能な技術を探索します。モデルの目的、追加のデータ要件のリスト、モデルの説明

初期設定

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt

CSV_PATH = './data/dataset_X.csv'

データ読み込み

# 先頭5行
df = pd.read_csv(CSV_PATH)
df.head()

データの確認

データセットの概要、各カラムの中身を確認する。

• 概要
  • トータルレコード
    • 10,000件
  • 正解ラベルの割合（契約あり/なし）
    • 91,248 / 8,720 件
  • 欠損率
    • 約1%未満
  • 各カラム
  • 担当者
    • 244 (名)
  • 製品タイプ数
    • 10 (種類)
  • 地域数
    • 3310 (エリア)
  • 年齢
    • 25～80 (歳) 平均52歳
  • コール回数
    • 1~55 (回) 平均3回
  • タイムゾーン
    • 2 (種)
  • Phoneコード
    • 1 (種)

# 各カラムのSummary
summary = pd.DataFrame(df.fillna('nan').apply(lambda x:set(x)), columns=['value'])
summary['unique'] = summary['value'].apply(lambda x: len(x))
summary['null_cnt'] = df.isnull().sum().T
summary['null_rate'] = summary['null_cnt'].apply(lambda x: x/len(df))
summary

# 各カラムの統計情報
df.describe()

sns.countplot(data=df, x='Sale')
for value in df['Sale'].unique():
    value_len = len(df[df['Sale'] == value])
    print(f'{value}: {value_len}')

前処理

モデルに入れるための前処理する。

• 不必要なカラムの削除
  • 意味のないカラムを削除する
• 欠損があるレコードの削除
  • レコード数が十分にあり欠損率も小さい為、欠損は補完せずに削除する
• 値のラベル化
  • Stringを数値に変換する
• ダミー変数化
  • カテゴリごとにカラムにする

不要なカラム削除

df_processed = df
df_processed = df_processed.drop('Call_ID', axis=1)
df_processed = df_processed.drop('Phone_code', axis=1)
df_processed = df_processed.drop('Timezone', axis=1)
df_processed.head()

欠損レコードの削除

契約の有無とエージェントIDが欠損しているレコードを削除

df_processed = df_processed.dropna(subset=['Sale', 'Agent_ID'])
df_processed = df_processed[df_processed['Gender']!='Others']
print(f'before:{len(df)} after:{len(df_processed)}')

値のラベル化

Stringの値を数値に変換

df_processed['Sale'] = df_processed['Sale'].map({True:1, False:0})
df_processed['Gender'] = df_processed['Gender'].map({'Male':1, 'Female':0})

change_columns = ['Agent_ID','Product_ID', 'First_Name','Last_Name','Area_Code']
for column in change_columns:
    labels, uniques = pd.factorize(df_processed[column])
    df_processed[column] = labels

df_processed.head()

変数相関の確認

変数同士の相関を確認する

• Saleと強い相関があるパラメータはない
  • パラメータ直接答えになっているパラメータがないことが確認できた
• FirstNameには、LastNameと性別に相関がみられる
  • 名前と性別に相関があることは推測できるが、苗字との相関についてはなぜ？
• 電話回数には、製品IDとエージェントIDとの相関がみられる。
  • 製品や人によって電話回数が異なる

sns.heatmap(df_processed.corr(), annot=True, cmap='coo[l')

ダミー変数化

複数カテゴリのカラムをダミー変数化する。
※ 今回は使用せず。精度は上がると思われるが、パラメータが多すぎて解釈が煩雑になる。

df_dummy = pd.get_dummies(df_processed, columns=['Agent_ID', 'Product_ID', 'First_Name', 'Last_Name', 'Area_Code'])
df_dummy.head()

ダウンサンプリング

契約の有無の数に偏りがある為、少ない方にダウンサンプリングする

df_sales = df_processed.sample(frac=1) # シャッフル
mim_count = df_sales['Sale'].value_counts().iat[-1] # 一番少ないクラス数
df_sales = df_sales.groupby('Sale')
df_sales = df_sales.head(mim_count)
sns.countplot(data=df_sales, x='Sale')

データセット作成

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# 標準化
sc = StandardScaler() 
x = sc.fit_transform(np.array(df_sales.iloc[:, 1:]))
y = np.array(df_sales.iloc[:, 0])

# データセットの分割
x_train, x_test, y_train, y_test = train_test_split(x, y, train_size=0.8, random_state=0, stratify=y)

モデル作成

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()
lr.fit(x_train, y_train)

評価

モデルの判別結果を各指標を用いて評価する。
■ 結果
正解率: 約56％
AUC: 0.58

■ 考察
説明変数と目的変数との相関がほとんどなく、単純なモデルだったのでこれくらいの精度で妥当だと考察される。
Recallの割合が高かったので、契約獲得時の特徴はとらえられている可能性がある。
精度を上げる手段としては、カテゴリのダミー変数化、上位モデルの使用（xgBoostやLightGBMなどのブースティング系やDeepLearning）が挙げられる。

混同行列

from sklearn.metrics import confusion_matrix, accuracy_score, precision_score, recall_score, f1_score

y_pred = lr.predict(x_test)
print('confusion matrix = \n', confusion_matrix(y_true=y_test, y_pred=y_pred))
print('accuracy = ', accuracy_score(y_true=y_test, y_pred=y_pred))
print('precision = ', precision_score(y_true=y_test, y_pred=y_pred))
print('recall = ', recall_score(y_true=y_test, y_pred=y_pred))
print('f1 score = ', f1_score(y_true=y_test, y_pred=y_pred))

AUC_ROC曲線

from sklearn.metrics import roc_curve, auc

y_score = lr.predict_proba(x_test)[:, 1] # 検証データがクラス1に属する確率
fpr, tpr, thresholds = roc_curve(y_true=y_test, y_score=y_score)

plt.plot(fpr, tpr, label='roc curve (area = %0.3f)' % auc(fpr, tpr))
plt.plot([0, 1], [0, 1], linestyle='--', label='random')
plt.plot([0, 0, 1], [0, 1, 1], linestyle='--', label='ideal')
plt.legend()
plt.xlabel('false positive rate')
plt.ylabel('true positive rate')
plt.show()

スコア分布

import plotly.graph_objs as go
from plotly.offline import * 
init_notebook_mode()

def hist_plot(pred_score, y_true, title='likely hood graph', thresh_hold=0.5):
    df_score = pd.DataFrame({'pred_score': pred_score,'y_true':y_true})
    true_score = df_score[df_score['y_true']==1]['pred_score']
    false_score = df_score[df_score['y_true']==0]['pred_score']

    # plot
    fig = go.Figure()
    fig.add_trace(go.Histogram(name='False', x=false_score, xbins=dict(start=0,end=1,size=0.01), histnorm='percent'))
    fig.add_trace(go.Histogram(name='True', x=true_score, xbins=dict(start=0,end=1,size=0.01), histnorm='percent'))
    fig.update_layout(width=800, height=600, xaxis=dict(range=[0,1]), title=f'{title} (thresh_hold={thresh_hold})', barmode='overlay')
    fig.update_traces(opacity=0.5)
    fig.add_shape(type='line', x0=thresh_hold, x1=thresh_hold, y0=0, y1=10)
    fig.show()

hist_plot(y_score, y_test)

解釈

SHAPを使用して、どのパラメータが判定結果に影響したのかを可視化す

SHAP

import shap
shap.initjs()

explainer = shap.LinearExplainer(lr, x_train, feature_dependence="independent")
shap_values = explainer.shap_values(x_test)

モデル全体

モデル全体で判定結果の解釈を確認する。

• 電話回数が多いと契約がとれにくい
• 商品が契約の有無に影響している
• 女性より男性の方が契約をとりやすい
• 年齢が若い方が契約を取りやすい

shap.summary_plot(shap_values, x_test, feature_names=df_sales.columns[1:])    

個別の結果

個別のレコードで解釈結果を確認する。

• 1番スコアが高かったレコード（契約がとれると判断）

  • 担当者、製品、電話回数の順で結果に大きく影響

• 1番スコアが低かったレコード（契約がとれないと判断）

  • 電話回数、名前、製品の順で結果に大きく影響

# モデルのスコアが一番高かったレコードの結果
ind = np.argmax(y_score)
shap.force_plot(
    explainer.expected_value, shap_values[ind,:], x_test[ind,:],
    feature_names=df_processed.columns[1:]
)
# モデルのスコアが一番低かったレコードの結果
ind = np.argmin(y_score)
shap.force_plot(
    explainer.expected_value, shap_values[ind,:], x_test[ind,:],
    feature_names=df_processed.columns[1:]
)

おわりに

だいたいどの分析でも似たような流れになります。