（六）线性回归、逻辑回归

统计学习=机器学习，即根据指标找到最优超参数，确定模型。
统计学≠统计学习
线性回归、 逻辑回归是统计学和机器学习共用的两个算法。

1 线性回归

相关：检验X、Y是否独立
回归：构建Y与X的关系（线性关系、抛物线关系、对数关系） y=f(x)

1.1 简单线性回归

方法：最小似然估计

import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt          #可视化
import statsmodels.api as sm       
from statsmodels.formula.api import ols  #线性回归

• 导入数据

raw=pd.read_csv('creditcard_exp.csv',skipinitialspace=True)

• 对数据进行清洗

exp = raw[raw['avg_exp'].notnull()].copy().iloc[:, 2:].drop('age2',axis=1) #去除平均支出为0的客户，去除前两列和[age2]
exp

• 相关性分析

exp.describe(include='all')
exp.plot('Income', 'avg_exp', kind='scatter')  #两个连续变量->散点图 ‘scatter’
plt.show()

• 线性回归算法

lm_s = ols('avg_exp ~ Income', data=exp).fit() #拟合平均支出与收入的线性关系
lm_s.summary()

#p值 => β=0的概率，p=0代表β1=0的概率为0，非常显著

• 产生预测值
  在原始数据集上得到预测值和残差
  predict() 预测
  resid 残差=实际值-预测值

pd.DataFrame([lm_s.predict(exp), lm_s.resid], index=['predict', 'resid']).T.head()

在待预测数据集上得到预测值 ->对未开卡人做预测，如果未来开卡的消费支出值

exp_new= raw[raw['avg_exp'].notnull()].copy().iloc[:, 2:].drop('age2',axis=1)
lm_s.predict(exp_new).head()

1.2 多元线性回归

• 定义拟合模型

lm_m = ols('avg_exp ~ Age + Income + dist_home_val + dist_avg_income',
           data=exp).fit()     
lm_m.summary()

• 变量筛选
  向前法、向后法、逐步法
• 多元线性回归假设：使得模型满足某些假设

对x和y之间的关系的假设：
1.Y的平均值能够准确地被由X组成的线性函数建模出来 （根据理论、看散点图）
2. *解释变量X和扰动项不存在线性关系 （根据理论或常识判断，无法检验）
（如果存在线性关系，应将与X相关的尽量多的变量纳入模型中）
3. 解释变量之间不存在线性关系 （膨胀系数）
对y、€的假设：
4. 随机误差项€是一个均值为0的正态分布 （QQ检验）
5. 随机误差项方差恒为σ^2
6. 误差是独立的 （异方差检验、DW检验）

• 线性回归诊断：检验模型与假设是否一致，不一致调整至一致

1. 残差分析

首先画图，用散点图和残差图检查假设
观察，
残差中是否有离群值；
残差散点图是否和某个解释变量X有曲线关系（抛物线）
残差离散程度是否和某个解释变量有关
残差图：

非线性关系：加入X的高阶形式，一般加X2已经足够了
异方差：横截面数据，对Y取自然对数，弹性关系。
自相关：时间序列数据，周期性波动，加入Y的一阶滞后项进行回归

ana1=lm_s

exp['pred']=anal.predict(exp)   #查看残差图
exp['resid']=anal.resid
exp.plot('pred','resid',kind='scatter')  #p=βr   ‘income’和‘pred’结果一样
plt.show()

#散点图呈喇叭状，异方差，下一步优化应取对数

ana1 = ols('avg_exp ~ Income', exp).fit()   #拟合模型
ana1.summary()

#R^2为0.454

将lny放入模型

ana2 = ols('avg_exp_ln ~ Income', exp).fit()
exp['Pred'] = ana2.predict(exp)
exp['resid'] = ana2.resid
exp.plot('Income', 'resid',kind='scatter')

#散点图没有异方差，R^2为0.403，但是此时模型不符合常理，再对X取对数

