这节的内容是整个回归分析的核心,后面几节看似是不同的模型,其实都可以转化成 MLE.

这一系列的文章里面很多乱七八糟的矩阵公式,看不懂很正常,多动手就好了。按照文中的思路把每个公式都亲自推导一遍,就能写出来比这个更完善的教程。

基本模型

多元线性回归(MLE)模型是 $$ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta}+\boldsymbol\varepsilon. $$

各矩阵大小为:

  • $\mathbf{Y}: n\times1$
  • $\mathbf{X}: n\times p$
  • $\boldsymbol\beta: p\times n$
  • $\boldsymbol\varepsilon:p\times n$

误差项满足 $\boldsymbol\varepsilon\sim N(0, \sigma^2\mathbf{I})$.

这里面有两个参数:$\boldsymbol{\beta}$ 和 $\sigma^2$。

注意参数的个数,$\boldsymbol\beta$ 的维度是 $p$,意味着有 $p-1$ 个解释变量和 $1$ 个截距项。可能有的教材习惯把解释变量个数记作 $p$,此时 $\boldsymbol\beta$ 就是 $p+1$ 维。

参数估计

$\hat{\boldsymbol\beta}$

回忆高斯-马尔可夫定理(Gauss-Markov Theorem):在线性回归模型中,如果误差满足零均值、同质方差且互不相关,则回归系数的最佳线性无偏估计就是普通最小二乘法估计。

使用最小二乘法,容易算出其最佳线性无偏估计是 $$ \hat{\boldsymbol{\beta}}=(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1}\mathbf{X}^{\mathrm T}\mathbf{Y}. $$

它服从正态分布 $N(\boldsymbol\beta, \sigma^2(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1})$。

$\hat\sigma^2$

$\sigma^2$ 的估计则类似于“样本方差”: $$ \hat\sigma^2=\dfrac{(\mathbf{Y}-\hat{\mathbf{Y}})^{\mathrm T}(\mathbf{Y}-\hat{\mathbf{Y}})}{n-p}. $$

我们知道,独立标准正态分布的平方和是卡方分布,即 $$ \chi^2_k\sim \sum_{i=1}^kN(0,1). $$

$k$ 称为卡方分布的自由度,并且有 ${\mathrm E}(\chi^2_k)=k$。

关于多元正态分布和卡方分布,我们有一个重要结论:设 $X\sim N_p(0,\sigma^2I_p)$,$A$ 是对称幂等矩阵,且 $\mathrm{rank}(A)=r$,则二次型 $X^{\mathrm T}AX/\sigma^2\sim \chi^2_r$。

计算可得,$(\mathbf{Y}-\hat{\mathbf{Y}})^{\mathrm T}(\mathbf{Y}-\hat{\mathbf{Y}})=\boldsymbol\varepsilon^{\mathrm T}(\mathbf I - \mathbf{X}(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1}\mathbf{X}^{\mathrm T})\boldsymbol\varepsilon$。中间那个矩阵是对称幂等矩阵,秩是 $n-p$,而 $\boldsymbol\varepsilon$ 是一个多元正态分布,所以

$$ (n-p)\hat\sigma^2/\sigma^2\sim \chi^2_{n-p}. $$

从而 $\hat\sigma^2$ 是 $\sigma^2$ 的一个无偏估计。

决定系数

预备知识:方差分解

我们来观察这个式子: $$ y = \beta_0+\beta_1x_1+\cdots+\beta_{p-1}x_{p-1}+\varepsilon, $$ 我们得到了一堆 $y_i$,它们互不相同。有两个原因造就了 $y_i$ 之间的差异:

  • $x$ 不同导致 $y$ 不同(由自变量解释的部分)
  • 随机误差 $\varepsilon$

这里的记号非常恶心。有的教材里把 E 当作 explain,R 当作 residual;有的地方又把 E 当作 error,R 当作 regression。这导致 R 和 E 在不同场合表达完全相反的含义。并且有的地方把 SS 放前面,称作“SSE,SSR”,又有的地方放后面,叫“ESS,RSS”。建议谁来规范一下这里的乱象。

骂完了,接着看。$y_i$ 之间的差异性可以用总平方和 $S^2$ 来衡量: $$ S^2=\sum_i (y_i-\bar y)^2. $$ 回归模型做的事情就是让自变量来解释因变量的变化,所以由自变量解释的部分可以用拟合的 $\hat y$ 的均方和 $S_1^2$ 来衡量。由于 $\bar{\hat y}=\bar y$,我们有 $$ S_1^2=\sum_i (\hat y_i-\bar y)^2. $$ 随机误差的影响就是剩下的部分 $S_2^2$: $$ S_2^2=\sum_i(y_i-\hat y_i)^2. $$ 我们有一个非常漂亮的结论: $$ S^2=S_1^2+S_2^2. $$ 证明留给读者作为练习。这一结论成功将 $y_i$ 的差异性分成了两部分:可被回归模型解释的部分和随机误差。

