scikit-learnでロジスティック回帰分析を試す

とにかく試して見るシリーズ第二弾。

ロジスティック回帰分析とは

概要

ロジスティック回帰分析 (Logistic regression) は、発生確率を予測する手法。
2値しかとりえない値を目的変数とし、説明変数を用いてその発生確率を説明する。

ロジスティック回帰分析でできること

• 予測値を算出する
• 説明変数の目的変数に対する貢献度を算出する

参考にしたサイト*1

ビジネスでの活用例

• 顧客がある特定のキャンペーンに対して反応を示すか否か、を他の変数を利用して予測する
  → マーケティングの潜在顧客推定に使用されている
• 金融機関では、リスク予測に使用されている
• 医療分野では、複数の因子からある症状を発症する可能性を予測するのに使われていた
  → これがマーケティングに応用されるようになったようだ

取り組んだ課題

ロジスティック回帰に適していそうなデータセットを自力で探し、予測モデルを作ってみる

試行過程と結果

使えそうなデータセットを探す

目的変数が0または1の2値になるデータセットを探した結果、過去に利用したサイトにちょうどいいデータセットがあった。
台湾で実施された、６ヶ月間のクレジットカードの支払履歴と翌月の支払状況の調査結果のようだ。
目的変数となる支払状況は0,1の２値で示せるため、ロジスティック回帰に適していそうだ。
（のちにKaggleに出ているのを発見した）

データセットを確認する

Datasetはここから入手した

データの入手元にあった各項目の説明

default payment next month：翌月滞納 (Yes=1,No=0)  
  
   LIMIT_BAL  ：利用可能枠 (台湾新ドル、家族カードの利用可能枠を含む)  
   SEX        ：性別 (1=male,2=female)  
   EDUCATION  ：最終学歴 (1=graduate school,2=university,3=high school,4=others)  
   MARRIAGE   ：既婚/未婚/その他 (1=married,2=single,3=others)  
   AGE        ：年齢 (year)  
   PAY_0-9    ：支払歴 (-1=遅延なく支払,1=1ヶ月延滞,...,8=8ヶ月延滞,9=9ヶ月以上延滞)  
   BILL_AMT1-6：請求額 (BILL_AMT1=請求額(2005年9月),BILL_AMT2=請求額(2005年8月),...BILL_AMT6=請求額(2005年4月))  
   PAY_AMT1-6 ：支払額 (PAY_AMT1=支払額(2005年9月),PAY_AMT2=支払額(2005年8月),...PAY_AMT6=支払額(2005年4月))  

# .xlsファイルを読込んで、rawを見てみる
import pandas as pd
dframe_in = pd.read_excel('default of credit card clients.xls', sheetname='Data')
dframe_in.head(5)

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10	...	X15	X16	X17	X18	X19	X20	X21	X22	X23	Y
ID	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_0	PAY_2	PAY_3	PAY_4	PAY_5	...	BILL_AMT4	BILL_AMT5	BILL_AMT6	PAY_AMT1	PAY_AMT2	PAY_AMT3	PAY_AMT4	PAY_AMT5	PAY_AMT6	default payment next month
1	20000	2	2	1	24	2	2	-1	-1	-2	...	0	0	0	0	689	0	0	0	0	1
2	120000	2	2	2	26	-1	2	0	0	0	...	3272	3455	3261	0	1000	1000	1000	0	2000	1
3	90000	2	2	2	34	0	0	0	0	0	...	14331	14948	15549	1518	1500	1000	1000	1000	5000	0
4	50000	2	2	1	37	0	0	0	0	0	...	28314	28959	29547	2000	2019	1200	1100	1069	1000	0

5 rows × 24 columns

# 確認できなかった列を参照
from pandas import DataFrame
dframe_in[['X11','X12','X13','X14']].head(5)

	X11	X12	X13	X14
