数据科学实例教程之Logistic

前言

对于因变量为分类变量的分析常常使用逻辑回归模型。逻辑回归模型历史悠久，运算速度快，模型可以输出连续的概率预测值用于排序，常常用于信用评级等领域。由于计算高效，逻辑回归也常与其他模型融合，提高分类准确率。

本文整体是一个客户初始信用评级的案例，案例中的企业从事个人汽车金融服务，向购车的个人提供信用贷款。该公司的风控部门根据贷款申请者的基本属性，信贷历史，历史信用情况，贷款标的物的情况等信息构建贷款违约预测模型，根据模拟的预测做出相应的调整，从而来避免因大面积用户违约而造成公司大量的损失。

数据

本章使用汽车违约贷款数据集accepts.csv进行代码演示。

后文会用到显著性水平的判断，一般显著性水平定为0.05,即如果值小于0.05说明两者水平显著，反之即没有显著性。

下面是词条对显著性水平的解释

部分原始数据：

数据描述：

使用statsmodel分析包。案例所使用的包为

首先使用pandas读取数据，并进行简单清洗

f=open('./accepts.csv')

accepts=pd.read_csv(f,skipinitialspace=True)

accepts=accepts.dropna(axis=,how='any')

其中是否违约bad_ind是因变量

我们考虑曾经破产标识与是否违约之间是否有相关关系

cross_table=pd.crosstab(accepts.bankruptcy_ind,accepts.bad_ind,margins=True)

输出的交叉表如下

以曾经破产的样本为例，总计有310个，其中违约用户为67个，不违约的有243个。为了深入分析，将该表转换为列联表。

#转化为列联表

defpercConvert(ser):

returnser/float(ser[-1])

cross_table.apply(percConvert,axis=1)

输出结果如下

可以看到曾经破产的样本中，违约率为21.6%，而不破产的样本中违约率为18.9%。虽然在不同破产状况下违约率有差异，但这种差异是否显著，还需要使用卡方检验来判断。

#使用卡方检验

print('''chisq=%6.4f

p-value=%6.4f

dof=%i

expected_freq=%s'''%stats.chi2_contingency(cross_table.iloc[:2,:2]))

两变量卡方检验如下：

可以看到p值为0.28，说明差异可能并不显著，也即曾经是否破产与用户是否违约没有显著相关性。

接下来我们将数据集随机划分为两个部分，训练集和测试集，其中训练集用于模型的训练，测试集用于检验模型的泛化能力。

两数据集的样本量如下

训练集样本量: 2874

测试集样本量: 1231

经过抽样，训练集样本与测试集样本大致比例为7：3，这里因为是简单随机抽样，因此训练集和测试集当中的违约比例不一定是一样的，如果要保证训练集和测试集中的违约率相等，可以使用分层抽样，即在正例（所有违约的样本）和负例（所有未违约的样本）中分别抽取固定比例的样本，再合并成训练集和测试集。感兴趣的读者可以自行尝试。

我们使用训练集建立模型并且使用广义线性回归和logit变换对数据进行处理。

使用多元逻辑回归建模，并且使用summary查看模型的一些信息

在所有的参数当中，可以看到rev_util和veh_mileage的p系数不显著，因此可以删除，我们可以用变量筛选法来对变量进行筛选。例如使用AIC准则进行向前法变量筛选。

#向前法筛选变量

def forward_select(data, response):

remaining =set(data.columns)

remaining.remove(response)

selected = []

current_score, best_new_score =float(2397.2),float('inf')

whileremaining:

aic_with_candidates=[]

forcandidateinremaining:

formula ="{} ~ {}".format(

response,' + '.join(selected + [candidate]))

aic = smf.glm(

formula=formula,data=data,

family=sm.families.Binomial(sm.families.links.logit)

).fit().aic

aic_with_candidates.append((aic, candidate))

aic_with_candidates.sort(reverse=True)

best_new_score, best_candidate=aic_with_candidates.pop()

ifcurrent_score

remaining.remove(best_candidate)

selected.append(best_candidate)

print ('aic is {},continuing!'.format(best_new_score))

else:

print ('forward selection over!')

break

formula ="{} ~ {} ".format(response,' + '.join(selected))

