五、逻辑回归

统计模型综述

变量		预测变量
		分类	连续	连续和分类
反应变量	连续	方差分析	普通最小二乘法	协方差分析
	分类	列联表分析、逻辑回归	逻辑回归	逻辑回归

1. 引言——汽车贷款信用违约预测案例

- 1.导入数据：汽车金融数据
##bad_ind---是否违约
##vehicle_year---汽车购买时间
##vehicle_make---汽车制造商
##bankruptcy_ind---曾经破产标识
##tot_derog---五年内信用不良事件数量(比如手机欠费消号)
##tot_tr---全部帐户数量
##age_oldest_tr---最久账号存续时间(月)
##tot_open_tr---在使用帐户数量
##tot_rev_tr---在使用可循环贷款帐户数量(比如信用卡)
##tot_rev_debt---在使用可循环贷款帐户余额(比如信用卡欠款)
##tot_rev_line---可循环贷款帐户限额(信用卡授权额度)
##rev_util---可循环贷款帐户使用比例(余额/限额)
##fico_score---FICO打分
##purch_price---汽车购买金额(元)
##msrp---建议售价
##down_pyt---分期付款的首次交款
##loan_term---贷款期限(月)
##loan_amt---贷款金额
##ltv---贷款金额/建议售价*100
##tot_income---月均收入(元)
##veh_mileage---行使历程(Mile)
##used_ind---是否二手车
##weight---样本权重
import os
import numpy as np
from scipy import stats
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

#pd.set_option('display.max_columns', None)
os.chdir(r"E:\脚本")
accepts = pd.read_csv(r'accepts.csv').dropna()
accepts

1.1 分析衍生变量

• 拿到数据之后，可以考虑是否做衍生变量。

# 定义除法函数
def divMy(x,y):
    import numpy as np
    if x==np.nan or y==np.nan:
        return np.nan
    elif y==0:
        return -1
    else:
        return float(x)/float(y)

- 衍生变量
- 1.历史负债收入比:tot_rev_line/tot_income
accepts["dti_hist"]=accepts[["tot_rev_line","tot_income"]].apply(lambda x:divMy(x[0],x[1]),axis = 1)
- 2.本次新增负债收入比:loan_amt/tot_income
accepts["dti_mew"]=accepts[["loan_amt","tot_income"]].apply(lambda x:divMy(x[0],x[1]),axis = 1)
- 3.本次贷款首付比例:down_pyt/loan_amt
accepts["fta"]=accepts[["down_pyt","loan_amt"]].apply(lambda x:divMy(x[0],x[1]),axis = 1)
- 4.新增债务比:loan_amt/tot_rev_debt
accepts["nth"]=accepts[["loan_amt","tot_rev_debt"]].apply(lambda x:divMy(x[0],x[1]),axis = 1)
- 5.新增债务额度比:loan_amt/tot_rev_line
accepts["nta"]=accepts[["loan_amt","tot_rev_line"]].apply(lambda x:divMy(x[0],x[1]),axis = 1)

1.2 判断是否相关

Y(是否违约)是二分类变量，判断是否相关：计算频次并且频次求行百分比，若百分比相差不大则无相关性。

步骤一：
原假设：是否破产与是否违约没有关系
备择假设：是否破产与是否违约是有相关关系的。
步骤二：
判断相关性，分类变量计算行百分比
步骤三：
检验相关性：分类变量采用卡方检验

1.2.1 分类变量相关性：交叉表

• Y：是否违约；X：是否二手车
  行变量：解释变量(X)；列变量：被解释变量(Y)

- 1.分类变量的相关关系：交叉表 
cross_table = pd.crosstab(accepts.used_ind,accepts.bad_ind,margins=True)
cross_table

• 2.计算行百分比，判断相关性

- 2.计算行百分比
def percConvert(ser):
    return ser/float(ser[-1])

cross_table.apply(percConvert, axis=1)

百分比相差不大，认为是相互独立的。

• 3.数据独立满足的公式
  X、Y独立，则满足$P(X=0,Y=0)=P(X=0)\ast P(Y=0)$

1.3 检验相关性

