基于机器学习的股市投资避雷指南分析实战
玄武约 46 分钟
本实战案共分为六个部分：商业背景、指标设计、部署环境、数据准备、探索分析，回归分析，其中回归分析包括：筛选解释变量，子集法筛选，拟合回归模型，评估模型性能，探测交互效应，评估交互模型，模型结果呈现。 本文着实有点长，多少需要忍下，建议 点赞收藏 ！
©️ 数据STUDIO投稿 · 作者｜ 玄武
1.商业背景
我们知道，特别处理是我国股市特有的一项旨在保护投资者利益的政策，它第一次出现在人们视野是在1998年4月22日，当时，沪深交易所宣布，将对财务状况或其它状况出现异常的上市公司股票交易进行特别处理，在其股票简称前冠以“ST”标记，这类股票称为ST股，这项制度叫做ST制度。
根据ST退市最新规则，在财务方面：
-
第一，年报扣非前或后净利润取低者为负且营收低于1亿元，将被戴上*ST，如果已经是*ST的，直接退市； -
第二，2020年起第一年审计资产为负，将被载上*ST，连续两年触及该指标，将被终止上市。
于是，公司被ST就意味着投资风险增大，为了降低投资风险，本案例研究的目的是寻找出上市公司被特别处理的影响因素，并预测上市公司被特别处理的可能性。本研究根据相关文献，围绕着公司的财务指标，总共设计了7个解释变量。
2.指标设计
第一个指标：应收账款占比。
该指标是应收账款占总资产的比例，当然你也可以用应收账款占营业收入的比例，应收账款占比反映了盈利质量。
应收账款表示应该收回来但是还没有收回来的货款，存在以下问题：第一，流动性差，显然应收账款无法满足一般性的支付需求；第二，高风险，如果到了待付款期限，对方无力付款，这会给公司造成巨大的损失。
因此，从企业运营角度看，应收账款应越少越好。如果应收账款占比低，这说明该企业的收入质量高；相反，如果应收账款高，这说明该企业的收入质量低。
第二个指标：对数资产规模。
该指标是对数变换后的资产规模，用于反映公司规模。这不是反映公司规模的唯一指标，甚至也不是最好的指标。例如，还可以考虑净资产规模，当然，净资产规模也有它的问题，因为对于资不抵债的企业而言，净资产是负的。
对于金融类企来说，如基金公司，它们所管理的资产规模也许是一个不错的指标。但是，此类指标对于制造类企业也许不适用。一个资金规模不大的制造业企业，也许就有成百上千的员工。对于此类企业，也许员工数量是一个更好的选择。对于一个互联网企业，这两种方法可能都不好，这类企业往往资产很轻，因此资产规模很小，它们往往强调技术，因此人员规模也不大，但是，有可能它们的产品有过亿用户，因此，对这类企业，也许用户规模是一个更合理的选择。
第三个指标：资产周转率。
该指标是衡量企业资产管理效率的重要财务比率，是一个企业在一定时期内的销售收入净额除以资产平均总额而得。
假设两个不同的企业A和B都有一单位的平均资产。在某一年内，A企业总共做了10单生产，B企业只做了1单。因此，A企业的营业额会高于B企业，A企业的资产周转率也会高于B企业。这说明A企业对资产的利用率要高于B企业。
因此，资产周转率在财务分析指标体系中具有重要地位。通过对该指标的对比分析，能够在一定程度上了解企业总资产的运营效率。这个指标数值越高，表明企业总资产周转速度越快，资产利用效率就越高，经营状况越好，被ST的风险应该越低。
第四个指标：资产回报率。
该指标是企业在一定时期内的利润总额除以总资产，它反映的是每单位资产能够给企业带来的利润。因此，该指标可以看作对企业盈利能力的反映。但它不是反映企业盈利的唯一指标，甚至不是最好的指标。例如，在净资产不是负数时还可以考虑净资产收益率。
对于某些特定行业，人们还会关注销售收益率等。对于很多处在快速扩张的互联网企业，这也不是一个好的指标，因为它们常常是亏损的，但是，资本市场对这类企业仍然非常认可，一个关键原因就是看好它们未来的长期价值。
第五个指标：销售收入增长率。
该指标是企业在一定时期内的销售总额除以前一个时期的销售总额，它反映的是企业的增长速度。
对于高科技企业，在其成立之初很难实现盈利，但并不妨碍企业高速增长。企业的高速增长会反映在什么指标？是销售收入？还是资产或净资产？还是市场占有率等其它非财务指标？企业的高速成长会如何影响其盈利和被ST的概率？这不是一个简单问题。如果企业的销售是盈利的，那么高速成长的销售带来的应该是更好的盈利和更小的ST概率，反之，高速成长的销售带来的应该是更大的亏损和更大的ST概率。
值得注意的是，销售收入增长率也不是唯一最好的关于企业成长速度的指标，如新兴互联网企业，在疯狂扩张初期，不但没有利润连销售收入也很少，但他们的用户规模却在飞速增长。
第六个指标：杠杆水平。
该指标是债务资产比率，是企业债务在其总资产中所占比率，反映的是企业的总资产中来自债权人的比率。
企业的债务资产比率是如何影响企业盈利，进而影响其被ST的概率呢？这是一个颇有争议的问题。首先，过高的债务资产比显然不好，这会使企业背上沉重的债务负担，使得企业每年盈利中的一大部分将用于偿还利息，因此损伤盈利。但是，另外一方面，几乎没有企业不举债，适当举债能给企业带来很多好处，如很多创业初期的企业，发展势头很好，但是缺乏资金，那么合理举债能够帮助企业迅速成长，占领终端市场，确立优势。
第七个指标：第一大股东持股比率。
该指标反映的是企业的股权结构。如果企业的第一大股东持股比例很高，比如大于50%，说明该企业一股独大，其持有者对企业的方方面面具有绝对话语权；如果企业的第一大股东持股比例很低，如小于10%，说明该企业股权非常分散。企业的股权结构会如何影响盈利？过度分散的股权结构是不好的，因为这使得所有人都不会真正地关心企业，承担责任。过度集中的股权结构也不好，因为这使得第一大股东有能力侵害小股东的利益。大股东占款被认为是中国股市的一个重大顽疾，被广为诟病，会影响到上市企业的正常经营行为。但是，大股东占款会如何影响上市企业的ST风险，却并不为人熟知。
我们这个研究，就是想看看在这7个指标当中，哪些指标会影响上市公司被特别处理的概率，并对上市公司是否会被ST进行预测。
3.部署环境
3.1导入依赖库
import numpy as npimport pandas as pdimport statsmodels.api as smfrom scipy import statsfrom itertools import combinationsfrom sklearn.metrics import accuracy_score,confusion_matrix,classification_report,roc_curve,roc_auc_score # 预测准确率、混淆矩阵、分类报告、ROC曲线、AUC值from IPython.core.display_functions import display3.2 设置可视化
import seaborn as snssns.set_style('white') # 当出现pycharm风格错乱时，可使用此语句修正from matplotlib import pyplot as pltplt.rcParams['xtick.bottom']=True # 设置横轴刻度线可见plt.rcParams['ytick.left']=True # 设置纵轴刻度线可见plt.rcParams['figure.autolayout']=True #自动tight布局，解决多个plots出现重叠问题plt.rcParams['axes.unicode_minus']=False # 解决负号显示问题plt.rcParams['font.sans-serif']=['Microsoft YaHei'] # 解决中文显示问题plt.rcParams['figure.dpi']=80 # 设备图形分辨率%matplotlib inline3.3 忽略警告信息
详情可以参见这篇： Python 警告控制模块：warnings
import warningswarnings.filterwarnings('ignore')3.4 定义函数
上下滑动查看更多
绘制柱状图(百分比)
def bar_plot_percent(series_data,series_data_label): """ 绘制柱状图(百分比) """ # 确定柱子在画布纵坐标上的位置，即各类别的频数 y_axes=series_data.value_counts(normalize=True).sort_index() # 确定柱子在横坐标上的位置 x_axes=np.arange(len(y_axes))+1 # x代表柱子在图形中的位置，y_axes代表柱子在图形中的高度，facecolor代表柱子颜色 plt.bar(x_axes,y_axes) plt.xlabel(series_data_label) # 设置横轴的名称 plt.ylabel('百分比') # 设置纵轴的名称 plt.xlim(0,len(y_axes)+1) # 设置横坐标的区间 plt.ylim(0,y_axes.max()*1.1) # 设置纵坐标的区间 plt.xticks(x_axes,y_axes.index) # 设置柱子标签 # 在柱子顶端显示柱子高度的数据 for bar_num in range(len(y_axes)): # 设置柱子高度数据在横轴上的位置 plt.text(x_axes[bar_num], # 设置柱子高度数据在纵轴上的位置，必须丢弃索引， # 否则作为数字的bar_num将无法与原索引匹配 y_axes.reset_index(drop=True)[bar_num], # 设置柱子高度的具体显示数据，必须丢弃索引， # 否则作为数字的bar_num将无法与原索引匹配 '{:.1%}'.format(y_axes.reset_index(drop=True)[bar_num]), # 垂直对齐方式，参数：[‘center’|‘right’|‘left’] horizontalalignment='center', # 水平对齐方式，参数：[‘center’|‘top’|‘bottom’|‘baseline’] verticalalignment='bottom' ) # plt.gca().spines[['right','top']].set_visible(False) # 设置右框、上框不可见，gca='get current axis'绘制分组箱线图
def grouping_boxplot( grouping_variable, grouping_variable_label, box_variable, box_variable_label, data):# 绘制分组箱线图 boxplot_data = [data[data[grouping_variable] == value][box_variable].dropna() for value in data[grouping_variable].value_counts().index] # 构造用于制作箱形图的数据 plt.boxplot(boxplot_data, # 数据结构为列表套列表，一个列表绘制一个箱线图 labels=data[grouping_variable].value_counts().index, # 每个箱线图在横轴的标签 patch_artist=False, # 设置是否往箱体填充颜色，默认不填充(False) showmeans=True, # 在箱子内标出均值所在位置，箱体内部的横线为中位数 widths=0.5, # 设置箱子的宽度，默认为0.5 showfliers=False, # 设置是否显示异常值，默认显示异常值(True) medianprops={'linestyle': '-', 'color': 'black'} # 设置中位数线的类型和颜色 ) plt.xlabel(grouping_variable_label) plt.ylabel(box_variable_label)步进法(stepwise)选择进入模型的自变量
def stepwise_select_variable( x_var, y_var, pvalue_in=0.05, pvalue_out=0.1):"""步进法(stepwise)选择进入模型的自变量：如果p值小于pvalue_in，则往模型中放入此变量；如果p值大于pvalue_out，则从模型中移出此变量""" variables_in=[] # 用于存储模型中的变量 while True: # 直接设置为True，则让while永远循环，直到遇到break才跳出循环 flag=False # 挑选p值最小且小于pvalue_in的自变量 variables_out=list(set(x_var.columns)-set(variables_in)) # 用于存储未进入模型的变量，其值为矩阵x的所有变量列减去模型内的变量列（set是一个集合，无序且内容不能重复） pvalues=pd.Series(index=variables_out,dtype='float64') # 创建一个series（索引为未进入模型的变量的名称），用于存储自变量回归系数(不含截距)的p值 for new_column in variables_out: # 遍历未进入模型的变量 x_var_cons=sm.add_constant(x_var[variables_in+[new_column]]) # 确定自变量，并在最左边加上一列全为1的数据，使得模型矩阵中包含截距 model_stepwise=sm.Logit(y_var,x_var_cons).fit() # 用未标准化数据拟合模型：如自变量为x_var_cons，则拟合含截距模型；如自变量为x_var，则拟合不含截距模型 pvalues[new_column]=model_stepwise.pvalues[new_column] # 将刚进入模型的变量的回归系数的p值存储到pvalues列表（不含截距的p值） sig_pvalue=pvalues.min() # 获取pvalues列表中的最小值 if sig_pvalue<pvalue_in: flag=True # 由于存在小于pvalue_in的p值，所以循环得继续，不可跳出循环，因此，设置flag为True sig_variable=pvalues.index[pvalues.argmin()] # 获取最小p值对应的自变量（argmin返回最小值的坐标，如遇多个最小值，则返回第1个最小值的坐标） variables_in.append(sig_variable) # 将最小p值对应的自变量添加到用于储存模型中的变量的variables_in列表 # 制作系数表格（含变量名称、非标准化回归系数、t值、p值） x_var_cons=sm.add_constant(x_var[variables_in]) model_stepwise=sm.Logit(y_var,x_var_cons).fit() # 用未标准化数据拟合模型：如自变量为x_var_cons，则拟合含截距模型；如自变量为x_var，则拟合不含截距模型 print('variable in: {:10} Pseudo R-squ： {:5} 模型系数的Omnibus检验p值： {:5}'.format(sig_variable,model_stepwise.prsquared,model_stepwise.llr_pvalue)) coefficient=( pd.DataFrame({'未标准化回归系数':model_stepwise.params}) .assign(t值=model_stepwise.tvalues) # 增加“t值"列 .assign(显著性=model_stepwise.pvalues) # 增加”显著性“列 ) display(coefficient) # 剔除p值最大且大于pvalue_out的自变量 x_var_cons=sm.add_constant(x_var[variables_in]) # 确定自变量，并在最左边加上一列全为1的数据，使得模型矩阵中包含截距 model_stepwise=sm.Logit(y_var,x_var_cons).fit() # 用未标准化数据拟合模型：如自变量为x_var_cons，则拟合含截距模型；如自变量为x_var，则拟合不含截距模型 pvalues=model_stepwise.pvalues.iloc[1:] # 除截距p值外的所有自变量的p值 notsig_pvalue=pvalues.max() # 获取pvalues列表中的最大值 if notsig_pvalue>pvalue_out: flag=True # 由于存在大于pvalue_out的p值，所以循环得继续，不可跳出循环，因此，设置flag为True notsig_variable=pvalues.index[pvalues.argmax()] # 获取最大p值对应的自变量（argmax返回最大值的坐标，如遇多个最大值，则返回第1个最大值的坐标） variables_in.remove(notsig_variable) # 将最大p值对应的自变量移出用于储存模型中的变量的variables_in列表 # 制作系数表格（含变量名称、非标准化回归系数、t值、p值） x_var_cons=sm.add_constant(x_var[variables_in]) model_stepwise=sm.Logit(y_var,x_var_cons).fit() # 用未标准化数据拟合模型：如自变量为x_var_cons，则拟合含截距模型；如自变量为x_var，则拟合不含截距模型 print('variable in: {:10} Pseudo R-squ： {:5} 模型系数的Omnibus检验p值： {:5}'.format(notsig_variable,model_stepwise.prsquared,model_stepwise.llr_pvalue)) coefficient=( pd.DataFrame({'未标准化回归系数':model_stepwise.params}) .assign(t值=model_stepwise.tvalues) # 增加“t值"列 .assign(显著性=model_stepwise.pvalues) # 增加”显著性“列 ) display(coefficient) if not flag: # 当flag为False时（既没有需要进入模型的变量，也没有需要剔除的变量），跳出循环 break return variables_in# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
4.数据准备
4.1原始数据导入
df=( pd.read_csv('数据STUDIO/上市公司被特别处理(ST)的影响因素研究.csv') # 数据获取：公众号：数据STUDIO，后台回复【data】 # 将字段名称全部统一成小写，以免由于大小写问题而造成代码报错 .rename(columns=lambda z:z.lower()) # 0表示”非ST“，1表示”ST“ .replace(['非ST','ST'],[0,1]) )display(df)4.2数据记录去重
if df.duplicated().sum()>0: # 显示重复数据记录 display(df[df.duplicated()]) print('数据集中存在以上{}条重复数据记录，现已删除。'.format(df[df.duplicated()].shape[0])) # 删除所有字段都完全重复的数据并立即生效 df.drop_duplicates(inplace=True) display(df)else: print('数据集中不存在重复数据记录，无需去重。')数据集中不存在重复数据记录，无需去重。5.探索分析
5.1数据的描述分析
df.describe(percentiles=(0.01,0.25,0.5,0.75,0.99)).T结果解读：
每个变量都有684条数据记录，不存在缺失值，未发现明显的极端值，所有字段都不存在违背业务逻辑的数据记录。
5.2 数据的分布形态
查看数据分布情况，有助于特征工程根据数据分布选择合适的数据处理办法（包括缺失值、异常值处理，连续特征离散化），还有助于深入了解用户行为。
对于连续数据，当偏度系数等于0时，数据呈左右对称分布；当偏度系数绝对值大于等于1时，数据呈严重偏斜分布；当偏度系数绝对值大于等于0.5并且小于1时，数据呈中等偏斜分布；当偏度系数绝对值大于0并且小于0.5时，数据呈轻微偏斜分布。
对于分类数据，主要观察柱状图的左右对称性。
本数据集全为连续数据，没有分类数据。
# 偏度系数分析skw_analysis=( # 只计算数值型字段的偏度系数，两个日期型字段不计算 pd.DataFrame({'偏度系数':df.skew(numeric_only=True)}) .assign(偏斜程度=df.skew(numeric_only=True) \ .to_frame()[0] \ .apply( lambda z:'严重' if abs(z)>=1 else '中等' if abs(z)>=0.5 else '轻微') ) # 新增“偏斜程度”列 .assign(偏斜方向=df.skew(numeric_only=True) \ .to_frame()[0] \ .apply( lambda z:'左偏' if z<0 else '对称' if z==0 else '右偏') ) # 新增“偏斜方向”列) .sort_values(by='偏度系数', ascending=True) # 按“偏度系数”变量升序排序)display(skw_analysis)plt.figure(figsize=(16,8)) # 确定画布尺寸大小# 为各自变量绘制直方图for i in range(241,248): # 将画布分为2*4=8个子图区域（2表示分2行，4表示分4列）， # 每个图的编号依次为1、2、3、4、5、6、7、8，在这里只使用1-7号子图 plt.subplot(i) plt.hist(df[['ara', 'asset', 'ato', 'roa', 'growth', 'lev', 'share'][i-241]]) plt.xlabel(['应收账款占比(ara)','对数资产规模(asset)', '资产周转率(ato)','资产回报率(roa)', '销售收入增长率(growth)','杠杆水平(lev)', '第一大股东持股比率(share)'][i-241]) plt.ylabel('频数')# 为因变量绘制柱状图plt.subplot(248)bar_plot_percent(df.st,'是否ST')plt.xticks([1,2],['否','是']) 结果解读：
-
第一，应收账款占比(ara)的偏度系数为1.94，呈严重右偏分布，表示大多数样本的应收账款占比取值偏小，说明绝大多数企业的应收账款占比不高。但是，数据分布的严重右偏同时说明，有不少企业的应收账款占比很高，直方图显示有不少企业的应收账款占比超过50%，最大值高达63.5%。 -
第二，资产周转率(ato)、资产回报率(roa)的偏度系数分别为2.22、1.70，都呈严重右偏分布，后续建模时可考虑对资产周转率(ato)、资产回报率(roa)这两个变量数据取自然对数。 -
第三，对数资产规模(asset)的偏度系数为0.55，呈中等右偏分布，后续建模时可直接进入模型，无须进一步进行转变。 -
第四，销售收入增长率(growth)、杠杆水平(lev)、第一大股东持股比率(share)的偏度系数分别为0.01、0.08、0.09，呈轻微右偏分布，后续建模可直接进入模型，无须进一步进行转变。 -
第五，目标变量st取值为1的占比为5.3%，存在不平衡问题，这将导致模型预测较难将正样本预测出来(因为过度的照顾了负样本)，使得AUC值过低。对于样本分布不平衡的数据集，一般需先进行过采样或欠采样，然后再进行建模。
5.3数据的正态检验
检验数据是否服从正态分布的方法比较多，不同检验方法对样本量的敏感度不一样：样本量n<50时，优先使用Shapiro-Wilk检验；50≤样本量n<5000时，酌情使用W检验及K-S检验；5000≤样本量n时，建议使用K-S检验。
不过，我更推荐使用D’Agostino and Pearson omnibus normality test，因为这是一种通用和强大的正态性检验方法，其基本思想：首先，计算偏度和峰度以便在不对称和形状方面量化分布离正态分布的距离：然后，计算这些值中的每一个与正态分布的预期值之间的差异，并基于这些差异的总和，计算各P值。
"""不能对字符串型数据检验是否服从正态分布，axis=1表示删除列，axis=0表示删除行。"""stats.normaltest(df.drop('st',axis=1)) NormaltestResult( statistic=array( [281.34589097, 38.95314649, 348.78258169, 270.09432854, 18.36113533, 6.52235294, 116.62349139]), pvalue=array( [8.06335948e-62, 3.47881805e-09, 1.83157495e-76, 2.23749339e-59, 1.03022035e-04, 3.83432618e-02, 4.73729701e-26]))结果解读：
所有p值都小于0.05，都拒绝原假设，认为数据都不服从正态分布。
5.4数据的相关性分析
sns.pairplot(df.drop('st',axis=1)) df.corr(method='pearson',numeric_only=True)结果解读：
除杠杆水平(lev)与资产回报率(roa)的皮尔逊相关系数为-0.41外，其它所有相关系数的绝对值都小于0.3，说明自变量之间的相关性不强。
5.5数据的箱线图分析
"""绘制分组箱线图"""# 确定画布尺寸大小plt.figure(figsize=(16,8)) for i in range(241,248): """ 将画布分为2*4=8个子图区域， 每个图的编号依次为1、2、3、4、5、6、7、8， 在这里只使用1-7号子图 """ plt.subplot(i) grouping_boxplot( 'st', '是否ST', ['ara', 'asset', 'ato', 'roa', 'growth', 'lev', 'share'][i-241], ['应收账款占比(ara)','对数资产规模(asset)', '资产周转率(ato)','资产回报率(roa)', '销售收入增长率(growth)','杠杆水平(lev)', '第一大股东持股比率(share)'][i-241], data=df) plt.xticks([1,2],['否','是'])结果解读：
-
第一，“应收账款占比(ara)”分组箱线图显示，非ST组的中位数和均值都要明显低于ST组，说明应收账款占比高的公司更容易被ST。 -
第二，“对数资产规模(asset)”分组箱线图显示，非ST组的中位数和均值与ST组差异不明显T。 -
第三，“资产周转率(ato)”分组箱线图显示，非ST组的中位数和均值都要明显高于ST组，说明资产周转率低的公司更容易被ST。 -
第四，“资产回报率(roa)”分组箱线图显示，非ST组的中位数和均值都要明显高于ST组，说明资产回报率低的公司更容易被ST。 -
第五，“销售收入增长率(growth)”分组箱线图显示，非ST组的中位数和均值都要明显高于ST组，说明销售收入增长率低的公司更容易被ST。 -
第六，“杠杆水平(lev)”分组箱线图显示，非ST组的中位数和均值都要明显低于ST组，说明杠杆水平高的公司更容易被ST。 -
第七，“第一大股东持股比率(share)”分组箱线图显示，非ST组的中位数和均值都要明显高于ST组，说明第一大股东持股比率低的公司更容易被ST。
由此可见，各个自变量似乎跟是否ST都有一些相关性，但是，这种相关性是否具有显著的统计意义，需要后续使用正式模型进行分析验证。
6.逻辑回归
6.1筛选解释变量
6.1.1输入法筛选
# 确定自变量数据x_enter=df[['ara','asset','ato','roa','growth','lev','share']] # 确定因变量数据y_enter=df.st # 加上一列全为1的数据，使得模型矩阵中包含截距X_enter=sm.add_constant(x_enter) # 用未标准化数据拟合模型：X大写则拟合含截距模型，x小写则拟合不含截距项模型，# 标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致model_enter=sm.Logit(y_enter,X_enter).fit() # 和使用sklearn的结果不一样，# 因为和sklearn的损失函数不一样，sklearn的损失函数加入了正则项pd.DataFrame({'pvalue':model_enter.pvalues}).T Optimization terminated successfully. Current function value: 0.183849 Iterations 8结果解读：
最终只有
ara(p=0.001077)
和
lev(p=0.050050)
进入模型，其它自变量被剔除出模型。
6.1.2步进法筛选
# 确定自变量数据x_step=df[['ara','asset','ato','roa','growth','lev','share']] # 确定因变量数据y_step=df.st """ 使用步进法筛选出的自变量， stepwise_select_variable 函数的功能是使用步进法筛选自变量， 并返回最终筛选出的自变量名称"""x_step=df[stepwise_select_variable(x_step,y_step)] # 加上一列全为1的数据，使得模型矩阵中包含截距X_step=sm.add_constant(x_step) """ 用未标准化数据拟合模型： X大写则拟合含截距模型， x小写则拟合不含截距项模型， 标准化与未标准化模型的x大小写需一致， 否则检验统计量取值将不一致"""model_step=sm.Logit(y_step,X_step).fit() 上下滑动查看更多
结果解读：
最终只有ara和lev进入模型，其它自变量被剔除出模型，和使用“输入法”选择自变量的结果一致。
6.1.3子集筛选
上下滑动查看更多源码
# 确定包含所有自变量的数据x_all=df[['ara','asset','ato','roa','growth','lev','share']] # 确定因变量数据y_all=df.st subset_score=pd.DataFrame( columns=['independent_variable','AIC','BIC'])# 用于确定往subse_score数据集的第几行添加数据i=0# 确定自变量组合中自变量的数量for x_quantity in range(1,x_all.shape[1]+1): # 遍历自变量数量为x_quantity的各种自变量组合for x_com in combinations(x_all.columns,x_quantity): x=df[list(x_com)] X=sm.add_constant(x) model_all=sm.Logit(y_all,X).fit() subset_score.loc[i]=[x_com,model_all.aic, model_all.bic] i+=1display(subset_score)print('AIC最小值为{}，对应的自变量组合为{}。' .format(subset_score.AIC.min(), subset_score.loc[subset_score.AIC.argmin(),'independent_variable'])) # argmin()函数返回最小值对应的索引号print('BIC最小值为{}，对应的自变量组合为{}。' .format(subset_score.BIC.min(), subset_score.loc[subset_score.BIC.argmin(),'independent_variable'])) # argmin()函数返回最小值对应的索引号# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
结果解读：
此表显示了自变量的所有组合对应的AIC值、BIC值：根据AIC取值最小自变量筛选原则，所筛选出的自变量为“应收账款占比(ara)”、“销售收入增长率(growth)”、“杠杆水平(lev)”；根据BIC取值最小自变量筛选原则，所筛选出的自变量为“应收账款占比(ara)”。
6.2拟合回归模型
6.2.1全模型
print(model_enter.summary())6.2.2step模型
print(model_step.summary())6.2.3 AIC模型
# 确定自变量数据x_aic=df[['ara','growth','lev']] # 确定因变量数据y_aic=df.st # 加上一列全为1的数据，使得模型矩阵中包含截距X_aic=sm.add_constant(x_aic) """用未标准化数据拟合模型： X大写则拟合含截距模型，x小写则拟合不含截距项模型，标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致"""model_aic=sm.Logit(y_aic,X_aic).fit()"""和使用sklearn的结果不一样，因为和sklearn的损失函数不一样，sklearn的损失函数加入了正则项"""print(model_aic.summary()) 结果解读：
-
第一，常数项，取值为-4.6022，表示当ara(应收账款占比)、growth(销售收入增长率)、lev(杠杆水平)同时为0时，公司被ST概率是不被ST概率的exp(-4.6022)=0.01倍。 -
第二，回归系数，以ara为例，ara回归系数为5.1301，表示在其它条件保持不变的情况下，ara每增加1单位(即100个百分点)，公司被ST概率的发生比是原来的exp(5.1301)=169.034倍。那么，如果应收账款占比ARA每增加1个百分点，公司被ST概率的发生比是原来的多少倍？这时，可能有人会觉得是169.034除以100，也就是1.69倍，但是，我要告诉大家的是，这个答案是错误的。因为，当自变量ara增加1个百分点(即1%个单位)时，公司被ST概率的发生比是原来的exp(5.1301/100)=1.05倍，而不是1.69倍。
6.2.4BIC模型
# 确定自变量数据x_bic=df[['ara']] # 确定因变量数据y_bic=df.st # 加上一列全为1的数据，使得模型矩阵中包含截距X_bic=sm.add_constant(x_bic) """用未标准化数据拟合模型： X大写则拟合含截距模型，x小写则拟合不含截距项模型，标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致"""model_bic=sm.Logit(y_bic,X_bic).fit() """和使用sklearn的结果不一样，因为和sklearn的损失函数不一样，sklearn的损失函数加入了正则项"""print(model_bic.summary()) 6.3评估模型性能
6.3.1预测准确度
# 计算预测准确率(全模型)# 确定自变量数据x_full=df[['ara','asset','ato', 'roa','growth','lev','share']]# 确定因变量数据 y_full=df.st # 加上一列全为1的数据，使得模型矩阵中包含截距X_full=sm.add_constant(x_full) """用未标准化数据拟合GLM广义线性回归模型(连接函数为logit)： X大写则拟合含截距模型，x小写则拟合不含截距项模型，标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致"""model_full=sm.GLM( y_full, X_full, family=sm.families.Binomial(sm.families.links.logit()) ).fit()# 输出将来被ST的预测概率y_full_proba_pred=model_full.predict(X_full) # 输出将来被ST的预测分类y_full_pred=y_full_proba_pred.apply(lambda z:0 if z<0.5 else 1) # 计算预测准确率accuracy_score(y_full,y_full_pred) 0.9473684210526315结果解读：
全模型的预测准确率为94.7368%，但是，这个预测准确率并不可靠，因为如果把所有公司都预测为非ST，预测准确率都有94.7368%，当结果变量数据分布越不均衡时，全部预测为某种类型时的预测准确率越高。显然，这种高准确率没有任何商业意义，因为它一个ST公司都没有预测出来。在商业实战中，我们更关心两个指标：命中率（真正率）、假警报率（假正率）。
6.3.2ROC曲线
1、ROC曲线的基本原理
目的：用于选择模型，也就是选择AUC值最大的模型。
# 混淆矩阵定义( pd.DataFrame(index=['实际未发生_0','实际发生_1']) .assign(预测未发生_0=['True Negative(TN)正确否定', 'False Positive(FN)漏报']) .assign(预测发生_1=['False Positive(FP)虚报', 'True Positive(TP)正确肯定']))结果解读：
-
TP：预测发生且实际发生的数量，即正确肯定数量。 -
FP：预测未发生但实际发生的数量，即漏报数量。 -
FP：预测发生但实际未发生的数量，即虚报数量。 -
TN：预测未发生且实际未发生的数量，即正确否定数量。 -
命中率(TPR,True Positive Rate)，计算的是所有实际发生的客户中被预测为发生的客户所占的比例，也称为真正率和召回率，等于TP/(TP+FN)。 -
假警报率(FPR,False Positive Rate)，计算的是所有实际未发生的客户中被预测为发生的客户所占的比例，也称为假正率，等于FP/(FP+TN)。
# 绘制混淆矩阵(全模型)pd.DataFrame(confusion_matrix(y_full,y_full_pred), index=['实际非ST','实际ST'], columns=['预测非ST','预测ST'])一个优秀的预测模型，命中率(TPR)应尽可能高，即能尽量识别潜在被ST的公司，同时假警报率(FPR)也很低，即不要把不会被ST的公司误判为会被ST的公司。然而，TPR和FPR往往正相关，因为如果调高阈值，假如认为被ST的概率超过90%才被认定为将来会被ST，那么会导致假警报率很低，而如果调低阈值，例如认为被ST的概率超过10%就认定为将来会被ST，那么命中率就会很高，但是假警报率也会很高。
因此，为了衡量一个模型的优劣，统计学家根据不同阈值下的命中率和假警报率绘制了ROC曲线，该曲线以假警报率(FPR)为横轴，以命中率(TPR)为纵轴。
在不同阈值条件下，都希望命中率尽可能高，假警报率尽可能低。例如，某样本量为1000，其中ST公司为200家，非ST公司800家，当阈值为20%时，也就是说，被ST的概率超过20%时就预测其将来会被ST，模型A和模型B预测出来的ST公司都是150家。
如果模型A预测被ST的150家公司中有100家确实被ST，有50家属于误判，那么命中率等于100/200=50%，此时假警报率等于50/800=6.25%；如果模型B预测被ST的150家公司中有50家确实被ST，有100家属于误判，那么命中率等于50/200=25%，此时假警报率等于100/800=12.5%。此时，模型A的命中率是模型B的2倍，假警报率是模型B的一半，因此，可认为模型A是一个较优的模型。
如果把假警报率理解为代价的话，那么命中率就是收益，所以也可以说在阈值相同的情况下，希望假警报率(代价)尽可能小，命中率(收益)尽可能大，该思想反应在图形上就是ROC曲线尽可能陡峭。曲线越靠近左上角，说明在相同的阈值条件下，命中率(收益)越高，假警报率(代价)越低，模型越完善。换一个角度理解，一个完善的模型是在不同的阈值条件下，假警报率(代价)都接近于0，而命中率(收益)都接近于1，该特征反映在图形上就是ROC曲线非常接近(0,1)这个点，即曲线非常陡峭。
# 输出精准率、命中率(召回率)、f1-score、样本量、准确率(全模型)print(classification_report(y_full,y_full_pred,digits=5)) precision recall f1-score support 0 0.94868 0.99846 0.97293 648 1 0.50000 0.02778 0.05263 36 accuracy 0.94737 684 macro avg 0.72434 0.51312 0.51278 684weighted avg 0.92507 0.94737 0.92450 684结果解读：
-
support：样本数量，取值为0(非ST)的样本数量为648家，取值为1(ST)的样本数量为36家，684为总样本量。 -
accuracy：预测准确率=(TP+TN)/(TP+TN+FP+FN)=(647+1)/(647+0+36+1)=94.737%。 -
recall：命中率(召回率)，"ST"类别的召回率=预测为ST且实际为ST的公司数量/实际为ST的公司数量=TP/(TP+FN)=1/(35+1)=2.778%；“非ST”类别的召回率=预测为非ST且实际为非ST的公司数量/实际为非ST的公司数量=TN/(TN+FP)=647/(647+1)=99.846%。 -
precision：精准率(不重要)，“ST”类别的精准率=预测为ST且实际为ST的公司数量/预测为ST的公司数量=TP/(TP+FP)=1/(1+1)=50.00%；“非ST”类别的精准率=预测为非ST且实际为非ST的公司数量/预测为非ST的公司数量=TN/(TN+FN)=647/(647+35)=94.868%。 -
f1-score(简单了解即可，不重要)：混合的度量，对不平衡类别比较有效，计算公式为2TP/(2TP+FP+FN)。
2、绘制ROC曲线
理解ROC曲线的原理后，接下来就要利用ROC曲线来评估刚刚搭建的全模型。
在商业实践中，我们希望在阈值相同的情况下，假警报率(代价)尽可能的小，命中率(收益)尽可能大，该思想反应在图形上，就是ROC曲线非常接近(0,1)，即曲线非常陡峭。然而，用曲线来描述会显得比较抽象，于是统计学家使用AUC值来衡量模型的好坏。AUC值(Area Under Curve)指ROC曲线下方的面积，该面积的取值范围通常为0.5~1，0.5表示随机判断，1表示完美模型。由于扰动因素的存在，在商业实践中，AUC值能达到0.75就可接受，如果能够达到0.85，就是非常不错的模型。
# 输出不同阈值(thres)条件下的命中率TPR、假警报率FPR(全模型)fpr_fulls,tpr_fulls,thres_fulls=roc_curve(y_full,y_full_proba_pred) # y_full表示是否被ST的真实结果，y_full_proba_pred表示将来被ST的预测概率，thres_fulls表示阈值( pd.DataFrame() .assign(阈值=thres_fulls) .assign(假警报率FPR=fpr_fulls) .assign(命中率TPR=tpr_fulls))结果解读：
可见，随着阈值的下降，命中率和假警报率都在上升，对以上表格，有如下几点说明(KS值在后面再介绍)。
-
第一，表格第1行的阈值表示，只有当公司被ST的概率≥167.2906%，才预测其将来会被ST，但因为概率不会超过100%，所以此时所有公司都不会被预测为ST，此时命中率和假警报率都为0。可见，这个阈值其实没有什么实际意义，那么，为什么还要设置它呢？这个阈值是roc_curve()函数的默认设置，官方文档中对它的介绍是“thresholds[0] represents to instances being predicted and is arbitrarily set to max(y_score)+1"。此话意思是第一个阈值没有含义，其往往设置为最大阈值(本案例中为67.2906%)+1，以保证没有任何记录被选中。 -
第二，表格第2行的数据表示，只有当公司被ST的概率≥67.2906%时，才预测其将来会被ST，这个条件比较苛刻，此时被预测为被ST的公司很少，命中率为2.7778%，即”预测为被ST且实际被ST的公司数量/实际被ST的公司数量“这一比率为0.027778，具体到本案例可计算出，预测为被ST且实际被ST的公司数量=36*0.027778=1。此时假警报率为0，即”预测为被ST但实际未被ST的公司数量/实际未被ST的公司数量“这一比率为0，也就是说，实际未被ST的公司中不会有公司被误判为流失。 -
第三，表格后续数据的含义可参照第2行数据的含义来理解。
# 绘制ROC曲线(全模型)plt.plot(fpr_fulls,tpr_fulls)plt.plot([0,1],[0,1],color='red')plt.title('ROC')plt.xlabel('FPR')plt.ylabel('TPR')# 计算AUC值(全模型)"""y_full表示是否被ST的真实结果，y_full_proba_pred表示将来被ST的预测概率，y_full表示是否被ST的真实结果"""roc_auc_score(y_full,y_full_proba_pred) 0.7700188614540466AUC值的计算方法：
-
将坐标点按照横坐标FPR排序。 -
计算第i个坐标点和第i+1个坐标点的间距dx。 -
获取第i或者i+1个坐标点的纵坐标y。 -
计算面积微元ds=y*dx。 -
对面积微元进行累加，得到AUC。
6.3.3KS曲线
1、KS曲线的基本原理
目的：用于确定预测分类标准，即阈值thres，也就是选择KS值最大时对应的阈值。
KS曲线和ROC曲线本质上是相同的，同样关注命中率(TPR)和假警报率(FPR)，希望命中率(TPR)尽可能高，即可能找出潜在可能被ST的公司，同时也希望假警报率(FPR)尽可能低，即不要把未被ST的公司误判为被ST的公司。
KS曲线与ROC曲线的区别：ROC曲线将假警报率(FPR)作为横轴，将命中率(TPR)作为纵轴；KS曲线将阈值作为机轴，将命中率(TPR)与假警报率(FPR)之差作为纵坐标。
和ROC曲线一样，除了可视化的图形外，还需要一个可量化的指标来衡量模型预测效果，与ROC曲线相对应的是AUC值，而与KS曲线相对应的则是KS值。
KS值的计算公式：KS=max(TPR-FPR)。
KS值就是KS曲线的峰值。具体来说，每个阈值都对应着一个(TPR-FPR)值，那么，一定存在一个阈值，使得在该阈值条件下，(TPR-FPR)值最大，那么，此时的(TPR-FPR)值便称为KS值。一般来说，我们希望模型有较大的KS值，因为较大的KS值说明模型有较强的区分能力。不同取值范围的KS值的含义如下：
-
KS值小于0.2，一般认为模型的区分能力较弱； -
KS值在[0.2～0.3]区间内，模型具有一定的区分能力； -
KS值在[0.3～0.5]区间内，模型具有较强的区分能力；
但是，KS值也不是越大越好，如果KS值大于0.75，往往表示模型有异常。在商业实践中，KS值处于[0.2～0.3]区间内就已经挺不错。
2、绘制KS曲线
# 绘制KS曲线plt.plot(thres_fulls[1:],tpr_fulls[1:], label='tpr',color='green')plt.plot(thres_fulls[1:],fpr_fulls[1:], label='fpr',color='blue')plt.plot(thres_fulls[1:], tpr_fulls[1:]-fpr_fulls[1:], label='tpr-fpr', color='red')plt.xlabel('threshold')plt.legend()"""把阈值从小到大排序再绘制KS曲线，gca()函数获取坐标信息，invert_xaxis()函数反转x轴，纯属为了让可视化图形美观一些"""plt.gca().invert_xaxis() #plt.axvline(x=thres,color='black',linestyle='--',linewidth=0.5)# 计算全模型的KS值max(tpr_fulls-fpr_fulls)0.4367283950617284结果解读：
全模型的KS值为0.4367，处于[0.3～0.5]区间，认为全模型具有较强的区分能力。
# 计算全模型KS值对应的thres(阈值)、TPR、FPR、混淆矩阵、预测准确率for thres_full,tpr_full,fpr_full in zip(thres_fulls,tpr_fulls,fpr_fulls): if tpr_full-fpr_full==max(tpr_fulls-fpr_fulls): y_full_pred_thres=y_full_proba_pred.apply(lambda z:0 if z<thres_full else 1) # 全模型KS值对应的预测分类 print(""" 全模型KS值为{}， \n 其对应的阈值为{}、 \n 命中率TPR为{}、 \n 假警报率FPR为{}、 \n 预测准确率为{}。 """ .format(max(tpr_fulls-fpr_fulls), thres_full, tpr_full, fpr_full, accuracy_score(y_full, y_full_pred_thres))) display(pd.DataFrame(confusion_matrix(y_full, y_full_pred_thres), index=['实际非ST','实际ST'], columns=['预测非ST','预测ST'])) print(classification_report(y_full, y_full_pred_thres, digits=5)) # 输出全模型KS值取值最大时的精准率、命中率(召回率)、f1-score、样本量全模型KS值为0.4367283950617284，其对应的阈值为0.05538975001766095、命中率TPR为0.6944444444444444、假警报率FPR为0.25771604938271603、预测准确率为0.7397660818713451。 precision recall f1-score support 0 0.97764 0.74228 0.84386 648 1 0.13021 0.69444 0.21930 36 accuracy 0.73977 684 macro avg 0.55393 0.71836 0.53158 684weighted avg 0.93304 0.73977 0.81099 6846.3.4模型性能比较
在前面，我们只计算了全模型的各项性能指标，现在我们需要综合比较全模型、step模型、AIC模型、BIC模型的优劣。
上下滑动查看更多源码
model_dictory={'model_full':['ara','asset','ato', 'roa','growth','lev','share'], 'model_step':['ara','lev'], 'model_AIC':['ara', 'growth', 'lev'], 'model_BIC':['ara']}# 存储预测准确率、AUC值、KS值及其对应的阈值、TPR、FPR、预测准确率model_score=pd.DataFrame(columns=[ '阈值0.5时的预测准确率', 'AUC值', 'KS值','KS值对应的阈值', 'KS值对应的TPR命中率', 'KS值对应的FPR假警报率', 'KS值对应的预测准确率'])plt.figure(figsize=(16,8))i=241 # 2表示将画布分为2行，4表示将画布分为4列，1表示第1个子图# 遍历字典的键和值，IDV=indepenfor key,IDV in model_dictory.items(): dent variable x_temp=df[IDV] # 确定自变量数据 y_temp=df.st # 确定因变量数据 # 加上一列全为1的数据，使得模型矩阵中包含截距 X_temp=sm.add_constant(x_temp) """ 用未标准化数据拟合GLM广义线性回归模型(连接函数为logit)： X大写则拟合含截距模型， x小写则拟合不含截距项模型， 标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致 """ model_temp=sm.GLM(y_temp,X_temp, family=sm.families.Binomial(sm.families.links.logit()) ).fit() # 输出被ST的预测概率 y_temp_proba_pred=model_temp.predict(X_temp) # 输出将来被ST的预测分类(以0.5为判别标准) y_temp_pred=y_temp_proba_pred.apply(lambda z:0 if z<0.5 else 1) print('{}混淆矩阵(阈值为0.5,即被ST的预测概率大于等于0.5就判别为“ST”，否则判别为“非ST”)\n'.format(key)) display(pd.DataFrame(confusion_matrix(y_temp,y_temp_pred), index=['实际非ST','实际ST'], columns=['预测非ST','预测ST'])) # 混淆矩阵 # 输出精准率、命中率(召回率)、f1-score、样本量 print(classification_report(y_temp,y_temp_pred,digits=6)) """ "-"表示用于填充的字符，“<”表示左对齐，55表示填充满30个字符 往model_score数据集，添加预测准确率、AUC值、 KS值及其对应的阈值、TPR命中率、FPR假警报率、预测准确率等数据记录 """ print(format('-','-<80')) """ 返回不同阈值(thres_temp)条件下的TPR、FPR， y_temp表示是否被ST的真实结果， y_temp_proba_pred表示将来被ST的预测概率 """ fpr_temps,tpr_temps,thres_temps=roc_curve( y_temp,y_temp_proba_pred) for thres_temp,tpr_temp,fpr_temp in zip(thres_temps,tpr_temps,fpr_temps): if tpr_temp-fpr_temp==max(tpr_temps-fpr_temps): y_temp_pred_thres=y_temp_proba_pred.apply(lambda z:0 if z<thres_temp else 1) # 全模型KS值对应的预测分类 model_score.loc[key]=[accuracy_score(y_temp,y_temp_pred), # 阈值为0.5时的预测准确率 roc_auc_score(y_temp,y_temp_proba_pred), # AUC值 tpr_temp-fpr_temp, # KS值 thres_temp, # KS值对应的阈值 tpr_temp, # KS值对应的命中率、召回率 fpr_temp, # KS值对应的假警报率 accuracy_score(y_temp,y_temp_pred_thres) # KS值对应的预测准确率 ] # 可视化-绘制ROC曲线 plt.subplot(i) plt.plot(fpr_temps,tpr_temps) plt.plot([0,1],[0,1],color='red') plt.title('ROC') plt.xlabel('FPR({})'.format(key)) plt.ylabel('TPR({})'.format(key)) # 可视化-绘制KS曲线 plt.subplot(i+4) plt.plot(thres_temps[1:], tpr_temps[1:], label='tpr({})'.format(key), color='green') plt.plot(thres_temps[1:], fpr_temps[1:], label='fpr({})'.format(key), color='blue') plt.plot(thres_temps[1:], tpr_temps[1:]-fpr_temps[1:], label='tpr-fpr({})'.format(key),color='red') plt.xlabel('threshold({})'.format(key)) plt.legend() plt.gca().invert_xaxis() """ 把阈值从小到大排序再绘制KS曲线， gca()函数获取坐标信息， invert_xaxis()函数反转x轴，纯属为了让可视化图形美观一些 """ i=i+1display(model_score)# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
上下滑动查看更多
# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
结果解读：
由于AIC模型的AUC值和KS值都最大，所以最终选择AIC模型，自变量包括ara、growth、lev等3个。
6.3.5方差膨胀因子
多元线性回归模型一般需要进行多重共线性诊断、方差齐次性诊断、序列相关性诊断，但是，逻辑回归由于其目标变量取值只有0和1两种，必然违背方差齐假设，所以，逻辑回归无须进行方差齐次性诊断，至于序列相关性诊断则一般通过定性判断样本是否独立来完成，因此，逻辑回归只需进行多重共线性诊断即可。
多重共线性指自变量之间是否存在线性关系，与因变量无关，换句话说因变量不会影响自变量之间的多重共线性，因此，逻辑回归的多重共线性诊断所使用的方法与线性回归完全相同。
# 计算VIFIDV=['ara','growth','lev'] # 确定自变量名称全集，IDV=independent variablevifs=[]for variable in IDV: x_vif=df[list(set(IDV)-{variable})] X_vif=sm.add_constant(x_vif) y_vif=df[variable] """ 用未标准化数据拟合模型： X大写则拟合含截距模型，x小写则拟合不含截距项模型， 标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致 """ model_vif=sm.OLS(y_vif,X_vif).fit() vif=1/(1-model_vif.rsquared) vifs.append(vif)VIFS=pd.DataFrame(index=IDV).assign(VIF=vifs).sort_values(by='VIF',ascending=False) # 按各变量vif值降序排序display(VIFS)结果解读：
当VIF<5时，回归方程存在轻度多重共线性；当5≤VIF<10时，回归方程存在较严重的多重共线性；当10≤VIF时，回归方程存在严重的多重共线性。
由于所有VIF都小于5，因此，自变量之间不存在多重共线性问题。
6.4探测交互效应
经过前面的分析，我们确定了ara、growth、lev等3个自变量进入最终模型，但是，还有必要从定量和定性的角度探测自变量之间是否存在交互效应(调节作用)，在这里只从定量的角度对自变量间的交互效应进行探测，因此，需要探测3对双向交互效应(统计上还存在三向交互效应，但是，在商业实践中，很少考虑三向交互效应)。
上下滑动查看更多源码
# 穷尽所有的两项交互效应IDV=['ara','growth','lev']interaction_pvalues=pd.DataFrame(columns=['pvalue'])# 遍历自变量数量为2的各种自变量组合for interact_var in combinations(IDV,2): """ 需使用[:]将df数据集全部复制， 否则会报“A value is trying to be set on a copy of a slice from a DataFrame”错。 """ x_temp=df[:][IDV] x_temp['{}*{}'.format(interact_var[0], interact_var[1]) ]=x_temp.loc[:,interact_var[0]] * \ x_temp.loc[:,interact_var[1]] # 添加交互项# 加上一列全为1的数据，使得模型矩阵中包含截距 X_temp=sm.add_constant(x_temp) y_temp=df.st # 确定因变量数据""" 用标准化数据拟合模型： X大写则拟合含截距模型， x小写则拟合不含截距项模型， 标准化与未标准化模型的x大小写需一致，否则检验统计量取值将不一致 """ model_temp=sm.Logit(y_temp,X_temp).fit() interaction_pvalues.loc[model_temp.pvalues.index[-1] ]=model_temp.pvalues[-1]display(interaction_pvalues)# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
结果解读：
交互效应探测结果显示ara(应收账款占比)与growth(销售收入增长率)存在双向交互效应。
6.5评估交互模型
6.5.1构建交互模型
x_interact=df[['ara','growth','lev']].assign(ara_growth=df.ara*df.growth)X_interact=sm.add_constant(x_interact)y_interact=df.stmodel_interact=sm.Logit(y_interact,X_interact).fit()print(model_interact.summary())6.5.2绘制ROC曲线
# 输出不同阈值(thres)条件下的命中率TPR、假警报率FPR(交互模型)y_interact_proba_pred=model_interact.predict(X_interact) # 计算被ST的预测概率fpr_interacts,tpr_interacts,thres_interacts=roc_curve(y_interact,y_interact_proba_pred) """y_interact表示被ST的真实结果，y_interact_proba_pred表示被ST的预测概率，thres_interact表示阈值"""( pd.DataFrame() .assign(阈值=thres_interacts) .assign(假警报率FPR=fpr_interacts) .assign(命中率TPR=tpr_interacts))# 绘制ROC曲线(交互模型)plt.plot(fpr_interacts,tpr_interacts)plt.plot([0,1],[0,1],color='red')plt.title('ROC')plt.xlabel('FPR')plt.ylabel('TPR')# 计算AUC值(交互模型)"""# y_interact表示是否被ST的真实结果，y_interact_proba_pred表示将来被ST的预测概率，y_interact表示是否被ST的真实结果"""roc_auc_score(y_interact,y_interact_proba_pred) 0.7872942386831276结果解读：
交互模型的AUC值(0.787)大于AIC模型的AUC值(0.774)，说明交互效应模型的性能优于AIC模型。由于扰动因素的存在，在商业实践中，AUC值能达到0.75就可接受，如果能够达到0.85，就是非常不错的模型。
6.5.3绘制KS曲线
# 绘制KS曲线plt.plot(thres_interacts[1:],tpr_interacts[1:],label='tpr',color='green')plt.plot(thres_interacts[1:],fpr_interacts[1:],label='fpr',color='blue')plt.plot(thres_interacts[1:],tpr_interacts[1:]-fpr_interacts[1:],label='tpr-fpr',color='red')plt.xlabel('threshold')plt.legend()plt.gca().invert_xaxis() # 把阈值从小到大排序再绘制KS曲线，gca()函数获取坐标信息，# invert_xaxis()函数反转x轴，纯属为了让可视化图形美观一些#plt.axvline(x=thres,color='black',linestyle='--',linewidth=0.5)# 计算交互模型KS值max(tpr_interacts-fpr_interacts)0.5138888888888888结果解读：
交互模型的KS值(0.514)小于AIC模型的KS值(0.520)，说明说明交互效应模型的性能劣于AIC模型。但是，由于KS值都大于0.5，可见交互模型与AIC模型都具有较强的区分能力。
上下滑动查看更多源码
# 计算交互模型KS值对应的thres(阈值)、TPR、FPR、混淆矩阵、预测准确率for thres_interact,tpr_interact,fpr_interact in zip(thres_interacts,tpr_interacts,fpr_interacts):if tpr_interact-fpr_interact==max(tpr_interacts-fpr_interacts): y_interact_pred_thres=y_interact_proba_pred.apply(lambda z:0if z<thres_interact else1) # 交互模型KS值对应的预测分类 print('交互模型KS值为{}， 其对应的阈值为{}、 命中率TPR为{}、 假警报率FPR为{}、 预测准确率为{}。' .format(max(tpr_interacts-fpr_interacts), thres_interact, tpr_interact, fpr_interact, accuracy_score(y_interact, y_interact_pred_thres) ) ) display( pd.DataFrame(confusion_matrix( y_interact, y_interact_pred_thres), index=['实际非ST','实际ST'], columns=['预测非ST','预测ST']))""" 输出交互模型KS值取值最大时的精准率、 命中率(召回率)、f1-score、样本量 """ print(classification_report( y_interact, y_interact_pred_thres, digits=5)) # 往model_score数据集添加交互模型的性能指标得分 model_score.loc['model_interact']=[# 阈值为0.5时的预测准确率 accuracy_score(y_interact, y_interact_proba_pred.apply(lambda z:0if z<0.5else1)),# AUC值 roc_auc_score(y_interact,y_interact_proba_pred), tpr_interact-fpr_interact, # KS值 thres_interact, # KS值对应的阈值 tpr_interact, # KS值对应的命中率、召回率 fpr_interact, # KS值对应的假警报率 accuracy_score(y_interact, y_interact_pred_thres) # KS值对应的预测准确率 ]display(model_score)# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
交互模型KS值为0.5138888888888888，其对应的阈值为0.05092930650051774、命中率TPR为0.8055555555555556、假警报率FPR为0.2916666666666667、预测准确率为0.7134502923976608。 precision recall f1-score support 0 0.98498 0.70833 0.82406 648 1 0.13303 0.80556 0.22835 36 accuracy 0.71345 684 macro avg 0.55900 0.75694 0.52620 684weighted avg 0.94014 0.71345 0.79270 684结果解读：
交互效应模型、AIC模型在AUC值、KS值、TRP命中率、FPR假警报率、KS值对应的预测准确率等5个模型性能指标上表现更有千秋，其中交互效应模型在AUC值、FPR假警报率、KS值对应的预测准确率等三个指标上胜出，其中AIC模型在KS值、TPR命中率等两个指标上胜出。由于交互效应模型所有回归系数均在0.05显著水平下显著，而AIC模型在0.05显著水平下却存在一个自变量不显著(growth的p值为0.098，大于0.05)，因此，最终使用交互效应模型，阈值确定为0.051。
6.6模型结果呈现
6.6.1交互效应结果解释
pd.DataFrame({'coef':round(model_interact.params,4), 'Std.Error':round(model_interact.bse,4), 'tvalue':round(model_interact.tvalues,3), 'pvalue':round(model_interact.pvalues,3), 'lower 95% CI':round(model_interact.conf_int(0.05)[0],3), 'upper 95% CI':round(model_interact.conf_int(0.05)[1],3), 'exp(coef)':np.exp(model_interact.params), 'VIF':[None]+list(VIFS.VIF)+[None] })结果解读：
此模型为最终报告模型，包含交互项在内的所有回归系数都具有统计意义(p值都小于0.05)。
-
第一，常数项，取值为-4.6096，表示当ara(应收账款占比)、growth(销售收入增长率)、lev(杠杆水平)同时为0时，公司被ST概率是不被ST概率的exp(-4.6096)=0.01倍。 -
第二，ara回归系数为6.2503，表示在其它条件保持不变的情况下，growth(销售收入增长率)取值为0时，ara每增加1单位(即100个百分点)，公司被ST概率的发生比是原来的exp(6.2503)=518.155倍，ara增加1个百分点(即1%个单位)时，公司被ST概率的发生比是原来的exp(6.2503/100)=1.06倍，而不是5.18倍。 -
第三，growth回归系数为-2.5700，表示在其它条件保持不变的情况下，ara(应收账款占比)取值为0时，growth每增加1单位(即100个百分点)，公司被ST概率的发生比是原来的exp(-2.5700)=0.0765倍，growth增加1个百分点(即1%个单位)时，公司被ST概率的发生比是原来的exp(-2.5700/100)=0.9746倍。 -
第四，lev回归系数为2.3544，表示在其它条件保持不变的情况下，lev每增加1单位(即100个百分点)，公司被ST概率的发生比是原来的exp(2.3544)=14.692倍，lev增加1个百分点(即1%个单位)时，公司被ST概率的发生比是原来的exp(2.3544/100)=1.02倍。 -
第五，交互项回归系数为7.9980，说明伴随着应收账款占比(ara)增加的销售收入增长率(growth)增加会导致公司被ST风险的剧增，因为销售收入增长率增加带来的是应收账款增加，而不是实际现金流增加，这种伴随着高应收账款的销售增长，其增长质量很低甚至很可能是种做假账行为。当growth(销售收入增长率)为调节变量时，表示growth每增加1单位(即100个百分点)，ara(应收账款占比)自变量对logit(p)的影响额外增加7.9980个单位，即将ara对公司被ST概率的发生比的影响是原来的exp(7.9980)=2975.08倍，也就是说，同时具备高应收账款占比和高销售收入增长率的公司具有高ST风险。当ara(应收账款占比)为调节变量时，其解释与growth(销售收入增长率)为调节变量时的解释类别。
回归方程：
logit_p=-4.6096+6.2503*ara-2.57*growth+2.3544*lev+7.998*ara*growth
6.6.2交互效应可视呈现
上下滑动查看更多源码
# 交互效应可视化plt.figure(figsize=(16,6))plt.subplot(121) #growth为自变量，ara为调节变量for ara,label in zip(df.ara.quantile( [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]), ['min','10%','20%','30%','40%','50%','60%','70%','80%','90%','max'] ): logit_p=[]# 取min、25%、50%、75%、max等5个点for growth in df.growth.quantile([0,0.25,0.50,0.75,1]): logit_p.append(model_interact.params[0] \ +model_interact.params[1] \ *ara+model_interact.params[2] \ *growth+model_interact.params[4] \ *ara \ *growth) plt.plot(df.growth.describe()[3:8], logit_p, label=label)plt.xlabel('销售收入增长率(growth)')plt.ylabel('Logit(P)')plt.legend(ncol=4,loc='upper left')plt.subplot(122) #ara为自变量，growth为调节变量for growth,label in zip(df.growth.quantile( [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]), ['min','10%','20%','30%','40%','50%','60%','70%','80%','90%','max'] ): logit_p=[]# 取min、25%、50%、75%、max等5个点for ara in df.ara.quantile([0,0.25,0.50,0.75,1]): logit_p.append(model_interact.params[0] \ +model_interact.params[1] \ *ara+model_interact.params[2] \ *growth+model_interact.params[4] \ *ara \ *growth) plt.plot(df.ara.describe()[3:8],logit_p,label=label)plt.xlabel('应收账款占比(ara)')plt.ylabel('Logit(P)')plt.legend(ncol=4,loc='upper left')# 关注和星标@公众号：数据STUDIO，更多优质内容等你发现～
🏴☠️宝藏级🏴☠️ 原创公众号『 数据STUDIO 』内容超级硬核。公众号以Python为核心语言，垂直于数据科学领域，包括 可戳 👉 Python ｜ MySQL ｜ 数据分析 ｜ 数据可视化 ｜ 机器学习与数据挖掘 ｜ 爬虫 等，从入门到进阶！
长按👇关注- 数据STUDIO -设为星标，干货速递