基于Toad及逻辑回归算法的风控评分模型

说明：文章内容记录核心代码及主要流程，核心标签及数据分析结果只做示例展示，核心标签以中文展示。

Python数据分析及风控评分模型包导入

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import missingno as msno
import seaborn as sns 
import math 
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
from sklearn import metrics 
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.tree import DecisionTreeClassifier,_tree
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.filterwarnings('ignore')
from sklearn.tree import _tree
from sklearn.metrics import roc_curve

pip install toad # 安装Toad
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple toad
# 中文指引：https://toad.readthedocs.io/en/latest/tutorial_chinese.html

读取数据及探索性分析

读取标签

zxbq = pd.read_excel(r"E:\data\数据字典.xlsx",sheet_name = "数据字典")
zxbq[['字段名称','中文含义']]
zxbq.head()
zxbq[['字段名称', '中文含义']]
df_zxbq = zxbq[['字段名称', '中文含义']]
df_zxbq.loc[df_zxbq['字段名称'] == 'cnt_qry_nbank_ld720'] # 查询标签
print(zxbq[['字段名称','中文含义']]['字段名称'].value_counts)
zxbq[['字段名称','中文含义']]['字段名称'].values

读取数据

data = pd.read_excel(r"E:\data.xlsx",sheet_name = "Sheet")
data_copy = data
data_copy.head()
data_copy.info()
data_copy = data_copy.drop_duplicates() # 删除重复数据

查看数据集缺失值

# 缺失值分布
import missingno as msn
msn.matrix(data_copy, labels = False, label_rotation = 90)
msno.matrix(data_copy)

# 缺失值统计
count_missing = data.apply(lambda x:'{}%'.format(round(100*sum(x.isnull())/len(x),2)))
print(count_missing) # 统计缺失率

# 每个变量缺失率的计算
def missing_cal(df):
    """
    df :数据集
    
    return：每个变量的缺失率
    """
    missing_series = df.isnull().sum()/df.shape[0]
    missing_df = pd.DataFrame(missing_series).reset_index()
    missing_df = missing_df.rename(columns={'index':'col',
                                            0:'missing_pct'})
    missing_df = missing_df.sort_values('missing_pct',ascending=False).reset_index(drop=True)
    return missing_df
missing_cal(data_copy)

# 变量的缺失分布图
def plot_missing_var(df,plt_size=None):
    """
    df: 数据集
    plt_size :图纸的尺寸
    
    return: 缺失分布图（直方图形式)
    """
    missing_df = missing_cal(df)
    plt.figure(figsize=plt_size)
    plt.rcParams['font.sans-serif']=['Microsoft YaHei']
    plt.rcParams['axes.unicode_minus'] = False
    x = missing_df['missing_pct']
    plt.hist(x=x,bins=np.arange(0,1.1,0.1),color='hotpink',ec='k',alpha=0.8)
    plt.ylabel('缺失值个数')
    plt.xlabel('缺失率')
    return plt.show()
plot_missing_var(data_copy)

# 单个样本的缺失分布
def plot_missing_user(df,plt_size=None):
    """
    df: 数据集
    plt_size: 图纸的尺寸
    
    return :缺失分布图（折线图形式）
    """
    missing_series = df.isnull().sum(axis=1)
    list_missing_num  = sorted(list(missing_series.values))
    plt.figure(figsize=plt_size)
    plt.rcParams['font.sans-serif']=['Microsoft YaHei']
    plt.rcParams['axes.unicode_minus'] = False
    plt.plot(range(df.shape[0]),list_missing_num)
    plt.ylabel('缺失变量个数')
    plt.xlabel('sanples')
    return plt.show()
plot_missing_user(data_copy)

# 缺失值填充（类别型变量）
def fillna_cate_var(df,col_list,fill_type=None):
    """
    df:数据集
    col_list:变量list集合
    fill_type: 填充方式：众数/当做一个类别
    
    return :填充后的数据集
    """
    df2 = df.copy()
    for col in col_list:
        if fill_type=='class':
            df2[col] = df2[col].fillna('unknown')
        if fill_type=='mode':
            df2[col] = df2[col].fillna(df2[col].mode()[0])
    return df2

