Python数据科学:线性回归诊断

2022-08-24 16:21:26 浏览数 (1)

接上篇的线性回归文章,传送门如下。

Python数据科学:线性回归

多元线性回归的前提条件:

  • 因变量不能和扰动项有线性关系
  • 自变量与因变量之间要有线性关系
  • 自变量之间不能有太强的线性关系
  • 扰动项或残差独立且应服从均值为0、方差一定的正态分布

/ 01 / 残差分析

残差分析是线性回归诊断的重要环节。

残差应服从的前提条件有三个:

  • 残差方差齐性
  • 残差独立同分布
  • 残差不能和自变量相关(不能检验)

通过查看残差图来查看残差情况。

残差图可分为四类:

  • 残差正常分布:残差随机分布,上下界基本对称,无明显自相关,方差基本齐性
  • 残差曲线分布:残差与预测值呈曲线关系,说明自变量与因变量不是线性关系
  • 残差方差不齐:残差上下界基本对称,但随着预测值的增大,其上下幅度也会不断增大
  • 残差周期变化:残差随预测值增大而呈周期性变化,说明自变量与因变量可能是周期性变化

下面以之前线性回归文章里的模型为例。

代码语言:javascript复制
# 简单线性回归模型,平均支出和收入
ana1 = lm_s
# 训练数据集的预测值
exp['Pred'] = ana1.predict(exp)
# 训练数据集的残差
exp['resid'] = ana1.resid
# 绘制收入与残差的散点图
exp.plot('Income', 'resid', kind='scatter')
plt.show()

得到模型的残差情况,随着预测值增大,残差基本保持上下对称。

但残差正负的幅度有逐渐变大的趋势,即模型有方差不齐的问题。

异方差的处理方法,可以把数据取对数,所以这里把平均支出数据做对数处理。

代码语言:javascript复制
# 使用简单线性回归建立模型,平均支出对数数据
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()

得到模型的残差情况,发现残差正负的幅度有改善的趋势。

书中是说异方差现象得到了改善,但我感觉已经挺好的了...

上面只是针对平均支出数据取对数,下面对收入数据也取对数,使得二者增加的百分比值大致一样。

代码语言:javascript复制
# 对收入数据取对数,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²情况。

代码语言:javascript复制
# 获取三种模型的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)是指标准化后的残差。

代码语言:javascript复制
# 学生化残差计算
exp['resid_t'] = (exp['resid'] - exp['resid'].mean()) / exp['resid'].std()
# 样本量为几百时取2,样本量为上千时取3
print(exp[abs(exp['resid_t']) > 2])

输出结果如下。

共两个,从之前的散点图,大家应该也能知道是哪两个点了。

以去除强影响点的数据建立回归模型。

代码语言:javascript复制
# 使用简单线性回归建立模型,去除强影响点
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 / 多重共线性分析

自变量之间不能有强共线性,又称多重共线性。

本次使用方差膨胀因子去诊断及减轻多重共线性。

在之前的数据加入当地房屋均价、当地平均收入数据。

代码语言:javascript复制
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,说明存在多重共线性。

按道理此时应该删除其中一个变量的。

这里使用高出平均收入的比例代替收入数据列,能够较好的体现出信息。

代码语言:javascript复制
# 高出平均收入的比例
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一波,以后慢慢消化...

0 人点赞