ID	PAY_6	BILL_AMT1	BILL_AMT2	BILL_AMT3
1	-2	3913	3102	689
2	2	2682	1725	2682
3	0	29239	14027	13559
4	0	46990	48233	49291

データの説明について、誤っていることがわかったこと

PAY_0-9[支払歴]は、PAY_0,(飛んで),PAY_2,PAY_3,PAY_4,PAY_5,PAY_6しか存在しない。
欠損値も確認しておく

# 欠損値を確認する
dframe_in.isnull().any().any()

False

欠損値はない

# ヘッダー行、ID列を読み飛ばして再度読み込み
# 紛らわしいため、'PAY_0'は'PAY_1'に置換、'default payment next month'も長いため省略
dframe_e1 = pd.read_excel('default of credit card clients.xls', header=1, sheetname='Data')
dframe_e1 = dframe_e1.drop("ID",axis=1)
dframe_e1.rename(columns={'PAY_0':'PAY_Sep',
                          'PAY_2':'PAY_Aug',
                          'PAY_3':'PAY_Jul',
                          'PAY_4':'PAY_Jun',
                          'PAY_5':'PAY_May',
                          'PAY_6':'PAY_Apr',
                          'BILL_AMT1':'BILL_AMT_Sep',
                          'BILL_AMT2':'BILL_AMT_Aug',
                          'BILL_AMT3':'BILL_AMT_Jul',
                          'BILL_AMT4':'BILL_AMT_Jun',
                          'BILL_AMT5':'BILL_AMT_May',
                          'BILL_AMT6':'BILL_AMT_Apr',
                          'PAY_AMT1':'PAY_AMT_Sep',
                          'PAY_AMT2':'PAY_AMT_Aug',
                          'PAY_AMT3':'PAY_AMT_Jul',
                          'PAY_AMT4':'PAY_AMT_Jun',
                          'PAY_AMT5':'PAY_AMT_May',
                          'PAY_AMT6':'PAY_AMT_Apr',
                          'default payment next month':'DEF_PAY'}, inplace=True)
dframe_e1.head(5)

	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_Sep	PAY_Aug	PAY_Jul	PAY_Jun	PAY_May	...	BILL_AMT_Jun	BILL_AMT_May	BILL_AMT_Apr	PAY_AMT_Sep	PAY_AMT_Aug	PAY_AMT_Jul	PAY_AMT_Jun	PAY_AMT_May	PAY_AMT_Apr	DEF_PAY
0	20000	2	2	1	24	2	2	-1	-1	-2	...	0	0	0	0	689	0	0	0	0	1
1	120000	2	2	2	26	-1	2	0	0	0	...	3272	3455	3261	0	1000	1000	1000	0	2000	1
2	90000	2	2	2	34	0	0	0	0	0	...	14331	14948	15549	1518	1500	1000	1000	1000	5000	0
3	50000	2	2	1	37	0	0	0	0	0	...	28314	28959	29547	2000	2019	1200	1100	1069	1000	0
4	50000	1	2	1	57	-1	0	-1	0	0	...	20940	19146	19131	2000	36681	10000	9000	689	679	0

5 rows × 24 columns

PAY_Apr〜PAY_Sepに'-2'と'0'の説明にない値があるため、規則性を探って見る

# 請求額、支払歴、支払額の順に並べ替える（前月請求額が翌月支払われている模様）
dframe_e2 = dframe_e1[['PAY_AMT_Apr',
                       'PAY_Apr',
                       'BILL_AMT_Apr',
                       'PAY_AMT_May',
                       'PAY_May',
                       'BILL_AMT_May',
                       'PAY_AMT_Jun',
                       'PAY_Jun',
                       'BILL_AMT_Jun',
                       'PAY_AMT_Jul',
                       'PAY_Jul',
                       'BILL_AMT_Jul',
                       'PAY_AMT_Aug',
                       'PAY_Aug',
                       'BILL_AMT_Aug',
                       'PAY_AMT_Sep',
                       'PAY_Sep',
                       'BILL_AMT_Sep'
                     ]]

