如何用Python检验线性回归的假设是否满足

文章目录

• 一、导入数据
• 二、进行线性回归
• 三、假设检验
  • 假设1：线性假设（Linearity of the model）
  • 假设2：误差项（errors）的期望均值为0
  • 假设3：无（完全）多重共线性假设
  • 假设4：异方差检验
• 总结

提示：该系列文章侧重代码编程实践，不涉及理论基础。如代码无法正常运行或有任何疑问，请告诉我，我将设法解答。

线性回归作为数据分析的基础方法之一被广泛应用，但回归结果是否能够满足有效（effective）、无偏（unbiased）和一致性（consistent）呢？这篇文章将逐步详细介绍如何用Python（Jupyter Notebook）来检验我们的多元线性回归是否满足五大假设的要求。

   点此下载文章数据    提取码：0001
   若超链接不可用，请复制以下地址：
链接：https://pan.baidu.com/s/1Mv0TJ4GL2dehytFmbfzmog
提取码：0001

   做好准备，我们现在开始🚀

一、导入数据

   本文以受教育年限对个人工资这一经典研究问题为例进行讲解。为保证示例简明，线性回归仅适用了数据中的部分连续数据变量。

#导入数据 请预先下载好数据并上传至Jupyter Notebook
data=pd.read_csv('train.csv')
#检查一下数据吧，显示数据里的前4行
train.head()

    运行后，以看到数据如下
 

   接着我们定义定义自变量和因变量

#定义自变量集X包含educ（受教育年限）、hours（工作小时数）、IQ（智商得分）以及exper（工作总年限）四个变量
X = train[['educ','hours', 'IQ','exper']]
#看一下我们的自变量吧
X.head()
#定义因变量y是lwage（小时工资的对数）
y = pd.Series(train.lwage)

     我们可以看到在17个变量里，我们只选用了其中的四个作为自变量
 

二、进行线性回归

   Python里线性回归的方法不止一种，我们用以下方式进行：

import statsmodels.api as sm
from statsmodels.formula.api import ols

lin_reg = ols('y~educ+hours+IQ+exper', data=train).fit()
lin_reg.summary()

我们得到如下结果：

   对于结果的解释，此处不是重点。我们现在就可以开始对于假设条件是否满足的检验了👇

三、假设检验

   在你们的计量课程中，应该已经学过BLUE。什么意思呢？也就是说根据高斯-马尔可夫（Gauss-Markov ）定理，在线性回归模型中，普通最小二乘 (OLS)估计 是最佳线性无偏估计 (BLUE)，但前提是满足经典线性回归假设：

1. 所有参数均为常数且满足线性关系
2. 误差项（errors）的期望均值为0，且与自变量无关
3. 因变量不能有完全共线性
4. 误差项是同方差（homoscedasticity of errors），也即误差项的方差不受到自变量的影响，是一个固定不变的值

   此外，与其他无偏线性估计量相比，BLUE 中的“最佳”意味着估计量的方差最小。

假设1：线性假设（Linearity of the model）

   在线性回归中因变量y是自变量X的线性方程，因此我们的参数（parameters）必须是线性的。如果我们用线性模型去拟合非线性的数据，那么就会造成很严重的预测误差。

   那么，我们可以用两张图来检查是否符合线性：（1）实际值与预测值；（2）残差与预测值。如果在第一张图中，所有的点都围绕对角线对称分布；或者第二张途中所有的点都成水平线对称分布，那么我们认为是符合线性假设的。
为了检测非线性，可以检查观察值与预测值或残差与预测值的图。 期望的结果是点在前一个图中围绕对角线对称分布，在后一个图中围绕水平线对称分布。 在这两种情况下，方差大致恒定。

%matplotlib inline
%config InlineBackend.figure_format ='retina'
import seaborn as sns 
import matplotlib.pyplot as plt
import statsmodels.stats.api as sms
sns.set_style('darkgrid')
sns.mpl.rcParams['figure.figsize'] = (15.0, 9.0)
#在这里定义了一个名为linearity_test的函数用来检验线性假设，如果做了不止一线性回归，只要运行这部分代码一次，后面直接可以调用，不必重复写代码
def linearity_test(model, y):
    '''
    Function for visually inspecting the assumption of linearity in a linear regression model.
    It plots observed vs. predicted values and residuals vs. predicted values.
    
    Args:
    * model - fitted OLS model from statsmodels
    * y - observed values
    '''
    fitted_vals = model.predict()
    resids = model.resid

    fig, ax = plt.subplots(1,2)
    
    sns.regplot(x=fitted_vals, y=y, lowess=True, ax=ax[0], line_kws={'color': 'red'})
    ax[0].set_title('Observed vs. Predicted Values', fontsize=16)
    ax[0].set(xlabel='Predicted', ylabel='Observed')

    sns.regplot(x=fitted_vals, y=resids, lowess=True, ax=ax[1], line_kws={'color': 'red'})
    ax[1].set_title('Residuals vs. Predicted Values', fontsize=16)
    ax[1].set(xlabel='Predicted', ylabel='Residuals')

