手把手教你用Python构建logit、负二项回归、决策树与随机森林机器学习模型
本次更新的主要内容为利用Python中的statsmodels库构建logit与负二项回归模型，以及利用sklearn库构建决策树以及随机森林模型。内容源自同济大学研究生课程《高级数理统计》的第一次大作业，一共有五道题，这里我把题目和我的答案都贴上来，仅供参考，不一定对，大家可以通过公众号直接后台私信我，多多互相交流。
1. 请采用计数数据分析模型（Count Data Model），对Crash Frequency.xls文件的数据进行建模分析，并回答以下问题：
（1）所采用计数数据分析模型的类型及原因；
（2）事故发生频率显著相关的因素；
（3）模型整体拟合优度的评价。
其中，Grade_G为分类变量，代表路段坡度；其余变量为连续变量，分别代表交通周转量（Logvmt）、平均能见度（Visibility）、平均气温（Temperature）、小时降雨量（Precipitation）和平均速度（Speed）。
Answer:
import pandas as pdimport numpy as npimport seaborn as snsimport matplotlib.pyplot as pltimport statsmodels.api as smfrom statsmodels.stats.outliers_influence import variance_inflation_factorimport warningswarnings.simplefilter('ignore') sns.set_style('whitegrid')首先导入相关数据
df_1 = pd.read_excel('Crash frequency.xls')描述性统计
可以看到，0值在Crash_Freq中的分布大于50%，且均值（1.096）远小于方差（3.802），可以考虑建立零膨胀模型。
def summary_statistics(df): df_summary = df.describe().T df_summary['std'] = df.var() df_summary.columns = ['count', 'mean', 'var', 'min', '25%', '50%', '75%', 'max'] return df_summarysummary_statistics(df_1)
计算变量的方差膨胀因子，检查多重共线性问题
可以看到，所有变量的方差膨胀因子（VIF）均小于5，说明没有多重共线性问题存在。
def cal_vif(df,col): df = sm.add_constant(df) col.append('const') df_vif = df.loc[:,col] return pd.DataFrame({'variable':col,'VIF':[variance_inflation_factor(np.matrix(df_vif),i) for i in range(len(col))]}).iloc[:-1]cal_vif(df_1,['Grade_G','Logvmt', 'Visibility', 'Temperature', 'Precipitation', 'Speed'])
plt.figure(figsize=(10,6))sns.histplot(x='Crash_Freq',data=df_1,binwidth=0.5,color='blue')plt.xticks(range(0,df_1['Crash_Freq'].max()+1),fontsize=15)plt.yticks(fontsize=15)plt.xlabel('Crash_Freq',fontsize=16)plt.ylabel('Count',fontsize=16)plt.show()在stata中导入数据后输入下面代码建立零膨胀负二项模型
zinb Crash_Freq Grade_G-Speed, inflate(Grade_G Logvmt) vuong
可以看到，Vuong检验的结果只在10%水平下显著（轻微显著），为了保证模型的简洁性，我们采用负二项回归重新建模。
model_1_nb = sm.NegativeBinomial.from_formula('Crash_Freq ~ Grade_G + Logvmt + Visibility + Temperature + \ Precipitation + Speed ',data=df_1).fit(method='ncg',maxiter=1000)print(model_1_nb.summary2())模型结果表明
事故发生频率与交通周转量及路段坡度显著正相关（即交通周转量越大、坡度越大的路段，事故发生的频率越高）。
事故发生频率与平均能见度、平均气温以及平均速度显著负相关 （即平均能见度越低、气温越低以及平均速度越低的路段，事故发生率越高）。
模型拟合优度方面
离散系数alpha在0.1%水平下显著，表明有必要采用负二项回归。
模型的AIC与BIC分别为579.29和607.13，伪R方为0.180。
为了更好地评价负二项回归模型的拟合优度，我们利用相同的变量拟合泊松回归模型
model_1_poisson = sm.Poisson.from_formula('Crash_Freq ~ Grade_G + Logvmt + Visibility + Temperature + \ Precipitation + Speed ',data=df_1).fit(method='ncg',maxiter=1000)print(model_1_poisson.summary2())可以发现泊松回归的AIC与BIC都明显大于负二项回归的AIC与BIC，再次表明负二项回归的拟合优度更高。
2. Red light running.xls文件是研究人员对四个交叉口开展闯红灯调查的记录数据。
其中，Running为二进制变量，表明观测对象是否闯红灯（1-闯红灯，0-未闯红灯）；Intersection为分类变量，表明交叉口编号；Local为二进制变量，表明观测对象是否是沪牌车（1-沪牌车，0-外地牌照车）；Passenger为二进制变量，表明观测对象是否载有乘客（1-载有乘客，0-未载有乘客）；Male为二进制变量，表明驾驶人是否为男性（1-男性，0-女性）；Age为分类变量，表明驾驶人的年龄分组。请基于此数据，采用恰当的分析模型，回答以下问题：
（1）闯红灯行为显著相关的变量；
（2）4个交叉口间闯红灯行为是否有显著差异；
（3）计算在交叉口1，一位年龄为45岁的男性沪牌车驾驶员（未载客）其闯红灯的概率（给出计算公式）。
Answer:
首先导入相关数据
df_2 = pd.read_excel('Red light running.xls')生成交叉口和年龄的哑变量
def get_dummy_var(df,col): dummy = pd.get_dummies(df[col]) new_col = [] for i in dummy.columns: new_col.append(col + '_' + str(i)) dummy.columns = new_col return dummydf_2 = pd.concat([df_2,get_dummy_var(df_2,'Intersection'),get_dummy_var(df_2,'Age')],axis=1)
描述性统计
可以看到，三号交叉口和40岁年龄的比例最高，因此在建模时讲其设置为参照项。
summary_statistics(df_2)检查共线性
所有变量VIF均小于5，表明没有多重共线性存在。
cal_vif(df_2,['Local', 'Male', 'Passenger', 'Intersection_1', 'Intersection_2', 'Intersection_4', 'Age_20', 'Age_25', 'Age_30', 'Age_35', 'Age_38','Age_45', 'Age_50', 'Age_55'])
建立二项logistics模型
model_2_logit = sm.Logit.from_formula('Running ~ Local + Male + Passenger + Intersection_1 + Intersection_2 + Intersection_4 + \ Age_20 + Age_25 + Age_30 + Age_35 + Age_38 + Age_45 + Age_50 + Age_55',data=df_2).fit(method='ncg',maxiter=1000)print(model_2_logit.summary2())
通过模型结果可以发现，与闯红灯行为显著相关的变量有：
闯红灯行为与沪牌车、男性显著正相关，表明沪牌车与男性驾驶员更容易闯红灯。
闯红灯行为与有载客显著负相关，表明有载客的车辆更不容易闯红灯。
闯红灯行为与驾驶员年龄显著负相关，具体来讲，相比40岁的驾驶员，25和30岁的驾驶员更容易闯红灯，且25岁的驾驶员闯红灯的概率更高，而其余年龄的驾驶员在闯红灯行为方面与40岁的驾驶员没有显著区别。
4个交叉口间闯红灯行为是否有显著差异？回答：是
通过模型结果我们可以发现，Intersection_1, Intersection_2, Intersection_4三个变量均显著，这表明1，2，4号交叉口闯红灯行为与3号交叉口相比均有显著差异，进一步绘制这三个变量的Odds Ratio：
from math import e,logplt.scatter([-0.8803,-1.9788,0,0.4991],[4,3,2,1],color='k',marker='s')plt.hlines(4,-1.2981,-0.4625,color='k',linewidth=1.5)plt.hlines(3,-2.6077,-1.3500,color='k',linewidth=1.5)plt.hlines(1,-0.0146,1.0127 ,color='k',linewidth=1.5)for i in [log(0.1),log(0.5),log(1.5),log(2)]: plt.vlines(i,0,5,linestyle='--',color='grey',linewidth=0.5)plt.vlines(0,0,5,color='k',linewidth=0.5)plt.ylim(0.5,4.5)plt.grid(False)plt.yticks(range(1,5),['Intersection_4', 'Intersection_3', 'Intersection_2', 'Intersection_1'],fontsize=14)plt.xticks([log(0.1),log(0.5),log(1),log(1.5),log(2)],[0.1,0.5,1,1.5,2],fontsize=13)plt.xlabel('Odds Ratio',fontsize=14)plt.show()可以发现，四个交叉口发生闯红灯行为的概率大小为交叉口4>交叉口3>交叉口1>交叉口2。
（3）计算在交叉口1，一位年龄为45岁的男性沪牌车驾驶员（未载客）其闯红灯的概率（给出计算公式）。
3. 我觉得第三题的数据描述很不清楚，因此这道题跳过。
4. 请基于Severity.csv文件，针对事故严重程度（severity变量）构建决策树模型：
（1）给出Split Criterion计算依据
（2）给出决策树结构图；
（3）对模型分类规则进行解释（解读三条规则即可）。
变量解释：#### 事故严重程度（severity），PDO为仅财产损失事故、INJ为受伤事故；中央隔离带宽度（Med_Width）；路肩宽度（Inside_Shld）；限速（Speed_Limit）;货车比例（Truck_Per）；温度（temperature）；能见度（visibility）；1小时降雨量（1hourprecip）；运行速度方差（speed std）；日均流量对数（logaadt）；车道数（lane）；季节（snow_season）；路面结冰（ice）；路面湿滑（slush）；是否陡坡（steep grade）。
Answer:
首先导入相关数据
df_4 = pd.read_csv('severity binary.csv')### 变量编码df_4['severity'] = df_4['severity'].map({'PDO':0,'INJ':1})df_4['lane'] = df_4['lane'].map({'two':2,'thr':3})df_4['snow_season'] = df_4['snow_season'].map({'dry':0,'snow':1})df_4['ice'] = df_4['ice'].map({'no':0,'yes':1})df_4['slush'] = df_4['slush'].map({'no':0,'yes':1})df_4['steep grade'] = df_4['steep grade'].map({'no':0,'yes':1})描述性统计
可以看到，仅有7.3%的样本为受伤事故，在该场景下，我们希望尽可能预测出受伤事故，因此我采用召回率（Recall）作为模型性能的评价指标，以便更多的找到受伤事故样本。此外，对于所有样本，变量ice均为no，因此在建模过程中，我舍弃了该变量。
summary_statistics(df_4)开始构建决策树
from sklearn.tree import DecisionTreeClassifier,export_graphvizfrom sklearn.model_selection import GridSearchCV,train_test_splitimport graphviz
x = df_4.drop(columns=['severity','ice'])y = df_4['severity']name = x.columns利用网格搜索，标定决策树最优超参数
cv_params = {'min_samples_split':range(1,20), 'max_depth':range(1,20),'criterion':['gini','entropy']}ind_params = {'random_state': 10}optimized_GBM = GridSearchCV(DecisionTreeClassifier(**ind_params),cv_params,scoring='recall', cv=5, n_jobs=-1, verbose=10)optimized_GBM.fit(x, y)print('最优分数:', optimized_GBM.best_score_)获得最优超参数：split criterion为gini，最大深度为9，最小划分样本数量为2
print(optimized_GBM.best_params_){'criterion': 'gini', 'max_depth': 9, 'min_samples_split': 2}
构建决策树
clf = DecisionTreeClassifier(criterion='gini', max_depth=9, min_samples_split=2,random_state=10)clf.fit(x,y)画出决策树结构图
dot_data = export_graphviz(clf, out_file=None, feature_names=name, filled=True, rounded=True, class_names=['PDO','INJ'], special_characters=True) graph = graphviz.Source(dot_data) graph.render('决策树结构')graph解读决策树分类规则
对上述红、黄、蓝三个分类路径进行解读：
首先决策树根据运行速度方差进行样本划分，如果运行速度方差大于1.515，样本进入右侧分类路径，这同样符合常理，当道路运行速度方差大时，说明低速高速车辆同时存在，发生车祸时危险程度更大。
对于右侧样本，继续根据路肩宽度进行划分，如果路肩宽度小于等于2，样本进入左侧划分区域，如果路肩宽度大于2，样本则进入右侧划分区域。
对于速度方差大于1.515、路肩宽度小于等于2的样本，如果路面湿滑，则被预测为PDO样本，如果路面不湿滑且运行速度方差小于等于1.726，则被预测为PDO样本。
对于速度方差大于1.515、路肩宽度大于2的样本，如果温度小于-9，则被预测为INJ样本。
5. 请基于Severity.csv文件，构建随机森林模型对事故严重程度的影响因素进行重要度排序：
（1）给出随机森林设定参数；
（2）给出变量重要度排序结果。
根据网格搜索选择最优参数
通常，针对随机森林，有三个重要参数，分别为：弱学习器（决策树）数量，弱学习器（决策树）的深度（代表了弱学习器的复杂程度）以及每次进行节点划分时使用的特征数量。下面，我将利用袋外样本recall指标，对这三个参数的最优组合进行标定。其中，弱学习器数量设置为[100,1600,100]，最大特征数设置为[3,6,9]，弱学习器深度设置为[1,50,5]
from sklearn.metrics import recall_score,accuracy_scorefrom sklearn.ensemble import RandomForestClassifierfrom sklearn.inspection import permutation_importancegs = pd.MultiIndex.from_product([np.arange(100,1600,100),[3,6,9],np.arange(1,50,5)],names=['n_estimators','max_features','max_depth']).to_frame().reset_index(drop=True)gs['oob_score'] = 0### 开始网格搜索from tqdm import tqdmfor i in tqdm(gs.index):rf = RandomForestClassifier(n_estimators=gs.loc[i,'n_estimators'],max_features=gs.loc[i,'max_features'],max_depth=gs.loc[i,'max_depth'],oob_score=True,random_state=20,n_jobs=-1).fit(x,y.ravel())gs.loc[i,'oob_score'] = recall_score(y,[0 if i[0] >= 0.5 else 1 for i in rf.oob_decision_function_ ])
plt.figure(figsize=(10,6))plt.scatter(100,gs[gs['max_depth']==16]['oob_score'].max(),color='red',label='Optimal point',zorder=1)plt.legend()sns.lineplot(x='n_estimators',y='oob_score',data=gs[gs['max_depth']==16],hue='max_features',palette=['C0','C1','C2'],zorder=0)plt.xticks(range(100,1600,100),fontsize=12,rotation=45)plt.yticks(fontsize=12)plt.xlabel('n_estimators',fontsize=15)plt.ylabel('oob_score',fontsize=15)得到最优模型参数
100个弱学习器，最大特征数为3，决策树深度为16
rf = RandomForestClassifier(n_estimators=100,max_depth=16,max_features=3,random_state=20)rf.fit(x,y)
计算impurity-based特征重要度和permutation特征重要度
其中，impurity-based特征重要度根据每个节点根据特征划分后获得的基尼增益计算；permutation特征重要度为将特征随机打乱，计算打乱后模型性能的变化。
feature = pd.DataFrame(list(zip(name,rf.feature_importances_)),columns=['Variable','Relative Importance (impurity-based)']).sort_values(by='Relative Importance (impurity-based)')perm_imp = permutation_importance(rf,x,y,random_state=10,n_repeats=99)feature_perm = pd.DataFrame(list(zip(name,perm_imp['importances_mean'])),columns=['Variable','Relative Importance']).sort_values(by='Relative Importance')feature_perm['Relative Importance'] = feature_perm['Relative Importance']/feature_perm['Relative Importance'].sum()feature_perm.columns = ['Variable', 'Relative Importance (permutation)']feature = pd.merge(feature,feature_perm,on='Variable',how='left')plt.figure(figsize=(6,8))plt.barh(np.arange(13)+0.15,feature['Relative Importance (impurity-based)'],height=0.3,color='k')plt.barh(np.arange(13)-0.15,feature['Relative Importance (permutation)'],height=0.3,color='C3')plt.ylim(-0.5,12.5)plt.yticks(range(13),feature['Variable'],fontsize=13)plt.xticks(fontsize=14)plt.xlabel('Relative Importance',fontsize=16)plt.yticks(fontsize=16)plt.legend(['impurity-based','permutation'],fontsize=13)plt.show()特征变量重要度排序结果
可以看出，总体上无论是impurity-based还是permutation方法计算的特征重要度，速度的标准差和温度都是最重要的两个变量，其次是能见度、降雨量、中央隔离带宽度和是否陡坡；在两种方式计算的特征重要度中，限速、车道数、货车比例和日均流量对数均为四个特征重要度最小的变量。