X、Y是分类变量，所以用卡方检验来检验相关性。公式即：
$${\chi }^2=\sum _{i=1}^R\sum _{j=1}^C\frac{(Ob{s}_{i}j-Ex{p}_{i}j{)}^2}{Ex{p}_{i}j}$$
卡方值 = (样本量的观测 - 原假设成立时的期望)^2/期望值

- 4.卡方检验
print(''' chisq = %6.4f p-value=%6.4f dof=%i excepted_freq=%s''' %stats.chi2_contingency(cross_table.iloc[:2,:2]))

模型解释：
原假设是X与Y之间是独立的，样本量是8211，则对应的α是【0.001 - 0.1】，P值<α表示原假设不成立，P值>α表示原假设是成立的。结果可知，P值为0.0004，所以是有关系的。

2.逻辑回归

案例：最久账号存续时间和违约有什么影响？
反应变量(X)：违约
预测变量(Y)：最久账号存续时间

2.1 构建假设

原假设：借贷时间与是否违约没有关系
备择假设：借贷事件与是否违约有关系

2.2 绘制散点图

accepts.plot(x='age_oldest_tr',y='bad_ind',kind='scatter')

模型解释：
账户存续时间越长，不违约的人要多一点，违约的人相对少一点。

2.3 选择回归模型

2.3.1 普通最小二乘法模型(OLS)

OLS回归：$Y_i={\beta }_0+{\beta }_1\ast X_{1i}+{ϵ}_i$
Y的取值可能为[0,1]之间的小数，也可能为>1或<0的数。而原数据只有0和1两种情况，不符常理。

2.3.2 线性回归模型

若构建线性回归，$$P_i={\beta }_0+{\beta }_1\ast X_{1i}$$
P(Y = 1)的范围是[0,1]，而X变量的范围(-∞，+∞)，等号两侧是不匹配的。

2.3.3 逻辑回归模型

2.3.3.1优势比(Odd Ratios)

优势比 = 事件发生的概率/事件不发生的概率，即$$Odds=\frac{P_{event}}{1-P_{event}}$$
注：
原本P(概率)的范围为(0,1)，经过优势比(Odds Ratios),P的范围变为(-∞，+∞)，Odds的值域[0,+∞)。

2.3.3.2 逻辑回归公式

逻辑回归公式1：
$$logit(p_i)={\beta }_0+{\beta }_1\ast X_{1i}+...+{\beta }_k\ast X_{ki}$$
$logit(p_i)$：事件发生概率的logit
${\beta }_0$：回归式的截距
${\beta }_k$：第k个预测变量的参数估计

逻辑回归公式2：
$$logit(p_i)=ln(\frac{p_i}{1-p_i})$$
logit：是Odds的自然对数
i：是所有案例(观察值)
$p_i$：在第i个案例中一个事件发生的概率
ln：自然对数(底数为e)

逻辑回归总公式：
$$ln(\frac{p_i}{1-p_i})={\beta }_0+{\beta }_1\ast X_{1i}+...+{\beta }_k\ast X_{ki}$$

注：逻辑回归可以反求出概率

3.逻辑回归建模

3.1 划分训练集和测试集

- 1.数据分为训练集(70%)和测试集(30%)
train = accepts.sample(frac=0.7, random_state=1234).copy()
test = accepts[~ accepts.index.isin(train.index)].copy()  #~表示否定
print(' 训练集样本量: %i \n 测试集样本量: %i' %(len(train), len(test)))

3.2 建立回归模型

根据下图，回归模型中红线上面的内容没有实际意义，红线下面是变量每个水平的估计系数，这里的“估计”等同于回归结果。

- 2.建立回归模型
import statsmodels.formula.api as smf
# family:族，伯努利分布族(0-1分布)
# 连接函数(link)为logit或probit
lg = smf.glm('bad_ind ~ age_oldest_tr', data=train, 
             family=sm.families.Binomial(sm.families.links.logit)).fit()
lg.summary()

模型解释：
intercept：截距项
age_oldest_tr(即X)：只跟logit有关，历史账户存续时间每增加一个单位，违约的logit会下降0.0057个单位

3.3 模型进行预测

根据公式$P=\frac{1}{1+e^{-({\beta }_0+{\beta }_1\ast X)}}$，可以预测得到P值(每个样本违约的预测)。

