評分卡系列（二）：特征工程

時間：JSong 時間：2018.01.14

文章很長，理論和實現都講的很細，大家可以先收藏，有時間再看。

在上一篇文章中，我們對LendingClub的數據有了一個大致的了解，這次我將帶大家把10萬多條、145個字段的原始數據一步一步處理成建模所需輸入的數據。

我們先按照上次一樣導入數據，這里我將逾期15天以上的都當作正類

import pandas as pd
import numpy as np
import reportgen as rpt
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

# 字段的中文對照表，暫時只翻譯了一部分
datadict=pd.read_excel('.\\LendingClubData\\LCDataDictionary.xlsx')
columnnew=dict(zip(datadict.loc[datadict['中文名稱'].notnull(),'LoanStatNew'],datadict.loc[datadict['中文名稱'].notnull(),'中文名稱']))
data=pd.read_csv('.\\LendingClubData\\LoanStats_2016Q4.csv',skiprows=1)
data=data.rename(columns=columnnew)
data['target']=data['target'].replace({'Current':np.nan,'Fully Paid':0,'Charged Off':1,'Late (31-120 days)':1,'In Grace Period':np.nan,\
                                      'Late (16-30 days)':1,'Default':np.nan})
data=data.loc[data['target'].notnull(),:]

1、評分卡簡介

在進行下一步操作之前，我們先來解構一下評分卡。

貸款機構（含銀行、信用卡、互聯網金融等）的主要目標之一是讓貸款產生更多的利潤，由於違約的存在，這里就需要在市場份額和利潤之間做一個平衡。評分卡本質是在量化借款人的違約可能性，通過評分卡，我們能更好的把控這個平衡。

一般我們把逾期 n 天以上的稱為壞人，反之是好人。貸款機構在全面掌握信息后，對借款人按時和按約定還款的可能性作出有洞見的預測。我們希望這些預測有更多區分度，而不是好人可能性高或低這么簡單（每個機構的策略都不同，或保守或偏激，所以需要預測結果有足夠區分度）。

令 X 是所有的特征數據（先驗知識），G 是好人，B是壞人。給定相應數據 x ,有條件概率 p(G|x) 和 p(B|x)表示給定特征數據下好人和壞人的概率。

\[p(G|x)+p(B|x)=1,\,\,\,x \in X \]

有時候我們更多的考慮事件的發生比率 ：

\[o(G|x)=\frac{p(G|x)}{p(B|x)}\quad\,or\quad\,o(B|x)=\frac{p(B|x)}{p(G|x)} \]

由 Bayes 定理，我們可以得到:

\[p(G|x)=\frac{f(x|G)p_{G}}{f(x)},\quad p(B|x)=\frac{f(x|B)p_{B}}{f(x)} \]

其中f(x)表示申請者具有屬性向量 x 的概率，p_G 和 p_B 表示先驗知識中好人和壞人的概率, f(x|G) 和 f(x|B) 被稱為似然函數，描述屬性向量有多大可能性落在好和壞的群體中。把上述公式代入發生比率，還可得：

\[o(G|x)=\frac{p(G|x)}{p(B|x)}=\frac{p_G}{p_B}\frac{f(x|G)}{f(x|B)}\triangleq o_{pop} \times I(x) \]

其中 o_{pop}=p_G/p_B 稱為總體發生比率（population odds）,I(x)稱為信息比率，當其大於1時，表明屬性x的借款人更可能是好人。

上述表達式其實就是Bayes 評分卡構建的所有理論了，其假設條件弱，應用范圍廣，但准確率不高。接下來我們會介紹邏輯回歸等方法，這些方法只支持數值變量，不支持因子變量等，於是在建模之前，除了特征篩選外，我們還需要對數據進行編碼。

2、特征工程

我把整個特征工程分為四個步驟

1. 數據清洗(預處理)
2. 特征衍生
3. 特征編碼
4. 特征篩選

2.1 數據清洗

因為業務關系，我這里的預處理操作主要來源於知乎專欄 《風險狗的數據分析之路》。主要的預處理如下：