# -1から-2に転じているデータを見てみる
dframe_e3 = dframe_e2[dframe_e2['PAY_Jun']>-2]
dframe_e3 = dframe_e3[dframe_e3['PAY_Jul']<-1]
dframe_e3.head(5)

	PAY_AMT_Apr	PAY_Apr	BILL_AMT_Apr	PAY_AMT_May	PAY_May	BILL_AMT_May	PAY_AMT_Jun	PAY_Jun	BILL_AMT_Jun	PAY_AMT_Jul	PAY_Jul	BILL_AMT_Jul	PAY_AMT_Aug	PAY_Aug	BILL_AMT_Aug	PAY_AMT_Sep	PAY_Sep	BILL_AMT_Sep
65	1000	2	7918	0	2	8198	300	-1	8174	8222	-2	144076	0	-2	148751	0	-2	152519
68	13899	-1	7319	7319	2	10161	0	-1	10311	20161	-2	-9850	0	-2	-9850	0	1	-190
198	15816	-1	1151	1151	-1	1206	1206	-1	1251	1251	-2	2299	2299	-2	138	138	-2	412
232	2000	0	46557	1747	-1	45567	45567	-1	2624	2624	-2	0	0	0	0	0	0	102800
265	316	-1	316	316	-1	316	316	-1	316	316	-2	0	0	-1	0	0	-1	6156

請求額と支払額が一致しないレコードがあるが、分割かリボルビング払いだろうか？

わからないので検索して見る

そもそも推測なんて意味ないので、調べて見るべきだろう。
Kaggleでデータセットを作成した教授にe-mailして聞いて見たという投稿を発見
https://www.kaggle.com/uciml/default-of-credit-card-clients-dataset/discussion/34608

-2: No consumption(請求なし);
0: The use of revolving credit(リボルビング払い);

やはりリボ払いか…
請求額があって、支払歴が'-2(請求額なし)‘のレコードもあるがこっちは分割払いだろうか…
いずれにせよリボルビング払いを含むのであれば、特定６ヶ月のみを切り出した支払額、請求額は説明変数に適さないだろう

rawを見て再整理した各項目の説明

default payment next month：翌月滞納 (Yes=1,No=0)  
  
   LIMIT_BAL       ：利用可能枠 (台湾新ドル、家族カードの利用可能枠を含む)  
   SEX             ：性別 (1=male,2=female)  
   EDUCATION       ：最終学歴 (1=graduate school,2=university,3=high school,4=others)  
   MARRIAGE        ：既婚/未婚/その他 (1=married,2=single,3=others)   
   AGE             ：年齢 (year)  
   PAY_Sep-Apr     ：支払歴 (-2=請求なし,-1=遅延なく支払,0=リボルビング払い,1=1ヶ月延滞,...,8=8ヶ月延滞,9=9ヶ月以上延滞)  
   BILL_AMT_Sep-Apr：請求額 (BILL_Sep=請求額(2005年9月),BILL_AMT_Aug=請求額(2005年8月),...BILL_AMT_Apr=請求額(2005年4月))  
   PAY_AMT_Sep-Apr ：支払額 (PAY_AMT_Sep=支払額(2005年9月),PAY_AMT_Aug=支払額(2005年8月),...PAY_AMT_Apr=支払額(2005年4月))

データの特徴を確認する

# データセットの特徴量を見てみる
dframe_e1.describe()

	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_Sep	PAY_Aug	PAY_Jul	PAY_Jun	PAY_May	...	BILL_AMT_Jun	BILL_AMT_May	BILL_AMT_Apr	PAY_AMT_Sep	PAY_AMT_Aug	PAY_AMT_Jul	PAY_AMT_Jun	PAY_AMT_May	PAY_AMT_Apr	DEF_PAY