##到这里我们对 linearity_test 这个函数的设定就完成啦，以后每次只需要输入下面这一行，把y前面替换成你所回归的模型的名称就大功告成！
linearity_test(lin_reg, y) 

   我们发现在本示例中，所有点都大致围绕对角线均匀分布，在残差-预测值图中，所有点大致均匀分布在水平线两侧，因此，我们认为大致符合线性假设。

假设2：误差项（errors）的期望均值为0

   这个检验异常轻松，代码如下：

#是的，就这一句！
lin_reg.resid.mean()

我们得到的结果：

-4.901604968612456e-15

残差的均值非常的小近乎为0。

假设3：无（完全）多重共线性假设

   如果模型中自变量之间存在精确相关关系或者高度相关，那么模型就会难以准确估计。比如，有X1和X2两个自变量，如果X1 = 2 + 3 * X2，那么X1和X2就是精确相关，从而违反假定。

   当然了，一定程度的多重共线性是可以存在于模型中，只要它不是“完美的”。 在前一种情况下，估计效率较低，估计将不那么精确并且对特定数据集高度敏感，但仍然无偏。

   我们可以使用方差膨胀因子 (VIF) 检测多重共线性。 如果变量完全没有线性关系， VIF 的值将为 1。也即VIF值越接近于1，多重共线性越轻，反之越重。

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame()
vif_data["feature"] = X.columns

vif_data["VIF"] = [variance_inflation_factor(X.values, i)
                          for i in range(len(X.columns))]
  
print(vif_data)

 结果如下：
 

   通常而言，我们以VIF=10（5也很常用）为判断标准，如果大于10则认为存在多重共线性，需要通过迭代删除VIF值高的变量。
   小贴士：可以通过相关系数矩阵去识别相关性。

假设4：异方差检验

   当残差呈现异方差时，会很难确定真实标准差，通常会导致置信区间过宽/过窄。例如，如果残差方差随时间而增大， 那么对于非样本数据预测的置信区间会窄得超乎寻常。 异方差的另一个影响也可能是在估计系数时对误差方差最大的数据赋予过多的权重。

   为了研究残差是否同方差，我们可以查看残差（或者标准差）与拟合值的关系图。 当残差随拟合值增长时，就可能存在异方差。此外，我们还可以采用Breusch-Pagan和Goldfeld-Quandt来进行检验。这两个检验的原假设都是同方差。

%matplotlib inline
%config InlineBackend.figure_format ='retina'


import statsmodels.stats.api as sms
import numpy as np

sns.set_style('darkgrid')
sns.mpl.rcParams['figure.figsize'] = (15.0, 9.0)

#类似的，我们从这里开始定义homoscedasticity_test这个函数用来检验异方差
def homoscedasticity_test(model):
    '''
    Function for testing the homoscedasticity of residuals in a linear regression model.
    It plots residuals and standardized residuals vs. fitted values and runs Breusch-Pagan and Goldfeld-Quandt tests.
    
    Args:
    * model - fitted OLS model from statsmodels
    '''
    fitted_vals = model.predict()
    resids = model.resid
    resids_standardized = model.get_influence().resid_studentized_internal

    fig, ax = plt.subplots(1,2)

    sns.regplot(x=fitted_vals, y=resids, lowess=True, ax=ax[0], line_kws={'color': 'red'})
    ax[0].set_title('Residuals vs Fitted', fontsize=16)
    ax[0].set(xlabel='Fitted Values', ylabel='Residuals')

    sns.regplot(x=fitted_vals, y=np.sqrt(np.abs(resids_standardized)), lowess=True, ax=ax[1], line_kws={'color': 'red'})
    ax[1].set_title('Scale-Location', fontsize=16)
    ax[1].set(xlabel='Fitted Values', ylabel='sqrt(abs(Residuals))')

    bp_test = pd.DataFrame(sms.het_breuschpagan(resids, model.model.exog), 
                           columns=['value'],
                           index=['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value'])

    gq_test = pd.DataFrame(sms.het_goldfeldquandt(resids, model.model.exog)[:-1],
                           columns=['value'],
                           index=['F statistic', 'p-value'])

    print('\n Breusch-Pagan test ----')
    print(bp_test)
    print('\n Goldfeld-Quandt test ----')
    print(gq_test)
    print('\n Residuals plots ----')

#运行一次以后，直接调用就可以再次检验异方差啦
homoscedasticity_test(lin_reg)

结果如下：
 

   如果是同方差性，那么图中的点应该是随机没有特定趋势的，也就是说残差值不会随着拟合值增加或减少。 我们可以看到，我们的数据集并不完全如此，有一些波动，而且Breusch-Pagan检验也在1%显著性上拒绝原假设，所以该数据很可能存在异方差。

总结

   以上我们对线性假设、误差期望值为0、无多重共线性以及同方差的假设进行了检验。 如有任何问题，请跟我联系哦📧