Python数据科学：线性回归诊断
Python数据科学：线性回归诊断
接上篇的线性回归文章，传送门如下。
多元线性回归的前提条件：
- 因变量不能和扰动项有线性关系
- 自变量与因变量之间要有线性关系
- 自变量之间不能有太强的线性关系
- 扰动项或残差独立且应服从均值为0、方差一定的正态分布
/ 01 / 残差分析
残差分析是线性回归诊断的重要环节。
残差应服从的前提条件有三个：
- 残差方差齐性
- 残差独立同分布
- 残差不能和自变量相关(不能检验)
通过查看残差图来查看残差情况。
残差图可分为四类：
- 残差正常分布：残差随机分布，上下界基本对称，无明显自相关，方差基本齐性
- 残差曲线分布：残差与预测值呈曲线关系，说明自变量与因变量不是线性关系
- 残差方差不齐：残差上下界基本对称，但随着预测值的增大，其上下幅度也会不断增大
- 残差周期变化：残差随预测值增大而呈周期性变化，说明自变量与因变量可能是周期性变化
下面以之前线性回归文章里的模型为例。
# 简单线性回归模型,平均支出和收入
ana1 = lm_s
# 训练数据集的预测值
exp['Pred'] = ana1.predict(exp)
# 训练数据集的残差
exp['resid'] = ana1.resid
# 绘制收入与残差的散点图
exp.plot('Income', 'resid', kind='scatter')
plt.show()得到模型的残差情况，随着预测值增大，残差基本保持上下对称。
但残差正负的幅度有逐渐变大的趋势，即模型有方差不齐的问题。
异方差的处理方法，可以把数据取对数，所以这里把平均支出数据做对数处理。
# 使用简单线性回归建立模型,平均支出对数数据
ana2 = ols('avg_exp_ln ~ Income', data=exp).fit()
exp['Pred'] = ana2.predict(exp)
# 训练数据集的残差
exp['resid'] = ana2.resid
# 绘制收入与残差的散点图
exp.plot('Income', 'resid', kind='scatter')
plt.show()得到模型的残差情况，发现残差正负的幅度有改善的趋势。
书中是说异方差现象得到了改善，但我感觉已经挺好的了...
上面只是针对平均支出数据取对数，下面对收入数据也取对数，使得二者增加的百分比值大致一样。
# 对收入数据取对数,np.log(e)=1
exp['Income_ln'] = np.log(exp['Income'])
# 使用简单线性回归建立模型,平均支出对数数据
ana3 = ols('avg_exp_ln ~ Income_ln', data=exp).fit()
exp['Pred'] = ana3.predict(exp)
# 训练数据集的残差
exp['resid'] = ana3.resid
# 绘制收入与残差的散点图
exp.plot('Income_ln', 'resid', kind='scatter')
plt.show()书中是说异方差现象消除了，真的是没看出来和上一张图有什么大区别...
对比一下三种模型的R²情况。
# 获取三种模型的R²值
r_sq = {'exp~Income': ana1.rsquared, 'ln(exp)~Income': ana2.rsquared, 'ln(exp)~ln(Income)': ana3.rsquared}
print(r_sq)
# 输出结果
{'ln(exp)~Income': 0.4030855555329649, 'ln(exp)~ln(Income)': 0.480392799389311, 'exp~Income': 0.45429062315565294}发现都取对数的模型R²最大，模型算是三种里较好的。
残差的自相关关系判断可以使用DW检验(Durbin-Watson检验)。
当DW值趋近2时，可以认为残差无自相关关系。
下面是以都取对数的模型输出的判断指标。
发现都取对数的模型，其DW值为1.368。
这我也不知到底是什么相关了，可能是正相关吧...
/ 02/ 强影响点
当某个点离群太远时，拟合的回归线会受到这个点的强烈干扰，从而改变回归线的位置。
这便是强影响点。
这里我们可以使用预测值-学生化残差图来识别强影响点。
学生化残差(SR)是指标准化后的残差。
# 学生化残差计算
exp['resid_t'] = (exp['resid'] - exp['resid'].mean()) / exp['resid'].std()
# 样本量为几百时取2,样本量为上千时取3
print(exp[abs(exp['resid_t']) > 2])输出结果如下。
共两个，从之前的散点图，大家应该也能知道是哪两个点了。
以去除强影响点的数据建立回归模型。
# 使用简单线性回归建立模型,去除强影响点
exp2 = exp[abs(exp['resid_t']) <= 2].copy()
ana4 = ols('avg_exp_ln ~ Income_ln', data=exp2).fit()
exp2['Pred'] = ana3.predict(exp)
# 训练数据集的残差
exp2['resid'] = ana3.resid
# 绘制收入与残差的散点图
exp2.plot('Income_ln', 'resid', kind='scatter')
plt.show()得到结果。
这个残差结果还是不错的。
/ 03 / 多重共线性分析
自变量之间不能有强共线性，又称多重共线性。
本次使用方差膨胀因子去诊断及减轻多重共线性。
在之前的数据加入当地房屋均价、当地平均收入数据。
exp2['dist_home_val_ln'] = np.log(exp['dist_home_val'])
exp2['dist_avg_income_ln'] = np.log(exp['dist_avg_income'])
def vif(df, col_i):
# 获取变量
cols = list(df.columns)
# 去除因变量
cols.remove(col_i)
# 获取自变量
cols_noti = cols
# 多元线性回归模型建立及获取模型R²
formula = col_i + '~' + '+'.join(cols_noti)
r2 = ols(formula, df).fit().rsquared
# 计算方差膨胀因子
return 1. / (1. - r2)
# 获取自变量数据
exog = exp2[['Age', 'Income_ln', 'dist_home_val_ln', 'dist_avg_income_ln']]
# 遍历自变量,获取其VIF值
for i in exog.columns:
print(i, '\t', vif(df=exog, col_i=i))获取的结果。
发现收入和当地平均收入的方差膨胀因子大于10，说明存在多重共线性。
按道理此时应该删除其中一个变量的。
这里使用高出平均收入的比例代替收入数据列，能够较好的体现出信息。
# 高出平均收入的比例
exp2['high_avg_ratio'] = exp['high_avg'] / exp2['dist_avg_income']
# 获取自变量数据
exog2 = exp2[['Age', 'high_avg_ratio', 'dist_home_val_ln', 'dist_avg_income_ln']]
# 遍历自变量,获取其VIF值
for i in exog2.columns:
print(i, '\t', vif(df=exog2, col_i=i))输出结果。
发现各变量的方差膨胀因子均较小，说明不存在共线性。
当然上述方法只能减轻共线性对模型的干扰，并不能完全消除多重共线性。
/ 04 / 总结
建立一个合理的线性回归模型的步骤如下。
初始分析：确定研究目标，收集数据。
变量选择：找到对因变量有影响力的自变量。
验证模型假定：
- 设置模型，选择回归方法，选择变量，以及变量以何种形式放入模型
- 解释变量和扰动项不能相关
- 解释变量之间不能有强线性关系
- 扰动项独立同分布
- 扰动项服从正态分布
多重共线性与强影响点的诊断与分析：修正回归模型。
预测和解释：使用模型来预测与解释。
说实话，这部分内容都不太好理解。
大多数我都选择先Mark一波，以后慢慢消化...