data_copy['als_m6_id_nbank_orgnum'].value_counts()
fillna_cate_var(data_copy, 'als_m6_id_nbank_orgnum', fill_type = None)

填充数据集缺失值

data_copy['flag'] = data_copy['flag'].fillna(-999)
data_copy['flag']
data_copy['逾期天数'] = data_copy['逾期天数'].fillna(-999)
data_copy['逾期天数'] # 填充series缺失值

数据探索性分析

import toad as td
data_copy_eda = td.detect(data_copy)
data_copy_eda
data_copy.describe()

数据集划分

data_train = data_copy.loc[data_copy['数据时间']<='2022-09'].reset_index(drop = True)
print(data_train.shape)
data_test = data_copy.loc[data_copy['数据时间']>'2022-09'].reset_index(drop = True)
print(data_test.shape)

数据初步筛选及评估

data_select = td.selection.select(data_train, target = 'flag',
                           empty = 0.9, #剔除缺失率高于阈值的变量
                           iv = 0.02, #剔除IV低于阈值的变量
                           corr = 0.2, #剔除相关性高于阈值且IV较低的变量
                           return_drop = True, #是否返回剔除的变量名称，当为TRUE时返回数据集和剔除的变量字典 
                           exclude = ['数据时间']) #不考虑的变量,默认保留
data_select[0].head() # 数据集
data_select[1] # 通过td.selection.select删除的变量字典

to_drop = [ '数据时间']
td.quality(data_select[0].drop(to_drop, axis = 1), 'flag', iv_only = True)
# 数据量较大或高维度数据，需要使用iv_only = True
# 删除数据主键，评分、日期且不用于建模的特征

Combiner 单变量分箱

df_zxbq.loc[df_zxbq['字段名称'] == '近720天查询机构数'] # 查询标签

#Combiner分箱功能
c = td.transform.Combiner()
c.fit(data_select[0], y = 'flag',
      method = 'chi', #分箱方法，支持’chi’ (卡方分箱), ‘dt’ (决策树分箱), ‘kmean’ , ‘quantile’ (等频分箱), ‘step’ (等步长分箱)
      min_samples = None, #每箱至少包含样本量，可以是数字或者占比
      n_bins = 3, #箱数，若无法分出这么多箱数，则会分出最多的箱数
      empty_separate = False) #是否将空箱单独分开

c.export()['近720天查询机构数'] #查看分箱节点
rule = {'近720天查询机构数' : [0.0, 1669.0]}
c.update(rule) #手动更新分箱节点,load\update\set_rules

from toad.plot import bin_plot, badrate_plot # 查看变量单调性
bin_plot(c.transform(data_select[0][['近720天查询机构数', 'flag']], labels = True), 
         x = '近720天查询机构数', 
         target = 'flag',
         iv = True) # 不考虑时间因素，单纯比较变量在各个分箱中的单调性

badrate_plot(c.transform(data_select[0][['近720天查询机构数', 'flag', '数据时间']], 
             labels = True), 
             x = '数据时间', 
             target = 'flag',
             by = '近720天查询机构数')

WOE

transer = td.transform.WOETransformer() # 初始化woe transer

data_woe = transer.fit_transform(c.transform(data_select[0], labels = True),
                                 data_select[0]['flag'],
                                 exclude = 'flag') #注意删除特征

逐步回归

data_final = td.selection.stepwise(data_woe, 
             target = 'flag', 
             estimator = 'ols', #用于拟合的模型，支持'ols', 'lr', 'lasso', 'ridge'
             direction = 'both', #逐步回归的方向，支持'forward', 'backward', 'both' （推荐） 
             criterion = 'aic',  #评判标准，支持'aic', 'bic', 'ks', 'auc'
             max_iter = None,  #最大循环次数
             return_drop = False,  #是否返回被剔除的列名
             exclude = None) #不需要被训练的列名

逻辑回归建模

import statsmodels.api as sm
data_final = sm.add_constant(data_final) #增加截距项，默认无截距
model_lr = sm.Logit(data_final['flag'], data_final.iloc[:, :-1]).fit()
model_lr.summary()