1、貸后的相關變量除了target變量，其余直接剔除，因為貸后表現在客戶申請時是沒有的，如果進入模型實際上就成了未來變量；2、缺失率太高的變量直接剔除，本文是按65%的閾值來剔除的；3、數值變量中所有值方差太小接近常量的變量剔除，因為不能提供更多信息；4、按業務邏輯完全不可解釋的變量直接剔除，5、分類變量中unique值大於20的直接剔除。

先剔除與建模無關的變量（待最后一步來操作）

# 實際執行順序放在特征衍生的后面
var_drop=list(set(data.columns) & set(datadict.loc[datadict['類別'].isnull(),'LoanStatNew']))
data=data.drop(var_drop,axis=1)

剔除缺失率過高的變量

# 去除缺失率大於65%的字段
thr=0.65
missing_pct=data.apply(lambda x : (len(x)-x.count())/len(x))
data=data.loc[:,missing_pct[missing_pct<thr].index]
# 另也可以直接用下面這種方法
#data=data.dropna(thresh=len(data) * thr,axis=1)

剔除unique值過少或過多的變量

# 每個字段的unique數目
nunique=data.apply(pd.Series.nunique)
data=data.loc[:,nunique!=1]
#log=[is_numeric_dtype(data[c].dtype) or nunique[c]<20 for c in data.columns]
#data=data.loc[:,log]

缺失值處理

涉及到實現上的一些處理，我們放在第三步中處理

# 實際執行順序在 WOE 編碼前
# 缺失值處理
# 對於連續變量，這里我們隨意填入
for v in continuous_var:
    s=list(data[v].dropna().unique())
    data[v]=data[v].map(lambda x : np.random.choice(s,1)[0] if pd.isnull(x) else x)
# 也可以直接用
#from sklearn.preprocessing import Imputer
#data.loc[:,continuous_var]=Imputer(strategy='mean').fit_transform(data.loc[:,continuous_var])
# 對於因子變量，這里我們填入眾數
for v in categorical_var:
    s=data[v].value_counts().argmax()
    data[v]=data[v].fillna(s)

data=data.reset_index(drop=True)#便於后續數據處理，不解釋
data0=data.copy()#備份用

無量綱化

• 標准化：

\[\frac{x-\bar{x}}{s} \]

• 極差化

\[\frac{x-x_{min}}{x_{max}-x_{min}} \]

# 此處只作為演示用，不實際操作
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
tmp=StandardScaler().fit_transform(data.select_dtypes(include=['float64']))
tmp=MinMaxScaler().fit_transform(data.select_dtypes(include=['float64']))

2.2 特征衍生

data['放款日期']=pd.to_datetime(data['放款日期'])
data['信用報告最早日期']=pd.to_datetime(data['信用報告最早日期'])
data['days']=(data['放款日期']-data['信用報告最早日期']).map(lambda x:x.days)
data=data.drop(['放款日期','信用報告最早日期'],axis=1)

2.3 特征編碼的理論介紹

我把特征的類型分為五種：

類型	名稱	描述
number	連續變量	數值型，例如收入、貸款金額等
category	因子變量	數值或者字符串，例如學歷、年齡、基礎等級等
datetime	時間變量	日期和時間，例如放款日期等
text	文本變量	字符串，例如產品評論等
text_st	結構化文本變量	例如樣本的ID等，需要剔除或者再處理的變量

自己寫了一個函數 type_of_var ，放在工具箱 reportgen 中用來識別變量的類型

# 此處category_detection 設為False，所有數值類型全部當成連續變量
VarType=rpt.analysis.type_of_var(data,category_detection=False)
VarType=pd.Series(VarType)
continuous_var=list(VarType[VarType=='number'].index)
continuous_var.remove('target')
print('連續變量有：\n',' , '.join(continuous_var))
categorical_var=list(VarType[VarType=='category'].index)
print('因子變量有：\n',' , '.join(categorical_var))
datetime_var=list(VarType[VarType=='datetime'].index)
print('時間變量有：\n',' , '.join(datetime_var))
text_var=list(VarType[VarType=='text'].index)
print('文本變量有：\n',' , '.join(text_var))