- 3.作预测
# 基于训练集的预测
train['proba'] = lg.predict(train)
# 基于测试集的预测
test['proba'] = lg.predict(test)
# 测试集的预测
test['proba'].head(10)

4. 模型评估

模型评估最重要的指标：成对比较

4.1 成对比较

• 1.一致对：模型预测和实际结果一致

在最久账号存续时间和客户履约正相关，客户A履约20个月后违约，客户B履约36个月，账户时长越长，违约概率越小，履约概率越大。

• 2.不一致对：模型预测和实际结果不一致
  客户A履约20个月后违约，客户B 履约10个月，得到的是A履约的概率大于B履约的概率。
  

• 3.相等对：X的取值一样，Y的取值不一样
  客户A履约1个月，客户B履约1个月后违约，两客户履约的概率相等。
  

• 4.成对比较的结果(SAS提供)：
  
  c的结果很重要，因为python的statmodels中不提供，所以采用ROC曲线完成。

4.2 ROC曲线面积

AUC = C统计量(SAS)

• 1.绘制ROC曲线

import sklearn.metrics as metrics
#测试集中真实的Y值和测试集中Y的概率，并返回三个值
fpr_test, tpr_test, th_test = metrics.roc_curve(test.bad_ind, test.proba)

# 训练集真实的Y值和训练集Y的概率作散点图
fpr_train, tpr_train, th_train = metrics.roc_curve(train.bad_ind, train.proba)
plt.figure(figsize=[3, 3])
plt.plot(fpr_test, tpr_test, 'b--')
plt.plot(fpr_train, tpr_train, 'r-')
plt.title('ROC curve')
plt.show()

• 2.计算曲线面积

- AUC面积范围[0.5,1],面积越大越好
print('AUC = %.4f' %metrics.auc(fpr_test,tpr_test))

Out[]: AUC = 0.6541

4.3 混淆矩阵

混淆矩阵：给定一个阈值，就可以做出一个混淆矩阵

混淆矩阵		打分值		合计
		反应(预测=1)	未反应(预测=0)
真实结果	呈现信号(真实=1)	A(击中)	B(漏报)	A+B
	未呈现信号(真实=0)	C(虚报)	D(正确否定)	C+D
合计		A+C	B+D	A+B+C+D

• 正确率 = (A+D)/(A+B+C+D)：只有样本模型是1：1时才可用
• 灵敏度(recall) = A/(A+B)：实际违约人中，模型预测出了多少
• 命中率(PV+) = A/(A+C)：预测违约的人中，有多少人是真正的违约
• 特异度 = D/(C+D)：特异度越高越好
• 负命中率 = D/(D+B)

步骤一：预测情况

- 可以根据Y的预测情况，大于0.3的为1，小于0.3的为0
test['prediction'] = (test['proba']>0.3).astype('int')
test['prediction']

步骤二：混淆矩阵

pd.crosstab(test.bad_ind,test.prediction,margins=True)

步骤三：根据混淆矩阵计算准确度

Int[]:acc = sum(test['prediction'] == test['bad_ind']) /np.float(len(test))
	  print('The accurancy is %.2f' %acc)

Out[]:The accurancy is 0.78

• 计算

- 精确度：A/A+B
- 召回率：A/A+C
- 特异度：D/C+D
for i in np.arange(0.02, 0.3, 0.02):              # 阈值0.02 - 0.3  步长0.02
    prediction = (test['proba'] > i).astype('int')
    confusion_matrix = pd.crosstab(test.bad_ind,prediction,
                                   margins = True)
    if i == 0.02:
        precision = 0
        recall = confusion_matrix.iloc[0,0] / confusion_matrix.iloc[2, 0]
        Specificity = 0
        f1_score = 2 * (precision * recall) / (precision + recall)
    else: 
        precision = confusion_matrix.iloc[0,0] / confusion_matrix.iloc[0, 2]
        recall = confusion_matrix.iloc[0,0] / confusion_matrix.iloc[2, 0]
        Specificity = confusion_matrix.iloc[1,1] / confusion_matrix.iloc[1,2]
        f1_score = 2 * (precision * recall) / (precision + recall)
    print('threshold: %s, precision: %.2f, recall:%.2f ,Specificity:%.2f , f1_score:%.2f'%(i, precision, recall, Specificity,f1_score))