R 平方 (R-squared)

我们定义决定系数(coefficient of determination) $$ R^2=1-\dfrac{S_2^2}{S^2}. $$ 这种定义对于回归模型都适用(包括非线性回归)。我们认为回归之前,$y$ 的变动完全由随机误差造成,也就是随机误差占比为 $1$;而在回归之后,随机误差占比为 $\dfrac{S_2^2}{S^2}$,活活减少了 $R^2$ 这么多,这就是回归模型的功劳!

调整后的 R 平方(Adjusted R-squared)

$R^2$ 有一个缺陷,只要我们无脑增加解释变量的个数,哪怕增加的是随机误差,$R^2$ 也不会下降(通常都是增加),也就是所谓的“力大砖飞”。证明留作习题(只要注意线性回归的优化目标和 $R^2$ 的关系即可)。

调整后的 $R^2$ 增加了一个惩罚项,使得解释变量越多,它的值越小。 $$ R^2_{\rm adj}=1- \dfrac{S_2^2/(n-p)}{S^2/(n-1)}. $$

$S_2^2$ 就是 $(\mathbf{Y}-\hat{\mathbf{Y}})^{\mathrm T}(\mathbf{Y}-\hat{\mathbf{Y}})$,即 $\dfrac{S_2^2}{\sigma^2}\sim \chi^2_{n-p}$. 而 $S^2$ 是正态分布的样本方差,所以 $\dfrac{S^2}{\sigma^2}\sim\chi^2_{n-1}$. 这样一来,所谓“调整”其实就是除以了各自的自由度。

(再次强调一定要搞清楚究竟是 $p$ 还是 $p+1$,这是很多资料记号不一致,容易混淆的地方。)

决定系数和相关系数有啥关系

两组数据 $x_i,y_i$ 的相关系数定义为 $$ \rho(x_i, y_i)=\dfrac{{\rm Cov}(x_i,y_i)}{\sqrt{{\rm Var}(x_i){\rm Var}(y_i)}}. $$ 一般来说,决定系数和相关系数是 Java 和 Javascript 的关系。但对于线性回归模型,我们可以证明: $$ \rho(y,\hat y)=\sqrt{R^2}. $$ 关于相关系数和 R 方的更多讨论可以移步这里。由于不是核心内容,就不过多介绍了。

假设检验

常用的是 $t$ 检验和 $F$ 检验。

$t$ 检验

独立的标准正态分布除以卡方分布的平方根就得到 $t$-分布。

设 $X\sim N(0, 1),\ Y\sim\chi^2_p$,$X, Y$ 独立,则定义 $$ \dfrac{X}{\sqrt{Y/p}} $$ 是自由度为 $p$ 的 $t$-分布。

我们可以从样本均值和样本方差导出 $t$-分布。设 $\bar X, S^2$ 分别是从 $N(\mu, \sigma^2)$ 中抽样得到的 $n$ 个样本的样本均值和样本方差,则由正态分布的性质知:

  • $\bar X$ 和 $S^2$ 独立
  • $\bar X\sim N(\mu, \sigma^2/n)$
  • $(n-1)S^2/\sigma^2\sim \chi^2_{n-1}$

于是 $$ \dfrac{\bar X-\mu}{\sqrt{S^2/n}} $$ 就是自由度为 $n-1$ 的 $t$-分布。

言归正传,我们想要检验某一个回归系数 $\beta_j$ 是不是 $0$。构造统计量 $$ t_j=\dfrac{\hat\beta_j-\beta_j}{\sqrt{\hat\sigma^2(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1}_{jj}}}. $$ 线性回归中有一个结论:$\hat{\boldsymbol{\beta}}$ 和 $\hat\sigma^2$ 独立。因此 $t_j$ 服从自由度为 $n-p$ 的 $t$-分布。

接下来,把 $\beta_j=0$ 和估计出来的 $\hat\beta_j$ 带入,算出 $t$ 值,再查表得到 p-value 即可。

$F$ 检验

将两个独立的卡方分布分别除以自由度,然后再相除,可以得到 $F$ 分布: $$ F(d_1,d_2)=\dfrac{\chi^2_1/d_1}{\chi^2_2/d_2}. $$