輸出如下：

>連續變量有：
>bc_open_to_buy , bc_util , days , mo_sin_old_il_acct , mo_sin_old_rev_tl_op , mo_sin_rcnt_tl , mths_since_last_delinq , mths_since_recent_bc , num_accts_ever_120_pd , num_actv_rev_tl , num_bc_sats , num_il_tl , num_rev_accts , num_rev_tl_bal_gt_0 , num_sats , num_tl_30dpd , num_tl_90g_dpd_24m , num_tl_op_past_12m , pct_tl_nvr_dlq , percent_bc_gt_75 , tot_cur_bal , tot_hi_cred_lim , total_bal_ex_mort , total_bc_limit , total_rev_hi_lim , 沖銷數量 , 年收入 , 總貸款筆數 , 總貸款金額 , 月負債比 , 過去24個月的交易數目。 , 過去6個月內被查詢次數 , 銀行卡帳號數目
>因子變量有：
>房屋所有權

特征編碼的對象包含因子變量和連續變量，對於連續變量，主要有如下幾種方式

2.3.1 標准化

就是上文講的無量綱化

2.3.2 二值化

首先設定一個閾值，大於閾值的設為1，否則設為0

# 不做演示，不實際執行
from sklearn.preprocessing import Binarizer
Binarizer_var=[]
Binarizer(threshold=3).fit_transform(data.loc[:,Binarizer_var])

2.3.3 特征分箱

特征分箱包含將連續變量離散化和將多狀態的離散變量合並成少狀態。

這里我們只介紹卡方分箱（ChiMerge）

自底向上的(即基於合並的)數據離散化方法。它依賴於卡方檢驗:具有最小卡方值的相鄰區間合並在一起,直到滿足確定的停止准則。

基本思想:對於精確的離散化，相對類頻率在一個區間內應當完全一致。因此,如果兩個相鄰的區間具有非常類似的類分布，則這兩個區間可以合並；否則，它們應當保持分開。而低卡方值表明它們具有相似的類分布。

詳細的分箱算法如下：

它們的Python實現方法如下：

# 只做演示，不實際執行
# 卡方分箱
from reportgen.utils import Discretization
dis=Discretization(method='chimerge',max_intervals=20)
dis.fit_transform(x,y)
# 等距分箱
pd.cut(x,bins=10)
# 等頻分箱
pd.qcut(x,bins=10)

對於因子變量，主要有如下幾種方式

2.3.4 LabelEncoder

這個比較簡單，就是直接將整數來替換類別，例如在性別中用 1 替換男，用 2 替換女

# 只做演示，不實際執行
from sklearn.preprocessing import LabelEncoder
LabelEncoder().fit_transform(data.loc[:,categorical_var])

2.3.5 啞編碼

啞編碼是一種狀態編碼，屬於一對多的特征映射。簡單點講，它將性別映射為兩個變量：是否是男性、是否是女性。它解決了 LabelEncoder 中序的問題，比如在 LabelEncoder 中，女性用2表示，但明顯不可能是2倍的男性。這里提供兩種啞編碼的實現方法，pandas和sklearn。它們最大的區別是，pandas默認只處理字符串類別變量，sklearn默認只處理數值型類別變量(需要先 LabelEncoder )

# 只做演示，不實際執行
# pandas 的方法，drop_first 是為了防止共線性
pd.get_dummies(data.loc[:,categorical_var],prefix=categorical_var,drop_first=True)
# sklearn 的方法
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import LabelEncoder
Enc_ohe, Enc_label = OneHotEncoder(), LabelEncoder()
Enc_ohe.fit_transform(Enc_label.fit_transform(data.loc[:,categorical_var]))

2.3.6 WOE編碼

在第一節中我們介紹了 總體發生比率 和 信息比率：

\[o_{opp}=\frac{p_G}{p_B},\quad I(x)=\frac{f(x|G)}{f(x|B)} \]

