Python-Statsmodels–出行行为分析
Python-Statsmodels–出行行为分析
问题解析
可以看到，一共有四问（有的逻辑我也没太懂，大概解读一下）：
- 估计只有一个自变量和常数项的二项logit模型，并解释系数、讨论相关的统计量，并选择合适的统计检验方式比较其与equally-likely model的优劣度。
- 计算equally-likely model和market share model的log-likelihood
- 估计一个正常的二项logit模型，并与market share模型进行比较
- 设计属于自己的二项logit模型（男女分开估计）
估计只有一个自变量和常数项的二项logit模型、计算EL和MS模型的Log-likelihood
这里我们先按照师兄的方法先把数据集清洗了一下(感谢瑞祥师兄
)，最终用来建模的数据集N=293，名称为model_data.csv
首先导入相关的包
from statsmodels.formula.api import logit
import pandas as pd
from scipy.stats.distributions import chi2
from math import log
import copy
import random
import itertools
import ast
import matplotlib.pyplot as plt
import seaborn as sns
pd.set_option('display.max_columns',None)这里我们主要用到的包有：statsmodels（用来建立logit模型）和scipy（用来做相关的统计检验）
读入数据:
df = pd.read_csv('model_data.csv')构造计算EL和MS模型LL值的函数
### 构建equally-likely模型的似然函数值
def ll_0(df):
N = len(df)
return -N*log(2)
### 构建market share模型的似然函数值
def ll_ms(df):
N_i = len(df[df['C3H17M']==0])
N_j = len(df[df['C3H17M']==1])
N = len(df)
y_in = N_i/N
y_jn = N_j/N
return N_i*log(N_i/N) + N_j*log(N_j/N)先计算一下EL和MS模型的LL值
print('EL模型的Loglikelihood值为：',round(ll_0(df),3))
print('MS模型的Loglikelihood值为：',round(ll_ms(df),3))
然后估计一下只有截距的logit模型
logit_intercept_only = logit('C3H17M~1',data=df).fit(method='bfgs',maxiter=100,disp=False)
print('intercept only模型的似然函数值为：'+str(round(logit_intercept_only.llf,3)))
#### 可以看出，实际上MS模型就是只包含截距的模型
可以发现，只有截距的模型的LL与MS模型的LL是一样的，因为只有截距项的模型，就是MS模型（那天在实验室虎哥还带我一起手推了一下
）
下面就开始估计只带有一个变量的logit模型，我们选的变量是FEMALE，即性别（假设女性可能更愿意居家办公？瞎猜的）
### 利用statsmodels构建只包含FAMALE的logit模型
logit_q1 = logit('C3H17M~FEMALE',data=df).fit(method='bfgs',maxiter=100,disp=False)
### 看一下结果
print(logit_q1.summary2())
可以看到，FEMALE的系数并不显著。同时，statsmodels还会给出模型的对应统计量，包括Pseudo R-squared（就是rho-squared），Log-Likelihood，LL-Null（只包含截距模型的Log-Likelihood），LLR p-value（该模型与只包含截距模型（也就是MS模型）的似然比检验值 。 这里，LLR p-value为0.146，大于0.05，可以认为加入FEMALE的模型相比只包含截距的模型并没有明显的提升。
这里我们再自己构造一个计算似然比的函数，与statsmodels估计的结果对比一下，看看是否一致：
### 构建似然比检验的函数
def likelihood_ratio_test(llmin, llmax,df):
lr = 2*(llmax-llmin)
p = chi2.sf(lr,df)
return round(p,4)### 与ms模型做似然比检验
likelihood_ratio_test(ll_ms(df),logit_q1.llf,1)
可以看到，我们算出来的似然比检验的p值与statsmodels给出的是一样的，欧耶。
再写一个计算rho-squared的函数
### 计算equally-kikely模型为基准的rho-squared
def rho_el_base(model_beta,df):
return round(1 - model_beta.llf/ll_0(df),4)
### 计算market share模型为基准的rho-squared
def rho_ms_base(model_beta,df):
return round(1 - model_beta.llf/ll_ms(df),4)计算一下
### 计算el-based rho-squared
rho_el_base(logit_q1,df)
### 计算ms-based rho-squared
rho_ms_base(logit_q1,df)
我们算的MS-based模型的rho-squared跟statsmodels给出的一样（都是0.005）,实际上还需要算一个adjusted rho-squared，这里懒了没有算（在rho-squared公式中加入自变量个数k进行penalty即可）
以上已经把第一问和第二问回答了，下面开始构建属于我们自己的二项Logit模型。
构建属于自己的Logit模型
这里我们的逻辑如下：
这里比较水，实际上是参照了 Mokhtarian P L, Salomon I. Modeling the desire to telecommute: The importance of attitudinal factors in behavioral models[J]. Transportation Research Part A: Policy and Practice, 1997, 31(1): 35-50. 的做法。下面是原版。
下面开始筛选变量：
### 构建所有自变量的dataframe，用来储存检验的结果
factor = pd.DataFrame(columns=['Variable','unique_values'])
for i in df.columns[1:]:
factor = factor.append([{'Variable':i,'unique_values':len(df[i].unique())}])
factor.reset_index(drop=True,inplace=True)
### 在数据集中增加一个用来计数的列
df['COUNT'] = 1
### 对每一个变量都进行检验
for i,j in enumerate(factor['Variable']):
#### 如果是二分类变量，进行卡方检验
if len(df[j].unique())==2:
factor.loc[i,'test'] = 'chi_square'
factor.loc[i,'p_value'] = stats.chi2_contingency(pd.pivot_table(df,index='C3H17M',columns='{}'.format(j),aggfunc='count',values='COUNT').fillna(0))[1]
### 如果是连续变量，进行t检验
else:
factor.loc[i,'test'] = 't_test'
### 首先通过levene检验确定是否方差齐
levene_pvalue = stats.levene(df[df['C3H17M']==0][j],df[df['C3H17M']==1][j])[1]
### p>0.05，方差齐，用方差齐的t检验
if levene_pvalue > 0.05:
factor.loc[i,'p_value'] = stats.ttest_ind(df[df['C3H17M']==0][j],df[df['C3H17M']==1][j],equal_var=True)[1]
### 否则用方差不齐的t检验
else:
factor.loc[i,'p_value'] = stats.ttest_ind(df[df['C3H17M']==0][j],df[df['C3H17M']==1][j],equal_var=False)[1]最终得到的DataFrame如下：
第一列是变量名称，第二列是变量独特值的数量（这里我们是比较懒的做法，如果一个变量只有两个值，我们就认为其是二元离散变量，否则认为其为连续变量，这么做的方法主要是出于以下考虑：如果一个变量以1，2，3，4，5的形式进入方程，并没有做dummy的处理的话，在方程估计时，实际上也是将其视作连续变量进行处理的，因此我们也就将其视作连续变量），第三列是采用的检验的方式（t检验或卡方检验），第四列为检验得到的p值。
筛选一下那些显著的变量：
### 筛选出显著的因素
final_factor = factor[factor['p_value']<0.05]观察变量间的相关系数，检验共线性：
### 计算皮尔逊相关系数，观察变量间共线性
corr_matrix = df[final_factor['Variable']].corr()
corr_matrix[corr_matrix>0.7]
可以看到，有不少变量之间的相关系数都大于0.7，我们将那些高度相关的变量去掉一个，避免多重共线性影响模型的假设检验准确性。
### 去除那些相关系数大于0.7的变量，得到最终的候选变量
final_factor = final_factor[~final_factor['Variable'].isin(['C1CONT','C2CONT','C2CONTC','CSO9FT9'])]
### 可以看到，所有候选变量间的相关系数都小于0.7，没有共线性问题
corr_matrix = df[final_factor['Variable']].corr()
corr_matrix[corr_matrix>0.7]
最终的候选变量，相关系数均小于0.7，可以进入最终的模型选择阶段。
现在，我们已经有了17个候选变量，下一步就是在这17个变量中寻找一个最佳的组合，确定最终的模型。这里我们用到的方法是，找到这17个变量的所有组合方式，也就是C17取1一直到C17取17（高中学过的排列组合），然后估计每一个组合对应的logit模型，比较每一个模型的AIC和BIC，分别选出AIC和BIC最小的模型作为我们的final model。
首先，构建一个列表，储存变量所有可能的组合方式：
### 构建所有候选变量的组合
variables = list(final_factor['Variable'])
variable_list = []
for i in range(1,19):
variable_list.extend(list(itertools.combinations(variables, i)))
可以看到，共有131071个组合。
然后开始估计对应的模型：
### 开始计算所有组合对应的模型，并保存各个模型的Log-likelihood, AIC, BIC
%%time
model_results = pd.DataFrame(columns = ['Variables','LL','AIC','BIC'])
for i in variable_list:
formula = 'C3H17M ~ {}'.format('+'.join(i))
model = logit(formula,data=df).fit(method='bfgs',maxiter=100,disp=False)
model_results = model_results.append([{'Variables':i,'LL':model.llf,'AIC':model.aic,'BIC':model.bic}])
print(i)
### 加一行变量个数信息
model_results['Number of variables'] = model_results.apply(lambda x: len(ast.literal_eval(x['Variables'])),axis=1)
### 将模型结果保存一下
model_results.to_csv('model_results_all.csv',index=0)
第一列是变量的组合，第二列是LL值，后面是AIC和BIC，最后一个是变量组合中一共有多少个变量。
下面可视化一下，不同变量个数组合的模型的goodness-of-fit
def box_plot(df,y,file_name):
fig,ax = plt.subplots(1,1,figsize=(10,5))
#ax.boxplot()
linedict = dict(color='k',linewidth=1.3)
sns.boxplot(x="Number of variables", y=y,
data=df,ax=ax,boxprops = dict(linestyle='-', linewidth=1.3, facecolor='white', edgecolor='k')
,width=0.5,
flierprops = dict(marker='o', markersize=2),
whiskerprops=linedict,
capprops=linedict,
medianprops=linedict )
#sns.swarmplot(x="Number of variables", y="BIC",
# data=results_all, color=".25",size=0.08,ax=ax)
sns.despine(offset=10, trim=True)
plt.xticks(fontsize=14)
plt.xlabel('Number of variables',fontsize=16)
plt.yticks(fontsize=14)
plt.ylabel(y,fontsize=16)
plt.savefig(file_name,dpi=300,bbox_inches='tight')box_plot(model_results,'AIC','AIC_box_plot')
box_plot(model_results,'BIC','BIC_box_plot')
这里我们只看boxplot的最低点就好，因为我们需要的是AIC和BIC最小的模型（AIC与BIC越小，证明模型越好（用最少的变量解释了最多的信息））。可以看到，以AIC作为评价标准的话，最优的模型出现在变量个数为9，10，11的时候（看不出来哪个最小，过于相近），而以BIC作为评价标准，最优的模型出现在变量个数为5的时候。也可以看出，其实变量个数过多和过少都不好，变量过少模型解释能力差，变量过多模型太复杂（这也正是AIC和BIC背后的思想）。
看一下AIC最优模型的变量
### AIC最优模型的变量
model_results[model_results['AIC']==model_results['AIC'].min()]['Variables'][106160]
估计一下对应的logit模型
print(logit('C3H17M ~ CSO9FT2+NOSPACE+F8AGEGR+UNAWARE+MANCONST+JOBCONST+C2CONTH+EFACT6+C1CONTC+CSO9FT1',data=df).fit(method='bfgs',maxiter=100).summary2())
BIC最优模型对应的变量
### BIC最优模型的变量
model_results[model_results['BIC']==model_results['BIC'].min()]['Variables'][6338]
估计对应的logit模型
print(logit('C3H17M ~ CSO9FT2+MANCONST+JOBCONST+EFACT6+CSO9FT1',data=df).fit(method='bfgs',maxiter=100).summary2())
文章到这里就结束啦，分别估计男女的模型的流程跟上面是一样的，就不放上来了，如果大家发现了错误还请后台私信我，也欢迎各位老板转发，点在看以及打赏，谢谢！
腾讯云开发者