count	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	30000.000000	...	30000.000000	30000.000000	30000.000000	30000.000000	3.000000e+04	30000.00000	30000.000000	30000.000000	30000.000000	30000.000000
mean	167484.322667	1.603733	1.853133	1.551867	35.485500	-0.016700	-0.133767	-0.166200	-0.220667	-0.266200	...	43262.948967	40311.400967	38871.760400	5663.580500	5.921163e+03	5225.68150	4826.076867	4799.387633	5215.502567	0.221200
std	129747.661567	0.489129	0.790349	0.521970	9.217904	1.123802	1.197186	1.196868	1.169139	1.133187	...	64332.856134	60797.155770	59554.107537	16563.280354	2.304087e+04	17606.96147	15666.159744	15278.305679	17777.465775	0.415062
min	10000.000000	1.000000	0.000000	0.000000	21.000000	-2.000000	-2.000000	-2.000000	-2.000000	-2.000000	...	-170000.000000	-81334.000000	-339603.000000	0.000000	0.000000e+00	0.00000	0.000000	0.000000	0.000000	0.000000
25%	50000.000000	1.000000	1.000000	1.000000	28.000000	-1.000000	-1.000000	-1.000000	-1.000000	-1.000000	...	2326.750000	1763.000000	1256.000000	1000.000000	8.330000e+02	390.00000	296.000000	252.500000	117.750000	0.000000
50%	140000.000000	2.000000	2.000000	2.000000	34.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	19052.000000	18104.500000	17071.000000	2100.000000	2.009000e+03	1800.00000	1500.000000	1500.000000	1500.000000	0.000000
75%	240000.000000	2.000000	2.000000	2.000000	41.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	54506.000000	50190.500000	49198.250000	5006.000000	5.000000e+03	4505.00000	4013.250000	4031.500000	4000.000000	0.000000
max	1000000.000000	2.000000	6.000000	3.000000	79.000000	8.000000	8.000000	8.000000	8.000000	8.000000	...	891586.000000	927171.000000	961664.000000	873552.000000	1.684259e+06	896040.00000	621000.000000	426529.000000	528666.000000	1.000000

8 rows × 24 columns

# 翌月の延滞有無でグループ分けしてみる
dframe_e1.groupby('DEF_PAY').mean()

	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_Sep	PAY_Aug	PAY_Jul	PAY_Jun	PAY_May	...	BILL_AMT_Jul	BILL_AMT_Jun	BILL_AMT_May	BILL_AMT_Apr	PAY_AMT_Sep	PAY_AMT_Aug	PAY_AMT_Jul	PAY_AMT_Jun	PAY_AMT_May	PAY_AMT_Apr
DEF_PAY
0	178099.726074	1.614150	1.841337	1.558637	35.417266	-0.211222	-0.301917	-0.316256	-0.355633	-0.389488	...	47533.365605	43611.165254	40530.445343	39042.268704	6307.337357	6640.465074	5753.496833	5300.529319	5248.220296	5719.371769
1	130109.656420	1.567058	1.894665	1.528029	35.725738	0.668174	0.458258	0.362116	0.254521	0.167872	...	45181.598855	42036.950573	39540.190476	38271.435503	3397.044153	3388.649638	3367.351567	3155.626733	3219.139542	3441.482068

2 rows × 23 columns

グループ分けしてみた結果、翌月の延滞と高い相関関係がありそうなのは、
LIMIT_BAL ：利用可能枠
PAY_0-9 ：支払歴

グループ分けでは、平均値にそれほど違いのなかった項目の分布を見てみる

# 最終学歴分布を見てみる
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
sns.set_style('whitegrid')
sns.set_context("paper")
plt.figure(figsize=(10,6))
sns.countplot('EDUCATION',data=dframe_e1.sort_values(by='EDUCATION'),hue='DEF_PAY',palette='coolwarm')

