用Python底层编写进行计量经济分析（二）：多重共线性（原因、结果、检验：条件数/方差膨胀因子、补救：岭回归）

偏执的太偏执、 2023-07-22 09:09 12阅读 0赞

系列文章：

1.  [多元线性回归和显著性检验（参数估计、T检验、F检验、拟合优度）][T_F]
2.  [多重共线性（导致结果、检验——方差膨胀因子、补救措施——岭回归）][Link 1]
3.  [异方差（导致结果、检验——White、补救措施——广义线性回归）][White]
4.  [自相关（导致结果、检验——D-W、补救措施——广义线性回归）][D-W]

写的比较仓促，代码中如有错误欢迎指正！

上一次对[多元线性回归的估计以及参数和方程的显著性进行了python实现][T_F]。但是这些都是建立多元线性回归的几个假设基础之上的：

1.  模型符合线性模式
2.   X X X满秩（无多重共线）
3.  零均值价值： E ( ε i ∣ X i ) = 0 E(ε\_i|X\_i)=0 E(εi∣Xi)=0（自变量外生）
4.  同方差： V a r ( ε i ∣ X i ) = σ Var(ε\_i|X\_i)=σ Var(εi∣Xi)=σ
5.  无自相关： c o v ( ε i , ε j ) = 0 cov(ε\_i, ε\_j)=0 cov(εi,εj)=0
6.  球形扰动：ε\_i是正态分布

若果模型违反了相应的假设就会犯对应的错误，我们在计量经济分析中的检验就是检验出是否可能犯了某一类错误，若果极有可能犯了一种错误时，我们应该怎么修正它，才能保证分析的结果是有效的。

#### 一、多重共线性 ####

**先来考虑一个现实问题，我们在预测一个活动的KPI时，利用多个变量作为解释变量，历史KPI做被解释变量，如果解释变量中同时包含PV和UV，这样的回归是好还是坏？**

**1.1.完全多重共线性与不完全多重共线性**

上面对于多元线性回归的假设中有一条是要求 X X X满秩，即： R a n k ( X ) = k Rank(X)=k Rank(X)=k，当X中存在一个变量完全可以由另一个变量表示的时候，变量间就出现了**完全多重共线**，即：

λ 1 x 1 + λ 2 x 2 + . . . + λ k x k = 0 λ\_1x\_1+λ\_2x\_2+...+λ\_kx\_k=0 λ1x1\+λ2x2\+...\+λkxk=0

由秩的定义可知，此时 X X X不再满秩，即 R a n k ( X ) < k Rank(X)<k Rank(X)<k，这违反了多元线性回归的基本假定。

实际上当 x i x\_i xi与 x j x\_j xj之间存在很想的线性相关性时，对模型的估计（区间估计）和检验也是非常不利的。出现这种情况，我们称为**不完全多重共线性**

**1.2.多重共线性的原因**  
由于经济和现实的生产过程中，许多变量往往是强相关和有相同变化趋势的。这就导致了多重共线性有事在所难免。实际上，**多重共线性大多不是模型问题，而是数据问题**。模型中我们大量的采用滞后变量或者变量选取不当，也会导致共线性的产生。

**1.3.完全多重共线性的后果**  
当存在完全多重共线性时，多元线性回归可能存在一下后果：

1.  模型参数无法估计  
    比如 x 1 − 2 x 2 = 0 x\_1-2x\_2=0 x1−2x2=0，如果 y = x 1 + 2 x 2 y=x\_1+2x\_2 y=x1\+2x2，则y可以表示成 y = 3 x 1 − 2 x 2 y=3x\_1-2x\_2 y=3x1−2x2，此时参数是无法估计的
2.  参数估计的方差无限大  
    由于不满秩，导致参数的方差估计中的分母为0，方差无限大（参数的方差估计详见前一节）

**1.4.不完全多重共线性的后果**

1.  参数估计方差变大  
    这与上面的完全多重共线性的参数估计方差无限大类似。
2.  t检验失效  
    由于参数估计方差，我们构造的T统计量会变小（T统计量的分母就是参数估计方差除自由度开根号）。此时我们做的t检验，无法帮助我们判断参数是否显著。  
    （不完全多重共线性时，参数的估计仍然是无偏的）

**1.4.多重共线性的检验**