前者反映了還沒有任何關於借款人的已知信息時，我們對該借款人是好人的可能性認知。后者很好的度量了借款人信息對好人可能性的貢獻程度，其自然對數 ln(I(x)) 是評估向量x攜帶信息的一種有效途徑，我們將這個數值稱之為 x 提供的證據權重（weights of evidence，WOE）

\[w(x)=\ln\ I(x)=\ln\big(\frac{f(x|G)}{f(x|B)}\big) \]

如果一個特征有K個類別，且用\(g_k\)和\(b_k\)表示第k類中好人和壞人的數量，用\(n_G\)和\(n_B\)表示好人和壞人的數量，則 WOE 可以表示為：

\[WOE_k=\ln\frac{g_kn_B}{b_kn_G},\quad \forall k \in K \]

WOE的值越大代表對應的變量對“是好人”的貢獻就越大，反之，越小就代表對應的變量對“是壞人”的貢獻越大。所以WOE值可以作為特征的一種編碼方式。

from reportgen.utils import WeightOfEvidence
woe_iv={}# 用於之后的特征選擇
for v in categorical_var:
    woe=WeightOfEvidence()
    woe.fit(data[v],data['target'])
    data[v]=woe.transform(data[v])

2.4 特征編碼

LendingClub 的數據類型很多，理論上講它的特征編碼有很多種組合方式，比如

1. 因子變量啞變量,連續變量標准化
2. 因子變量WOE編碼，連續變量標准化
3. 連續變量離散后再WOE編碼，同時因子變量 WOE 編碼

不過為了節省篇幅，這里我們采用最后一種方式，全部 WOE 編碼。其代碼如下

from reportgen.utils import WeightOfEvidence 
from reportgen.utils import Discretization 
max_intervals=8
N=len(data)
woe_iv={}# 用於之后的特征選擇
for v in categorical_var:
    # 如果因子數過多，則先進行分箱
    if not(isinstance(data[v].iloc[0],str)) and len(data[v].unique())>20:       
        dis=Discretization(method='chimerge',max_intervals=max_intervals)
        dis.fit(data[v],data['target'])
        data.loc[:,v]=dis.transform(data[v])
    woe=WeightOfEvidence()
    woe.fit(data[v],data['target'])
    woe_iv[v]={'woe':woe.woe,'iv':woe.iv}
    data[v]=data[v].replace(woe.woe).astype(np.float64)

for v in continuous_var:
    if len(data[v].unique())>max_intervals:
        dis=Discretization(method='chimerge',max_intervals=max_intervals,sample=1000)
        dis.fit(data[v],data['target'])
        data.loc[:,v]=dis.transform(data[v])
    woe=WeightOfEvidence()
    woe.fit(data[v],data['target'])
    woe_iv[v]={'woe':woe.woe,'iv':woe.iv}
    data[v]=data[v].replace(woe.woe).astype(np.float64)

2.5 特征選擇的理論介紹

特征選擇的標准是該特征與目標的相關性，與目標相關性高的特征，應當優選擇。按照選擇形式可以將特征選擇的方法分為三種：

• 過濾法，按照相關性對各個特征進行評分，設定閾值或者待選擇閾值的個數，選擇特征。
• 包裝法，根據目標函數（通常是預測效果評分），每次選擇若干特征，或者排除若干特征。
• 嵌入法，先使用某些機器學習的算法和模型進行訓練，得到各個特征的權值系數，根據系數從大到小選擇特征。類似於Filter方法，但是是通過訓練來確定特征的優劣。

除此之外，還有一個需要處理的是共線性問題。

2.5.1 Filter：卡方統計量

卡方檢驗常用於兩個變量之間的顯著性檢驗，假定fo、fe分別為觀察頻數和期望頻數，則卡方統計量的計算公式為：

\[\mathcal{X}^2=\sum_{i,j}\frac{(fo-fe)^2}{fe} \]

當我們計算了所有變量的卡方統計量后，可以用p值來篩選變量，也可以用衍生的V相關系數來篩選：

\[V=\sqrt{\frac{\mathcal{X}^2}{n\times\min[(R-1),(C-1)]}} \]

