1、 请采用计数数据分析模型(Count Data Model),对Crash Frequency.xls文件的数据进行建模分析,并回答以下问题:

  1. 所采用计数数据分析模型的类型及原因;
  2. 事故发生频率显著相关的因素;
  3. 模型整体拟合优度的评价。

其中,Grade_G为分类变量,代表路段坡度;其余变量为连续变量,分别代表交通周转量(Logvmt)、平均能见度(Visibility)、平均气温(Temperature)、小时降雨量(Precipitation)和平均速度(Speed)。

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.simplefilter("ignore") 
sns.set_style('whitegrid')

1.1、首先导入相关数据

df_1 = pd.read_excel('Crash frequency.xls')

1.2、描述性统计

可以看到,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_summary
summary_statistics(df_1)

在这里插入图片描述

1.3、计算变量的方差膨胀因子,检查多重共线性问题

可以看到,所有变量的方差膨胀因子(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()

在这里插入图片描述

1.4、采用负二项回归建模

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。

1.5、拟合泊松回归模型

为了更好地评价负二项回归模型的拟合优度,我们利用相同的变量拟合泊松回归模型

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岁的男性沪牌车驾驶员(未载客)其闯红灯的概率(给出计算公式)。

2.1、首先导入相关数据

df_2 = pd.read_excel('Red light running.xls')

2.2、生成交叉口和年龄的哑变量

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 dummy

df_2 = pd.concat([df_2,get_dummy_var(df_2,'Intersection'),get_dummy_var(df_2,'Age')],axis=1)

2.3、描述性统计

可以看到,三号交叉口和40岁年龄的比例最高,因此在建模时讲其设置为参照项。

summary_statistics(df_2)

在这里插入图片描述

2.4、检查共线性

所有变量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'])

在这里插入图片描述

2.5、建立二项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

OR值(odds ratio)又称比值比、优势比,主要指病例组中暴露人数与非暴露人数的比值除以对照组中暴露人数与非暴露人数的比值,是流行病学研究中病例对照研究中的一个常用指标。

参考:https://www.douban.com/note/352258282/

举例:
共有10个male,10个female考大学, 被录取的结果为:7个male, 3个female
使用上面的结论:
从录取这个角度而言,而言odds(male) = 0.7/0.3 = 2.33 odds(female) = 0.3/0.7 = 0.428
所以录取的odds ratio为:OR = odds(male)/odds(female) = 2.37/0.42=5.44;
结论和意义:对一个male而言,录取的成功率比femal 高5.44倍。

from math import e,log

plt.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、 请基于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)。

3.1、首先导入相关数据

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

3.2、描述性统计

可以看到,仅有7.3%的样本为受伤事故,在该场景下,我们希望尽可能预测出受伤事故,因此我采用召回率(Recall)作为模型性能的评价指标,以便更多的找到受伤事故样本。此外,对于所有样本,变量ice均为no,因此在建模过程中,我舍弃了该变量。

summary_statistics(df_4)

在这里插入图片描述

3.3、开始构建决策树

from sklearn.tree import DecisionTreeClassifier,export_graphviz
from sklearn.model_selection import GridSearchCV,train_test_split
import graphviz

x = df_4.drop(columns=['severity','ice'])
y = df_4['severity']
name = x.columns

3.4、利用网格搜索,标定决策树最优超参数

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

在这里插入图片描述

3.5、获得最优超参数:split criterion为gini,最大深度为9,最小划分样本数量为2

print(optimized_GBM.best_params_)

结果:{'criterion': 'gini', 'max_depth': 9, 'min_samples_split': 2}

3.6、构建决策树

clf = DecisionTreeClassifier(criterion='gini', max_depth=9, min_samples_split=2,random_state=10)
clf.fit(x,y)

3.7、画出决策树结构图

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

在这里插入图片描述

3.8、解读决策树分类规则

在这里插入图片描述对上述红、黄、蓝三个分类路径进行解读:

  • 首先决策树根据运行速度方差进行样本划分,如果运行速度方差大于1.515,样本进入右侧分类路径,这同样符合常理,当道路运行速度方差大时,说明低速高速车辆同时存在,发生车祸时危险程度更大;
  • 对于右侧样本,继续根据路肩宽度进行划分,如果路肩宽度小于等于2,样本进入左侧划分区域,如果路肩宽度大于2,样本则进入右侧划分区域;
  • 对于速度方差大于1.515、路肩宽度小于等于2的样本,如果路面湿滑,则被预测为PDO样本,如果路面不湿滑且运行速度方差小于等于1.726,则被预测为PDO样本;
  • 对于速度方差大于1.515、路肩宽度大于2的样本,如果温度小于-9,则被预测为INJ样本。

4、 请基于Severity.csv文件,构建随机森林模型对事故严重程度的影响因素进行重要度排序:

  1. 给出随机森林设定参数;
  2. 给出变量重要度排序结果。

4.1、根据网格搜索选择最优参数

通常,针对随机森林,有三个重要参数,分别为:弱学习器(决策树)数量,弱学习器(决策树)的深度(代表了弱学习器的复杂程度)以及每次进行节点划分时使用的特征数量。下面,我将利用袋外样本recall指标,对这三个参数的最优组合进行标定。其中,弱学习器数量设置为[100,1600,100],最大特征数设置为[3,6,9],弱学习器深度设置为[1,50,5]

from sklearn.metrics import recall_score,accuracy_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance

gs = 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 tqdm

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

在这里插入图片描述

4.2、得到最优模型参数

100个弱学习器,最大特征数为3,决策树深度为16

rf = RandomForestClassifier(n_estimators=100,max_depth=16,max_features=3,random_state=20)
rf.fit(x,y)

4.3、计算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()

在这里插入图片描述

4.4、特征变量重要度排序结果

可以看出,总体上无论是impurity-based还是permutation方法计算的特征重要度,速度的标准差和温度都是最重要的两个变量,其次是能见度、降雨量、中央隔离带宽度和是否陡坡;在两种方式计算的特征重要度中,限速、车道数、货车比例和日均流量对数均为四个特征重要度最小的变量。

参考:

https://mp.weixin.qq.com/s/3DHEAumY0F0K31Pb1TjruQ

Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