from statsmodels.stats.outliers_influence import variance_inflation_factor
def Check_VIF(df):
    VIF_list = [variance_inflation_factor(np.matrix(data_final), i) for i in range(df.shape[1])]
    col_name = list(df.columns)
    VIF_data = pd.DataFrame({'col_name' : col_name, 'VIF_value' : VIF_list})
    return VIF_data
Check_VIF(data_final)

模型评估及模型评价

from toad.metrics import KS, AUC, PSI, roc_curve
tmp = PSI(test = data_final.iloc[1000:, :], 
          base = data_final.iloc[:1000, :], 
          return_frame = True)
tmp

from toad.metrics import KS, AUC, PSI, roc_curve, KS_bucket
data_final['y_pre'] = model_lr.predict(data_final.iloc[:, :-1])
auc = AUC(data_final['y_pre'], data_final['flag'])
ks = KS(data_final['y_pre'], data_final['flag'])

ks_data = KS_bucket(data_final['y_pre'], data_final['flag'], 
                    bucket = 10, #箱体数量
                    method = 'quantile') #分箱方法，可以选择等频quantile或等长step
ks_data

模型评估结果展示

def perf_eva(score, target, data_name, if_positive = True):
    
    if if_positive == False:
        score = score * -1
    
    fpr_dev, tpr_dev, _  = roc_curve(target, score)
    auc = AUC(score, target)
    
    fig, axes = plt.subplots(1, 2, figsize = [12, 6])
    ax1 = axes[0] 
    ax2 = axes[1] 
    
    xlim = np.linspace(0, 1, len(fpr_dev))
    ks_dev = tpr_dev - fpr_dev
    ks = ks_dev.max()
    ks_argmax = ks_dev.argmax()
    ax1.plot(xlim, fpr_dev)
    ax1.plot(xlim, tpr_dev)
    ax1.plot(xlim, ks_dev)
    ax1.plot([xlim[ks_argmax]] * 20, np.linspace(0, ks, 20), '--')
    ax1.set_xlabel('% of population')
    ax1.set_ylabel('% of total Good / Bad')
    ax1.set_title('K-S', fontsize = 14, fontweight = 'semibold')
    ax1.text(xlim[ks_argmax] - 0.15, ks, 'KS : {0}'.format(np.round(ks, 4)), fontsize = 14)
    
    ax2.plot(fpr_dev, tpr_dev)
    ax2.plot([0, 1], [0, 1], '--')
    ax2.set_xlabel('FPR')
    ax2.set_ylabel('TPR')
    ax2.set_title('ROC Curve',  fontsize = 14, fontweight = 'semibold')
    ax2.text(0.35, 0.5, 'AUC : {0}'.format(np.round(auc, 4)), fontsize = 14)
    
    fig.suptitle('{0} : Model Performance'.format(data_name), fontsize = 14, fontweight = 'black')
    plt.show()

    return fig

perf_eva(score = data_final['y_pre'], target = data_final['flag'],
         data_name = 'Train')

标准评分卡转化

pdo = 60 #步长，odds每翻rate倍，增加的分数
rate = 2 
base_odds = 20 #good / bad
base_score = 1000

B = pdo / np.log(rate)
A = base_score - B * np.log(base_odds)
# score = A - B * np.log(p / (1 - p))

data_final['score'] = data_final['y_pre'].apply(lambda x: A - B * np.log(x / (1 - x))) 
data_final['score'].hist(bins = 10)
perf_eva(score = data_final['score'], target = data_final['flag'], data_name = 'Train', if_positive = False)

lr_params = model_lr.params
col_name = []
col_value = []
col_score = []
base_score = base_score - B * lr_params['const'] 
col_name.append('const')
col_value.append(1)
col_score.append(base_score)
col_woes = transer.export()
for col in list(lr_params.index)[1: ]:
    col_param = lr_params[col]
    dic = col_woes[col]
    for key in list(dic.keys()):
        col_name.append(col)
        col_value.append(key)
        col_score.append(dic[key] * col_param * - B)
score_card = pd.DataFrame({'col_name' : col_name,
                           'col_value' : col_value,
                           'col_scoer' : col_score})

2023年12月27日
青海西宁