其中R代表列聯表的行數，C代表列聯表的列數。

from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2
#from reportgen.utils import chi2

#選擇K個最好的特征，返回選擇特征后的數據
# sklearn 的chi2只支持兩個類別的特征，建議換成reportgen中的chi2
SelectKBest(chi2, k=2).fit_transform(X, y)

2.5.2 Filter：信息量（Info Value, IV）

如果想考察某個特征區分好壞借款人的表現，我們可以用該特征的均值之差來表示

\[\int(f(x|G)-f(x|B))xdx \]

然而這個差並沒有考慮到某些x值的信息量遠高於其他的情況，於是我們可以用權重之差來判斷

\[\int(f(x|G)-f(x|B))w(x)dx \]

這被稱為散度，也類似於相對熵（進行了對稱處理）。將散度離散化便得到信息量（IV） . 如果一個特征有K個類別，且用\(g_k\)和\(b_k\)表示第k類中好人和壞人的數量，用\(n_G\)和\(n_B\)表示好人和壞人的數量，則 IV 可以表示為：

\[IV=\sum_{k=1}^{K}[\frac{g_k}{n_G}-\frac{b_k}{n_B}]\ln[\frac{g_{k}n_B}{b_kn_G}] \]

一般IV值越大，該特征越要保留。

from reportgen.utils import WeightOfEvidence
woe=WeightOfEvidence()
woe.fit(x,y)
iv=woe.iv

注：卡方統計量和熵的關系

注意到卡方統計量等價於：

\[\mathcal{X}^2=\frac{(fo-fe)^2}{fe}=(fo-fe)\times(\frac{fo}{fe}-1) \]

而 fo/N 可以看成是聯合分布的概率，fe/N 可以看成是乘積分布的概率，又

\[fo/fe-1\approx\log(fo/fe) \]

所有我們有

\[\frac{\mathcal{X}^2}{N}=\int[p(x,y)-p(x)p(y)]\times[\frac{p(x,y)}{p(x)p(y)}-1]dxdy\approx\,D \]

從這個角度來看，卡方距離是 p(x,y) 和 p(x)p(y) 之間信息量的某種近似。

2.5.3 Wrapper：遞歸特征消除法

遞歸消除特征法使用一個基模型來進行多輪訓練，每輪訓練后，消除若干權值系數的特征，再基於新的特征集進行下一輪訓練。使用feature_selection庫的RFE類來選擇特征的代碼如下：

from sklearn.feature_selection import RFE
from sklearn.linear_model import LogisticRegression
rfe=RFE(estimator=LogisticRegression(),n_features_to_select=2)
rfe.fit_transform(X, y)

2.5.4 Embedded: 基於分類模型的特征選擇法

使用基模型，除了篩選出特征外，同時也進行了降維。這里我們用sklearn中的SelectFromModel來完成。

from sklearn.feature_selection import SelectFromModel
from sklearn.linear_model import LogisticRegression

#帶L1懲罰項的邏輯回歸作為基模型的特征選擇
sf=SelectFromModel(LogisticRegression(penalty="l1",C=0.1))
sf.fit_transform(X, y)

#樹模型中GBDT也可用來作為基模型進行特征選擇，
from sklearn.ensemble import GradientBoostingClassifier
sf=SelectFromModel(GradientBoostingClassifier())
sf.fit_transform(X,y)

2.5.5 共線性問題

詳見我的文章

2.6 特征選擇

由於我們上一步采用的是 WOE 編碼，所以這里選用信息量 IV 作為准則。

iv_threshould=0.005
varByIV=[k for k,v in woe_iv.items() if v['iv'] > iv_threshould]
print('剩余%d個變量'%len(varByIV))
print(varByIV)

X=data.loc[:,varByIV]
y=data['target']

輸出結果如下：

