线性回归中的假设检验

0 引言

本文介绍一元线性回归的基本假设，推导方法和统计检验，然后介绍Python编程实现，最后结合Python中statsmodels包对比计算结果。

1 一元线性回归模型

对于n个自变量，$(x^1,x^2,x^3,\dots ,x^n)$，因变量 $y^i=a\ast x^i+b+{ϵ}^i,i=1,2,\dots n$，其中${ϵ}^i\sim N(0,{\sigma }^2)$是服从正态分布的随机变量，且彼此独立。一元线性回归的目标就是对模型参数a和b进行估计。这里需要注意$x^i$不是随机变量。

利用最小二乘法对模型进行求解，首先定义代价函数：
$$J=\sum _{i=1}^n(y^i-(a\ast x^i+b){)}^2$$
将代价函数分别对a和b进行求导，并令导数为0，得到：
$$\begin{array}{rl}\frac{\partial J}{\partial a}& =-\sum _{i=1}^{n}2\ast (y^i-(a\ast x^i+b))\ast x^i=0\\ \frac{\partial J}{\partial b}& =-\sum _{i=1}^{n}2\ast (y^i-(a\ast x^i+b))=0\end{array}$$
求解得到最小而成解
$$\begin{array}{rl}\hat{a}& =\frac{{\sum }_{i=1}^n(x^i-\bar{x})\ast y^i}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}\\ \hat{b}& =\bar{y}-\hat{a}\ast \bar{x}\end{array}$$
其中，
$$\begin{array}{rl}\bar{x}& =\sum _{i=1}^n{x}^i\\ \bar{y}& =\sum _{i=1}^n{y}^i\end{array}$$

2 对于回归方程的检验

在一元线性回归中（多元也一样），假设检验主要分为两类，F检验和T检验，这里分别进行介绍。

F检验

构建下面的统计量：
$$\begin{array}{rl}SST& =\sum _{i=1}^n(y^i-\bar{y}{)}^2\\ SSR& =\sum _{i=1}^n(\widehat{y^i}-\bar{y}{)}^2\\ SSE& =\sum _{i=1}^n(\widehat{y^i}-y^i{)}^2\end{array}$$

其中$\widehat{y^i}=\hat{a}\ast x^i+\hat{b}$，为模型的估计值。易证$SST=SSR+SSE$（证明过程比较繁琐，但是思路比较简单，这里不再赘述）。对于空假设：
$$\begin{gathered} H_0:a=0 \\ {H}_1:a\ne 0 \end{gathered}$$
在空假设$H_0$成立的情况下，$SST/{\sigma }^2$服从自由度为 $n-1$ 的卡方分布，这一点在一般的概率统计中都会有介绍，不再赘述。而对于$SSE/{\sigma }^2$的自由度为$n-2$的证明比较复杂，可以根据经验，如果模型中需要估计p个参数，那么$SSE/{\sigma }^2$的自由度就为$n-p$，这里因为要估计a和b两个参数，因而自由度为n-2。因为$SSR=SST-SSE$，从而$SSR$服从自由度为$1$的卡方分布。且SSR和SSE相互独立，从而有
$$F=\frac{SSR/1}{SSE/(n-2)}\sim F(1,n-2)$$
构成F统计量。在空假设成立的情况下，估计误差SSE为比较大，因而统计量F会比较小。当F取得较大值$F>F_{1-\alpha }$是，拒绝原假设。

T检验

T检验是直接对回归系数的检验。已知
$$\hat{a}=\frac{{\sum }_{i=1}^n(x^i-\bar{x})\ast y^i}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}$$
$\hat{a}$是$y=[y^1,y^2,\dots ,y^i{]}^T$的线性组合，因而$\hat{a}$也服从高斯分布。又因为：
$$\begin{array}{rl}E(\hat{a})& =\frac{{\sum }_{i=1}^n(x^i-\bar{x})\ast E(y^i)}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}\\ & =\frac{{\sum }_{i=1}^n(x^i-\bar{x})\ast E(a\ast x^i+b+{ϵ}^i)}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}\\ & =\frac{{\sum }_{i=1}^n(x^i-\bar{x})\ast (a\ast x^i+b)}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}\\ & =\frac{{\sum }_{i=1}^{n}a(x^i-\bar{x})\ast (x^i-\bar{x})}{{\sum }_{i=1}^n(x^i-\bar{x}{)}^2}\\ & =a\end{array}$$
为了计算方差，我们将上述表达式定义为向量形式，首先定义向量$x_{\mu }=[x^1-\bar{x},x^2-barx,\dots ,x^n-\bar{x}]$，则：
$$\hat{a}=x_{\mu }^{T}y/x_{\mu }^T{x}_{\mu }$$
则估计参数的方差为:
$$\begin{array}{rl}Var(\hat{a})& =Cov(x_{\mu }^{T}y/x_{\mu }^T{x}_{\mu },x_{\mu }^{T}y/x_{\mu }^T{x}_{\mu })\\ & =Cov(x_{\mu }^{T}y,x_{\mu }^{T}y)/(x_{\mu }^T{x}_{\mu }{)}^2\\ & =x_{\mu }^{T}Cov(y,y)x_{\mu }/(x_{\mu }^T{x}_{\mu }{)}^2\\ & =x_{\mu }^T{\sigma }^2{I}_n{x}_{\mu }/(x_{\mu }^T{x}_{\mu }{)}^2\\ & ={\sigma }^2/(x_{\mu }^T{x}_{\mu })\end{array}$$
因而，$\hat{a}\sim N(a,{\sigma }^2/(x_{\mu }^T{x}_{\mu }))$ 。又根$SSE/{\sigma }^2$服从自由度n-2的卡方分布，且这两个变量相互独立（这里没给出证明），可以得到统计量T：
$$\begin{array}{rl}T& =\frac{\hat{a}-a}{\sigma /\sqrt{x_{\mu }^T{x}_{\mu }}}/\sqrt{(SSE/{\sigma }^2/(n-2))}\\ & \sim t(n-2)\end{array}$$
$$\begin{gathered} H_0:a=0 \\ {H}_1:a\ne 0 \end{gathered}$$
在空假设$H_0$成立时，T应该接近于0。在$|\hat{a}|>t_{1-\alpha /2}$时拒绝原假设。