1.  简单相关系数矩阵法、  
    我们可以计算 X X X中每两个变量组合之间的相关系数，来确定变变量间是否存在很强的线性相关性，利用dataframe的`df.corr()`计算即可完成。
2.  条件数  
    计算 ( X ′ X ) (X'X) (X′X)的最大特征根与最小特征根比值的平方根，当做条件数：  
    条件数= m a x ( λ i ) / m i n ( λ i ) \\sqrt\{max(λ\_i)/min(λ\_i)\} max(λi)/min(λi)  
    如果条件数超过20，即认为存在比较严重的多重共线性。（回想一下PCA的原理，和特征值 λ i λ\_i λi的含义，其实条件数还是很直观的）  
    **在下面实际的python测试中，条件数达到了4573个！但是根据其他法则判断，共线性又不是十分的明显，个人感觉最好不要单一的用条件数去判断多重共线性。**
3.  方差膨胀因子  
    偷懒一下，直接粘贴网上的定义吧  
    ![在这里插入图片描述][20200405115256876.png]  
    一般 V I F i VIF\_i VIFi大于10，我们即认为变量 x i x\_i xi存在多重共线性。 V I F VIF VIF的优点是可以定位到存在共线性的变量，缺点是计算比较复杂。我们可以计算出方程的条件数，判断是否存在多重共线性，然后再利用X的相关系数矩阵定位存在多重共线性的变量（计算相关系数矩阵也是比较耗时的）。

**1.5.多重共线性的补救**

1.  增加样本量  
    增加样本量可以达到减小参数估计方差的目的，在存在共线性时提高估计的可靠性
2.  逐步回归法  
    每种变量组合都进行一次回归，选出最好的组合，这种做法实在是太麻烦和耗时了
3.  PCA  
    我们可以对变量数据 X X X做一次PCA，得到新的矩阵 Z Z Z，Z中的主成分时正交的，完全没有共线性。但是这样缺点也很明显，回归是没有经济意义的，我偏离了我们做计量经济分析的初衷！
4.  岭回归  
    在机器学习领域我们知道有两种添加正则化的回归方式：岭回归和Lasso回归两种方式，对应着L2正则和L1正则。岭回归让参数变得很小，Lasso回归让变量变得系数，可用于特征选择。。。。  
    **在解决多重共线性中的岭回归办法和机器学习中添加L1正则的岭回归其实就是一个东西！**。岭回归是最小二乘的一种改进方法，用岭回归时参数的估计变为：  
     β = ( X ′ X + k I ) − 1 X ′ y β=(X'X+kI)^\{-1\}X'y β=(X′X\+kI)−1X′y  
    其中 k k k为岭回归参数（在机器学习中回应惩罚的力度）**，需要注意的是，岭回归的参数估计是有偏的，参数 β β β的方差估计也需要跟着一起变**。此时 β β β的方差估计为：  
    $ V a r ( β ∣ k ) = σ 2 ( X ′ X + k I ) − 1 ∗ ( X ′ X ) ∗ ( X ′ X + k I ) Var(β|k) = σ^2(X'X+kI)^\{-1\}\*(X'X)\*(X'X+kI) Var(β∣k)=σ2(X′X\+kI)−1∗(X′X)∗(X′X\+kI) $  
    对应的t检验的统计量也会跟着一起改变。

#### 二、多重共线性检验和补救的python实现 ####

我们任沿用上一节的多元线性回归代码，再其基础上进行添加，**我们需要实现的有：计算被解释变量的相关系数矩阵、没个变量的反差膨胀因子VIF，方程拟合后的条件数、对方程惊醒岭回归**。先把框架定义好：

class mul_linear_model(object):
      #用于计算数据的相关系数矩阵
      def data_relation(self):
          pass
      #用于计算每个变量的VIF
      def VIF_cul(self, ind):
          pass
      #所有变量的分行差膨胀因子检验
      def VIF_test(self):
          pass
      #计算拟合的条件数
      def  condition_num(self):
      	  pass
      #岭回归
      def Ridge_fit(self, k = 0.1):
      	  pass

下面就是完整的多重共线性的检验和补救方法。

# -*- coding: utf-8 -*-
    """
    Created on Fri Apr  3 15:40:38 2020
    
    @author: nbszg
    """
    
    import numpy as np
    import pandas as pd
    import scipy.stats as st
    
    class mul_linear_model(object):
      #intercept定义是否有截距项，在计量经济分析中，基本都要有截距项，否则计算出的拟合优度是没有意义的
      def __init__(self, y, X, intercept=True):
          ...
      #拟合模型
      def fit(self):
         ...
          return self.b
      #预测函数
      def predict(self, X_p):
          ...
      #参数T检验
      def T_test(self, alpha):
         ...
      #拟合优度计算
      def R_square(self):
         ...
      #整体F检验
      def F_test(self, alpha):
          ...
      #用于计算数据的相关系数矩阵
      def data_relation(self):
          #需要区分有截距项和没有截距项的情况
          if self.intercept:
              return self.data_x.loc[:,self.data_x.columns!='intercept'].corr()
          else:
              return self.data_x.corr()
      #用于计算每个变量的VIF    
      def VIF_cul(self, ind):
          #使用最小二乘法(X'X)的逆成X'y
          yk = np.array(self.X[:, ind].T)
          #有截距项时：
          if self.intercept:
              Xk = self.X[:, [i for i in range(1,ind)]+[j for j in range(ind+1, self.K)]]
          #无截距项时：
          else:
              Xk = self.X[:, [i for i in range(ind)]+[j for j in range(ind+1, self.K)]]
          
          #直接调用自身的类方法记性计算其他Xj对Xk回归的可决系数
          vif_model = mul_linear_model(yk, pd.DataFrame(Xk), intercept=False)
          vif_model.fit(output=False)
          Rk_2 = vif_model.R_square()
          #计算方差膨胀因子
          VIF_k = 1/(1-Rk_2)
          return VIF_k
          
      #所有变量的分行差膨胀因子检验
      def VIF_test(self):
          #循环计算每个变量的方差膨胀因子
          for i in range(1, self.K):
              VIF_k = self.VIF_cul(i)
              print("variable {0}'s VIF estimate is : {1}".format(self.columns[i], VIF_k))
      
      #计算拟合的条件数
      def condition_num(self):
          #有截距项时
          if self.intercept:
              M = np.dot(self.XT[[i for i in range(1, self.K)],:], self.X[:, [j for j in range(1, self.K)]])
          else:
              M = np.dot(self.XT, self.X)
          #计算X'X的特征值与特征向量
          eigenvalue, eigenvector = np.linalg.eig(M)
          #计算条件数
          cond_num = np.sqrt(eigenvalue[0]/eigenvalue[-1])
          print("data X's condition number is: {0}".format(cond_num))
          
      
      #岭回归
      def Ridge_fit(self, k = 0.1, output=True):
          #使用岭回归(X'X+kI)的逆成X'y
          self.b = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K))),self.XT),self.y).T)
          
          if output:
              for i in range(self.K):
                  print("variable {0}'s cofe estimate is: {1}".format(self.columns[i], self.b[i][0]))
          
          #计算t检验的因子需要跟着一起变动！！！
          self.S = np.array(np.dot(np.dot(np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K))), np.dot(self.XT,self.X)), np.linalg.inv(np.dot(self.XT,self.X)+np.mat(np.eye(self.K)))))
          
          return self.b

利用[上一节的数据和例子][T_F]继续进行检验：

print('data columns related coefficient is: \n{0}'.format(linear_model.data_relation()))
    print('\n')
    linear_model.VIF_test()
    print('\n')
    linear_model.condition_num()
    print('\n')
    linear_model.Ridge_fit()

得到的检验结果和岭回归系数如下，可以看到条件数达到了4573个！但是根据其他法则判断，共线性又不是十分的明显，个人感觉最好不要单一的用条件数去判断多重共线性。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3UwMTE1MTcxMzI_size_16_color_FFFFFF_t_70]

[T_F]: https://blog.csdn.net/u011517132/article/details/105286113
[Link 1]: https://blog.csdn.net/u011517132/article/details/105323135
[White]: https://blog.csdn.net/u011517132/article/details/105392242
[D-W]: https://blog.csdn.net/u011517132/article/details/105581371
[20200405115256876.png]: /images/20230528/3e687eaa25f64018a4543fde3ba04518.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3UwMTE1MTcxMzI_size_16_color_FFFFFF_t_70]: /images/20230528/996bf06d9fe14bdcaa8d5d3b030665f2.png