对于模型 $y=\beta_0+\beta_1x_1+\cdots+\beta_{p-1}x_{p-1}+\varepsilon$,$F$ 检验零假设为 $$ H_0:\beta=\beta_1=\cdots=\beta_{p-1}=0. $$ 即各个回归系数都是 $0$(不包括 $\beta_0$)。如果 $H_0$ 成立就房倒屋塌完犊子了,满地飘零。我们构造的统计量是 $$ F=\dfrac{S_1^2/(p-1)}{S_2^2/(n-p)}. $$

$S_2^2$ 的分布前面讨论多次了。而 $S_1^2=(\hat{\mathbf{Y}}-\bar{\mathbf{Y}})^{\mathrm T}(\hat{\mathbf{Y}}-\bar{\mathbf{Y}})$. 记 ${\mathbf 1}_n$ 为 $n\times n$ 的全 $1$ 矩阵,在 $H_0$ 成立时有 $\hat{\mathbf{Y}}-\bar{\mathbf{Y}}=\left(\mathbf{X}(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1}\mathbf{X}^{\mathrm T}-\frac{1}{n}\mathbf{1}_n\right)\boldsymbol\varepsilon $ .

可以验证,$\mathbf{X}(\mathbf{X}^{\mathrm T}\mathbf{X})^{-1}\mathbf{X}^{\mathrm T}-\frac{1}{n}\mathbf{1}_{n}$ 是对称幂等矩阵,秩为 $p-1$(该命题的证明可作为经典的线性代数题目,请自行练习)。

于是 $S_1^2/\sigma^2\sim\chi^2_{p}$. 从而 $F\sim F(p-1, n-p)$.

这个统计量的含义便是看 $S_1^2$ 和 $S_2^2$ 谁占上风。$F$ 越大,表明能被回归解释的部分越大,所以线性模型是成功的,也就越有可能拒绝 $H_0$。

区别与联系

简单来说,$F$ 检验是判断线性模型是否正确,如果不显著,可能要考虑换其他模型;而 $t$ 检验用来判断模型中是否该有这个变量。

$t$ 检验判断单个系数是否为 0,$F$ 检验判断整体是否为 0. 看起来是重复的——我对每个系数都做一次 $t$ 检验不就可以代替 $F$ 检验了吗?其实不然。有以下两个原因:

  • 假设检验的个数越多,第一类错误概率越大
  • 多重共线性(Multicollinearity)会导致所有 $t$ 检验显著时,$F$ 检验不显著

残差分析

回归模型中的一个基本假设是 $\varepsilon\overset{i.i.d.}\sim N(0,\sigma^2)$. 如果这个假设不成立,后面所做的东西就有可能完犊子。残差分析就是在做完回归之后检查一下关于 $\varepsilon$ 的假定是否成立。

残差图

残差(residual)就是 $r_i=y_i-\hat y_i$. 它可以看成对 $\varepsilon_i$ 的估计。以残差为纵坐标,其他量为横坐标作散点图,就叫做残差图(residual plots)。残差图是进行残差分析的有力武器。

常用的是以 $\hat y$ 作为横坐标和以 $x$ 作为横坐标的两种图。由于 $\hat y$ 和 $x$ 之间有线性关系,所以这两者的分析是几乎一样的。如下图(图片来自 https://blog.csdn.net/mengjizhiyou/article/details/82216278):

(a)表示的是理想的残差图,均匀分布在平行带中,没有任何可预测的信息。(b)中残差的波动范围与 $x$ 有关系,意味着每个 $\varepsilon_i$ 的方差可能不一致。(c)虽然方差一致,但是均值在变化,说明模型没有很好地解释 $y$ 的变化,需要考虑引入高次项或更换模型。

Q-Q 图

Q-Q 图是 quantile-quantile 的意思,即横纵坐标都是分位数。在残差分析中,横坐标用正态分布的分位数,纵坐标用残差的分位数画散点图。具体来说,残差有 $n$ 个,将其递增排序为 $r_{(1)},\cdots,r_{(n)}$,然后算出正态分布在 $(i-0.5)/n$ 处的分位数 $z_1,\cdots,z_n$, 画出 $(z_i,r_{(i)})$ 的散点图就是 Q-Q 图。

显然,理想的图应该分布在 $y=x$ 上,偏了或者弯了都有毛病。

结语

本文只是快速浏览了一遍线性回归的相关知识,很多地方写的不详细。复习时可以以此为线索,遇到不熟悉的地方再去查找其他资料。里面提到的没有证明的结论也最好自己亲自推导一下,很有意思。