最終学歴分布

# 利用可能枠の分布を見てみる
plt.figure(figsize=(16,6))
sns.set_style('whitegrid')
sns.set_context("paper")
sns.countplot('AGE',data=dframe_e1.sort_values(by='AGE'),hue='DEF_PAY',palette='coolwarm')

利用可能枠

# 年齢分布を見てみる
plt.figure(figsize=(16,6))
sns.set_style('whitegrid')
sns.set_context("paper")
sns.countplot('LIMIT_BAL',data=dframe_e1.sort_values(by='LIMIT_BAL'),hue='DEF_PAY',palette='coolwarm')

年齢分布

# 説明変数Xを設定
X = DataFrame(dframe_e1, columns=['LIMIT_BAL','SEX','EDUCATION','MARRIAGE','AGE','PAY_Apr','PAY_May','PAY_Jun','PAY_Jul','PAY_Aug','PAY_Sep'])
X.head(5)

	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_Apr	PAY_May	PAY_Jun	PAY_Jul	PAY_Aug	PAY_Sep
0	20000	2	2	1	24	-2	-2	-1	-1	2	2
1	120000	2	2	2	26	2	0	0	0	2	-1
2	90000	2	2	2	34	0	0	0	0	0	0
3	50000	2	2	1	37	0	0	0	0	0	0
4	50000	1	2	1	57	0	0	0	-1	0	-1

# 説明変数Yを設定
Y = dframe_e1['DEF_PAY']
Y.head()

0    1
1    1
2    0
3    0
4    0
Name: DEF_PAY, dtype: int64

# 説明変数は１次元に変換する
Y = Y.values
Y

array([1, 1, 0, ..., 1, 1, 1])

ロジスティック回帰のモデルを作成する

from sklearn.linear_model import LogisticRegression

# LogisticRegressionクラスのインスタンスを作成
log_model1 = LogisticRegression()
# 目的変数、説明変数を使って、モデルを作成。たったのこれだけ。
log_model1.fit(X,Y)
# モデルの精度を確認
log_model1.score(X,Y)

0.77880000000000005

# 'DEF_PAY'の平均値から割り出した、翌月延滞する確率は？
1 - Y.mean()

0.77879999999999994

ロジスティック回帰によって、より高い精度で予測できるとは言えないようだ。

# 変数名とその係数を可視化
coeff_df1 = DataFrame([X.columns, log_model1.coef_[0]]).T
coeff_df1

	0	1
0	LIMIT_BAL	-5.26002e-06
1	SEX	-0.000514679
2	EDUCATION	-0.000644289
3	MARRIAGE	-0.000552566
4	AGE	-0.00907855
5	PAY_Apr	0.000406061
6	PAY_May	0.000428339
7	PAY_Jun	0.000445863
8	PAY_Jul	0.000478645
9	PAY_Aug	0.000547339
10	PAY_Sep	0.00069057

LIMIT_BAL[利用可能枠]は、予測結果に大きく影響する。
（大きくなるほど翌月の延滞の可能性は下がり、小さくなるほど翌月の延滞の可能性は上がる）
支払歴は、最近のものほど翌月の支払に影響を与えているようだ

# LIMIT_BALを除いて、新しい説明変数X2を作成
X2 = DataFrame(X, columns=['SEX','EDUCATION','MARRIAGE','AGE','PAY_Apr','PAY_May','PAY_Jun','PAY_Jul','PAY_Aug','PAY_Sep'])
X2.head(5)

	SEX	EDUCATION	MARRIAGE	AGE	PAY_Apr	PAY_May	PAY_Jun	PAY_Jul	PAY_Aug	PAY_Sep
0	2	2	1	24	-2	-2	-1	-1	2	2
1	2	2	2	26	2	0	0	0	2	-1
2	2	2	2	34	0	0	0	0	0	0
3	2	2	1	37	0	0	0	0	0	0
4	1	2	1	57	0	0	0	-1	0	-1