exp['Income_ln'] = np.log(exp['Income']  #对X取对数

ana3 = ols('avg_exp_ln ~ Income_ln', exp).fit()
exp['Pred'] = ana3.predict(exp)
exp['resid'] = ana3.resid
exp.plot('Income_ln', 'resid',kind='scatter')
plt.show()

#在模型当中，取对数的变量是百分比的变化，即年收入每增加一个百分点，支出会增加0.9个百分点 ，R²为0.480，R²提高，得到目的。
模型解释的越符合常理，R^2越高

2. 强影响点分析

强影响点的统计量：学生化残差 = 残差 / 标准误差
强影响点的识别：
样本量为几百，|SR| > 2 为相对大的影响点
样本量为几千，|SR| > 3 为相对大的影响点
强影响点对于R²的影响不定，但去除强影响点，R²一般都会提高。异常值一般只去一次

exp['resid_t'] = (exp['resid'] - exp['resid'].mean()) / exp['resid'].std()      #学生化残差
exp[abs(exp['resid_t']) > 2]  #样本量为几百，认为大于2为强影响点
exp2 = exp[abs(exp['resid_t']) <= 2].copy()  #去除异常值
ana4 = ols('avg_exp_ln ~ Income_ln', exp2).fit()
exp2['Pred'] = ana4.predict(exp2)
exp2['resid'] = ana4.resid
exp2.plot('Income', 'resid', kind='scatter')
plt.show()
ana4.rsquared  #查看R^2值   0.49397191385172456

#散点图已没有信息，R²提高。

3. 多重共线性分析
   方差膨胀因子：Vif指标
   
   方差膨胀因子>10表示某变量的多重共线性严重

#增加变量
exp2['dist_home_val_ln'] = np.log(exp2['dist_home_val'])
exp2['dist_avg_income_ln'] = np.log(exp2['dist_avg_income'])

ana5 = ols('''avg_exp_ln ~ Age + Income_ln + 
           dist_home_val_ln + dist_avg_income_ln''', exp2).fit()
ana5.summary()
#多重共线性分析
def vif(df, col_i):          #定义vif
    cols = list(df.columns)
    cols.remove(col_i)
    cols_noti = cols
    formula = col_i + '~' + '+'.join(cols_noti)
    r2 = ols(formula, df).fit().rsquared
    
    return 1. / (1. - r2)

exog = exp2[['Income_ln', 'dist_home_val_ln','dist_avg_income_ln']]
for i in exog.columns:
    print(i, '\t', vif(df=exog, col_i=i))  

#将Vif大于10的剔除，剔除不显著的，删除即可。

建立线性模型流程：

初始分析：X、Y的分布，正态/右偏
->
变量选择：X、Y之间的关系，保留显著的关系
->
验证模型假定：散点图、经验假定X、Y之间的具体函数关系
->
多重共线性与强影响点诊断与分析
->
满足要求，得到模型

2 逻辑回归

Odds Ratio（优势比）：一个事件在一个组发生的可能性相对于在另一个组的大小
Odds（优势）：事件发生的概率与不发生概率之比，值域[0,∞)

#Odds Ratio将概率[0,1]范围扩大至∞

2.1 logit转换

一元连续变量逻辑回归，X连续变量 Y二分类变量
logit代表隐变量Y’，用X预测Y’



Logit是Odds的自然对数，logit(Pi)=事件发生概率的logit

i表示所有案例，观测值
Pi在第i个案例中一个事件发生的概率

2.2 模型评估

1. 决策类模型评估——混淆矩阵

设置阈值—> 真实值和预测值做交叉表，即混淆矩阵

正确率=(A+D)/(A+B+C+D)，只有正负样本1：1时才有意义
灵敏度Sensitivity（覆盖率）=A/(A+B)，重要
命中率Precision=A/(A+C)，预测违约的人中真的违约的人数
特异度Specificity（负的覆盖率）=D/(C+D)
负命中率=D/(D+B)

2. ROC曲线 （SAS的输出结果）
   模型预测的三种结果：一致对（√）、不一致对、相等对（50%√）

一致对表：

statmodels不提供一致对表，可用ROC曲线下面积（AUC）代表，AUC就是C统计量
ROC曲线y轴为灵敏度，x轴为1-特异度，曲线越趋向左上方越好

• 随机抽样，建立训练集与测试集

train = accepts.sample(frac=0.7, random_state=1234).copy()  #抽样70%，初始值1234随意设置
test = accepts[~ accepts.index.isin(train.index)].copy()  #找出没在的train内的样本 
print(' 训练集样本量: %i \n 测试集样本量: %i' %(len(train), len(test)))

• 建立逻辑回归模型

lg = smf.glm('bad_ind ~ age_oldest_tr', data=train,           family=sm.families.Binomial(sm.families.links.logit)).fit()  #Binomial伯努利分布，0-1分布
lg.summary()

• 预测

train['proba'] = lg.predict(train)
test['proba'] = lg.predict(test)

• 模型评估

#设置阈值
test['prediction'] = (test['proba'] > 0.3).astype('int')
#混淆矩阵
pd.crosstab(test.bad_ind, test.prediction, margins=True)
#计算准确度  ？？计算混淆矩阵指标
acc = sum(test['prediction'] == test['bad_ind']) /np.float(len(test))
print('The accurancy is %.2f' %acc)  #The accurancy is 0.77
#ROC曲线
import sklearn.metrics as metrics    

fpr_test, tpr_test, th_test = metrics.roc_curve(test.bad_ind, test.proba)
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()

print('AUC = %.4f' %metrics.auc(fpr_test, tpr_test))   #AUC = 0.6136