一元线性回归的Python编程实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def linear_ols(x, y):
    '''
    实现一元线性回归
    返回参数： a, b, F检验统计量及p值， T检验统计量及p值
    '''
    mu = np.mean(x)
    xmu = x - mu
    a = xmu.dot(y)/xmu.dot(xmu)
    b = np.mean(y) - a*np.mean(x)
    
    SST = np.sum((y-np.mean(y))**2)
    y_pred = a*x + b
    SSE = np.sum((y-y_pred)**2)
    SSR = SST - SSE
    
    N = len(x)
    F = SSR/1 / (SSE/(N-2)) #F统计量
    pF = 1-stats.f.cdf(F, dfn = 1, dfd = N-2)
    
    T =np.sqrt(xmu.dot(xmu)) * a/np.sqrt(SSE/(N-2)) #T统计量
    if(a>0):
        pT = 2*(1-stats.t.cdf(T, df = N-2))
    else:
        pT = 2*stats.t.cdf(T, df=N-2)
    
    return a, b, F, pF, T, pT

生成仿真数据，输出结果。

N = 300
x = np.linspace(-2, 2, N)
y = 3*x + 2*np.random.randn(N)

plt.plot(x, y, 'o')
a, b, F, pF, T, pT = linear_ols(x, y)
y_pred = a*x + b
plt.plot(x, y_pred)
plt.xlabel('x')
plt.ylabel('y')

print('a=%.3f, b=%.3f' %(a, b))
print('F statistics = %.2f, p=%.3f'%( F, pF))
print('T statistics =%.2f, p =%.3f'%( T, pT))

得到的结果如下：

a=3.003, b=-0.100
F statistics = 983.34, p=0.000
T statistics =31.36, p =0.000

由于生成随机数的原因，不同机器运行结果可能略有差异。

与 statsmodels 包的对比

statsmodels提供了线性回归的一般解法，我们将其输出的结果与自己编写的函数进行对比。

from statsmodels import api as sm
X = sm.add_constant(x)#添加常数项
model = sm.OLS(y,X)
result = model.fit()
print(result.summary())

程序输出结果为：

OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.767
Model:                            OLS   Adj. R-squared:                  0.767
Method:                 Least Squares   F-statistic:                     983.3
Date:                Sun, 18 Apr 2021   Prob (F-statistic):           2.18e-96
Time:                        20:39:24   Log-Likelihood:                -620.69
No. Observations:                 300   AIC:                             1245.
Df Residuals:                     298   BIC:                             1253.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0999      0.111     -0.901      0.369      -0.318       0.118
x1             3.0035      0.096     31.358      0.000       2.815       3.192
==============================================================================
Omnibus:                        3.843   Durbin-Watson:                   1.807
Prob(Omnibus):                  0.146   Jarque-Bera (JB):                3.162
Skew:                          -0.147   Prob(JB):                        0.206
Kurtosis:                       2.592   Cond. No.                         1.16
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

从上述结果中可以看出，const也就是我们估计结果中的b值，x1也就是我们估计结果中的a值，两者一致。F-statistics和t-statistics983.3和31.358，也和我们代码的运行结果一致。

关于statsmodels结果的解读，可以参考statsmodels中的summary解读（使用OLS）

关于多元线性回归

推导方法类似，原理类似，个人认为能够通过一元线性回归加深对问题的理解就足够了，多元线性回归在使用时直接使用工具包就行，参数的解读也是相似的。