剩余28個變量
['月負債比', 'tot_hi_cred_lim', 'bc_util', 'mths_since_last_delinq', '年收入', 'num_actv_rev_tl', 'bc_open_to_buy', 'total_rev_hi_lim', 'num_tl_90g_dpd_24m', '銀行卡帳號數目', 'num_rev_tl_bal_gt_0', 'tot_cur_bal', '沖銷數量', 'mo_sin_old_rev_tl_op', 'total_bc_limit', 'percent_bc_gt_75', 'pct_tl_nvr_dlq', '房屋所有權', '總貸款金額', 'num_rev_accts', '過去24個月的交易數目。', 'num_tl_op_past_12m', '過去6個月內被查詢次數', 'mo_sin_rcnt_tl', 'mths_since_recent_bc', '總貸款筆數', 'mo_sin_old_il_acct', 'days']

至於共線性問題，簡單點可以線看一下條件數，用於檢測是否存在共線性，然后再用嶺回歸來確定需要剔除的變量

print('條件數 = %.3f'%np.linalg.cond(X))
# 輸出為：條件數 = 8.682

條件數很小，一般至少要大於100才說明存在共線性。如果不放心的話，還可以看看嶺回歸的嶺跡：

from sklearn.linear_model import Ridge
plt.figure()
n_alphas = 20
alphas = np.logspace(-1,4,num=n_alphas)
coefs = []
for a in alphas:
    ridge = Ridge(alpha=a, fit_intercept=False)
    ridge.fit(X2, y)
    coefs.append(ridge.coef_)
ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
handles, labels = ax.get_legend_handles_labels()
plt.legend(labels=labels)
plt.show()

可以看出所有的特征都很穩定，不存在波動。

3、建模

本文對模型部分不做詳細介紹，都用默認參數。詳細可參見 JSong 的公眾號文章一切從邏輯回歸開始

from sklearn.cross_validation import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import auc,roc_curve

X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3)

lr=LogisticRegression()
lr.fit(X_train,y_train)
y_train_proba=lr.predict_proba(X_train)
y_train_label=lr.predict(X_train)
y_test_proba=lr.predict_proba(X_test)
y_test_label=lr.predict(X_test)

from sklearn.metrics import accuracy_score
print('訓練集准確率：{:.2%}'.format(accuracy_score(y_train_label,y_train)))
print('測試集准確率：{:.2%}'.format(accuracy_score(y_test_label,y_test)))

# ROC曲線和KS統計量
fpr, tpr, thresholds=roc_curve(y_train,y_train_proba[:,1],pos_label=1)
auc_score=auc(fpr,tpr)
w=tpr-fpr
ks_score=w.max()
ks_x=fpr[w.argmax()]
ks_y=tpr[w.argmax()]
fig,ax=plt.subplots()
ax.plot(fpr,tpr,label='AUC=%.5f'%auc_score)
ax.set_title('Receiver Operating Characteristic')
ax.plot([0, 1], [0, 1], '--', color=(0.6, 0.6, 0.6))
ax.plot([ks_x,ks_x], [ks_x,ks_y], '--', color='red')
ax.text(ks_x,(ks_x+ks_y)/2,'  KS=%.5f'%ks_score)
ax.legend()
fig.show()


#密度函數和KL距離
kl_div_score1=rpt.utils.entropyc.kl_div(y_train_proba[y_train==0,0],y_train_proba[y_train==1,0])
kl_div_score2=rpt.utils.entropyc.kl_div(y_train_proba[y_train==1,0],y_train_proba[y_train==0,0])
kl_div_score=kl_div_score1+kl_div_score2
fig,ax=plt.subplots()
sns.distplot(y_train_proba[y_train==0,0],ax=ax,label='good')
sns.distplot(y_train_proba[y_train==1,0],ax=ax,label='bad')
ax.set_title('KL_DIV = %.5f'%kl_div_score)
ax.legend()
fig.show()

輸出結果如下

訓練集准確率：72.44%
測試集准確率：72.51%

可以看到結果不太好，沒事，后面我們會慢慢提高它的。本文的代碼，如果不想一段一段復制，也可以后台回復：數據集 下載。

評分卡文章系列：

1、評分卡系列（一）：講講評分系統的構建

2、評分卡系列（二）：特征工程

3、評分卡系列（三）：分類學習期的評估

---

個人公共號，文章首發於此。轉載請聯系。