print('final formula is {}'.format(formula))

model= smf.glm(

formula=formula,data=data,

family=sm.families.Binomial(sm.families.links.logit)

).fit()

return(model)

#只有连续变量可以进行变量筛选，分类变量需要进行WOE转换才可以进行变量筛选

candidates = ['bad_ind','fico_score','bankruptcy_ind','tot_derog','age_oldest_tr'

,'rev_util','ltv','veh_mileage']

data_for_select = train[candidates]

lg_m1 = forward_select(data=data_for_select, response='bad_ind')

lg_m1.summary().tables[1]

结果如下所示

可以看到不显著的变量被自动删除了，但是与线性回归相似，自变量的多重共线性会导致逻辑回归模型的不稳定，判断多重共线性可以使用方差膨胀因子，由于statsmodel中定义的方差膨胀因子计算函数的判别阀值与我们推荐的不一样(在此我们的判断阀值定位10)，所以我们自己定义一个方差膨胀因子的计算函数。如下所示。

defvif(df, col_i):

fromstatsmodels.formula.apiimportols

cols = list(df.columns)

cols.remove(col_i)

cols_noti = cols

formula = col_i +'~'+'+'.join(cols_noti)

r2 = ols(formula, df).fit().rsquared

return1./ (1.- r2)

candidates = ['bad_ind','fico_score','ltv','age_oldest_tr','tot_derog','nth','tot_open_tr','veh_mileage','rev_util']

exog = train[candidates].drop(['bad_ind'], axis=1)

foriinexog.columns:

print(i,'\t', vif(df=exog, col_i=i))

结果如下：

fico_score1.5423133089544319

tot_derog1.347832436613074

age_oldest_tr1.1399926313381807

rev_util1.084380320084259

ltv1.024624792276886

veh_mileage1.0105135995489773

可以看到VIF小于10这个阀值，说明自变量没有显著的多重共线性。我们写出fico_score的回归方程。

即，fico_score（信用评分，越高说明越好）每增加一个单位后的违约发生比是原违约发生比的0.987倍，也就是说明fico_score每增加一个单位后违约发生的可能性是原来的0.987倍，这与实际业务和常识是一致的。至于其它的变量分析，读者可以自行尝试，思路也与上面一样。

最后我们可以使用predict将违约的概率输出

#输出违约概率得分

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

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

test['proba'].head(10)

结果如下：

predict输出的是0~1之间的违约得分(也可以理解为违约概率)，即每个账户的违约概率我们已经得出，我们可以自己设置阀值来得到违约的预测标签，比如设置得分大于0.5的为违约，这个时候我们就可以重点‘关注’那些超过阀值的账户了。

# 设定阈值

test['prediction']=(test['proba']>0.5).astype('int')

这个时候我们的模型已经构建完毕，但是我们还需要对模型进行评估，我们计算模型的准确率Accuracy如下：

acc = sum(test['prediction'] ==test['bad_ind']) /np.float(len(test))

print('The accurancy is %.2f'%acc)

预测模型准确率如下：

The accurancyis0.82

可以看到，被正确预测的样本（包括正例与负例）占所有样本总数的82%，说明模型的效果还不错。但是要注意到，正例和负例的重要性是不同的，我们实际上需要更多地抓取正例(违约)，因为对于金融机构来说，违约造成的损失是远大于不违约带来的收益。我们使用ROC曲线来检测此次的模型，ROC曲线绘制方法如下：

#roc曲线绘制

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')

输出图如下：

我们使用sklearn.metrics模块来用于模型的评估，并且使用它自动生成不同阀值下的模型灵敏度、特异度，并绘制成曲线图。可以看到，训练集和测试集的预测结果比较接近（实线和虚线很接近），说明模型过拟合的可能很小。我们再来通过AUC来定量比较下模型，AUC的取值在0.5到1之间，数值越接近1说明模型效果越好。

#AUC定量比较

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

预测模型的ROC曲线下面积

AUC=0.7318

说明模型在测试集上的ROC曲线下的面积为0.7318。

至此我们的贷款违约预测模型算是大功告成，总体的效果还算不错。不过模型还是有许多可以优化的方面，感兴趣的读者可以自行尝试

If it works for you.Please,star.

自助者，天助之