# 新しいモデルの作成
log_model2 = LogisticRegression()
log_model2.fit(X2,Y)
# 新しいモデルの精度を確認
log_model2.score(X2,Y)

0.80983333333333329

精度が向上した！

# 変数名とその係数を可視化
coeff_df2 = DataFrame([X2.columns, log_model2.coef_[0]]).T
coeff_df2

	0	1
0	SEX	-0.104595
1	EDUCATION	-0.055019
2	MARRIAGE	-0.134139
3	AGE	0.00447237
4	PAY_Apr	0.0015857
5	PAY_May	0.0348254
6	PAY_Jun	0.0211476
7	PAY_Jul	0.0868766
8	PAY_Aug	0.0852904
9	PAY_Sep	0.613086

# 今度は支払歴だけで予測してみる
X3 = DataFrame(X, columns=['PAY_Apr','PAY_May','PAY_Jun','PAY_Jul','PAY_Aug','PAY_Sep'])
X3.head(5)

	PAY_Apr	PAY_May	PAY_Jun	PAY_Jul	PAY_Aug	PAY_Sep
0	-2	-2	-1	-1	2	2
1	2	0	0	0	2	-1
2	0	0	0	0	0	0
3	0	0	0	0	0	0
4	0	0	0	-1	0	-1

# 新しいモデルの作成
log_model3 = LogisticRegression()
log_model3.fit(X3,Y)
# 新しいモデルの精度を確認
log_model3.score(X3,Y)

0.80876666666666663

若干精度は下がった

# より精度の高かったモデルで、直近３ヶ月の支払履歴に絞って予測して見る
X4 = DataFrame(X, columns=['SEX','EDUCATION','MARRIAGE','AGE','PAY_Jul','PAY_Aug','PAY_Sep'])
X4.head(5)

	SEX	EDUCATION	MARRIAGE	AGE	PAY_Jul	PAY_Aug	PAY_Sep
0	2	2	1	24	-1	2	2
1	2	2	2	26	0	2	-1
2	2	2	2	34	0	0	0
3	2	2	1	37	0	0	0
4	1	2	1	57	-1	0	-1

# 新しいモデルの作成
log_model4 = LogisticRegression()
log_model4.fit(X4,Y)
# 新しいモデルの精度を確認
log_model4.score(X4,Y)

0.80936666666666668

試した中で、ベストな目的変数は下記だった
‘SEX’,‘EDUCATION’,‘MARRIAGE’,‘AGE’,‘PAY_Apr’,‘PAY_May’,‘PAY_Jun’,‘PAY_Jul’,‘PAY_Aug’,‘PAY_Sep’

データを学習用とテスト用に分けて、モデルの性能を確認

from sklearn.cross_validation import train_test_split
from sklearn import metrics
# train_test_splitを使う
X_train, X_test, Y_train, Y_test = train_test_split(X2, Y)
# 新しいモデルを作成
log_model5 = LogisticRegression()
# 学習用のデータだけ学習
log_model5.fit(X_train, Y_train)

LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)

# テスト用データを使って予測
class_predict = log_model5.predict(X_test)
# 予測の精度を確認
print(metrics.accuracy_score(Y_test,class_predict))

0.814

全データを使用した時と同程度の精度を出すことができた

感想

• 決定木分析に比べてすごく扱いやすく、わかりやすい印象
• 説明変数の設定は、どのように試行して検証するのが良いのかを理解していく必要がある
• データの収集方法として、目的変数を得る前に説明変数を収集する(顧客の属性や行動履歴)のと、目的変数と一緒に説明変数を収集する(アンケート)のとでは、なんとなくだがロジスティック回帰の有効性に影響がありそうな気がする

参考にしたサイト

*1:こちらのサイトを参考にしました

ロジスティック回帰 :: 株式会社アイスタット｜統計分析研究所