LRA系列的第一篇，介绍单变量线性回归模型

手敲\(\LaTeX\)难免出现纰漏，有任何疑似错误或者不清楚的地方请直接在下方评论区留言，谢谢各位读者。

模型假设

\[Y_i=\beta_0+\beta_1 X_i+\varepsilon_i\]

其中\(\varepsilon_i\sim\mathcal{N}(0,\,\sigma^2)\)。当然有时也可以弱化到\(\mathbb{E}\{\varepsilon_i\}=0,\,\sigma^2\{\varepsilon_i\}=\sigma^2\)。但出于简单考虑，我们后文均假设\(\varepsilon_i\sim\mathcal{N}(0,\,\sigma^2)\)。

一般地，我们称\(X\)为自变量(independent variable or predictor)，\(Y\)为因变量(dependent variable or response)，\(\varepsilon\)为误差(random error or noise)

此时我们的误差函数为

\[Q=\sum_{i=1}^{n}(Y_i-\beta_0-\beta_1 X_i)^2\]

模型求解

最小化该函数立刻模型的最小二乘估计(least squares estimator, LSE)为

\[b_1=\widehat{\beta_1}=\frac{\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\] \[b_0=\widehat{\beta_0}=\bar{Y}-b_1\bar{X}\]

可以证明\(b_0,\,b_1\)是\(\beta_0,\,\beta_1\)的无偏估计。[proof1]

Gauss-Markov theorem 在更弱的假设下，甚至证明了\((b_0,\,b_1)\)是\(\beta_0,\,\beta_1\)一致最小方差无偏估计(uniformly minimum-variance unbiased estimator, UMVUE)

此时，我们通常记拟合值(fitted value)为\(\widehat{Y}_i=b_0+b_1 X_i\)，残差(residual) \(e_i=Y_i-\widehat{Y}_i{}__{}\)，残差平方和(error sum of squares)\(\mathrm{SSE}=\sum_{i=1}^{n}e_i^2\)

可以证明有如下性质[proof2]

\[\sum_{i=1}^{n} e_i=\sum_{i=1}^{n}X_i e_i = 0\]

参数估计

方差的估计

首先考虑更简单的模型\(Y_i=\mu+\varepsilon_i\), 也就是给出一组数据\(Y_i\)求其方差。可以证明最小二乘估计

\[\widehat{\mu}=\bar{Y}\]

且此时\(\sigma^2\)的无偏估计为

\[s^2=\frac{\sum_{i=1}^{n}(Y_i-\widehat{Y}_i)^2}{n-1}\]

其证明见[proof3]。一个直观的理解是，我们使用样本方差估计分布的方差总是低估的。比如单个样本的方差为0，两个样本的方差计算时使用的是到二者均值的差的平方和而非分布的均值，这包含着\(Y\)的信息并使得方差被低估了。从自由度来说，计算\(\mathrm{SSE}\)的时候使用了\(\widehat{Y}\)这个信息，使得其自由度仅有\(n-1\)

回到我们的线性回归模型，可以证明\(\sigma^2\)的无偏估计为[proof4]

\[\mathrm{MSE}=\frac{\mathrm{SSE}}{n-2}=\frac{\sum_{i=1}^{n}(Y_i-\widehat{Y}_i)^2}{n-2}\]

从自由度来说，计算\(\mathrm{SSE}\)的时候使用了\(b_0,\,b_1\)这两个信息，使得其自由度仅有\(n-2\)

\(\beta_1\)的估计

可以证明[proof5], \(b_1\sim\mathcal{N}(\beta_1,\sigma^2\{b_1\})\), 其中

\[\sigma^2\{b_1\}=\sigma^2/\sum_{i=1}^{n}(X_i-\bar{X})^2\]

实际中由于\(\sigma^2\)未知，采用其无偏估计，即

\[s^2\{b_1\}=\frac{\mathrm{SSE}}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\]

此时

\[\frac{b_1-\beta_1}{s\{b_1\}}\sim t_{n-2}\]

由此可得\(\beta_1\)的置信区间

\(\beta_0\)的估计

可以证明[proof6], \(b_0\sim\mathcal{N}(\beta_0,\sigma^2\{b_0\})\), 其中

\[\sigma^2\{b_0\}=\sigma^2\left(\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\]

实际中由于\(\sigma^2\)未知，采用其无偏估计，即

\[s^2\{b_0\}=\mathrm{MSE}\left(\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\]

此时

\[\frac{b_0-\beta_0}{s\{b_0\}}\sim t_{n-2}\]

由此可得\(\beta_0\)的置信区间

值得注意的是，当\(|\bar{X}|\)较大时，\(s^2\{b_0\}\)较大，意味着\(\beta_0\)较难估计。

\(\widehat{Y}_h\)的估计

可以证明[proof7], \(\widehat{Y}_h\sim\mathcal{N}(\beta_0+\beta_i X_h,\sigma^2\{hat{Y}_h\})\), 其中

\[\widehat{Y}_h=b_0+b_1 X_h,\quad \sigma^2\{hat{Y}_h\}=\sigma^2\left(\frac{1}{n}+\frac{(X_h-\bar{X})^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\]

实际中由于\(\sigma^2\)未知，采用其无偏估计，即

\[s^2\{\widehat{Y}_h\}=\mathrm{MSE}\left(\frac{1}{n}+\frac{(X_h-\bar{X})^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\]

此时

\[\frac{\widehat{Y}_h-(\beta_0+\beta_1 X_h)}{s\{\widehat{Y}_h\}}\sim t_{n-2}\]

由此可得\(\widehat{Y}_h-\)的置信区间。

值得注意的是，当\(|X_h-\bar{X}|\)较大时，\(s^2\{\widehat{Y}_h\}\)较大，意味着\(\beta_0+\beta_1 X_h\)较难估计。

事实上，\(\beta_0\)的估计就是其当\(X_h=0\)的特殊情况。

\(Y_h\)的估计

和上一小节的区别是，上一小节预测的是在\(X_h\)处平均响应的情况，而这里是\(X_h\)处单次响应的情况，差别就是\(Y_h=\widehat{Y_h}+\varepsilon\)的差别。

可以证明[proof8]

\[s^2\{Y_h\}=s^2\{\widehat{Y}_h\}+\mathrm{MSE}\]

由此可得\(Y_h\)的置信区间。

通常，我们将上一小节求得的置信区间称为置信区间(Confidence Interval, CI)，而将这一小节的称为预测区间(Prediction Interval, PI)。

同样的，在上面的计算中，我们发现\(\sum_{i=1}^{n}(X_i-\bar{X})^2\)越大其预测值的方差越小，这提示我们当数据的分布范围较大时，其效果相对更好。

方差分析

事实上，方差分析会在后面单独用一篇笔记来介绍，在这里只是进行结论性的介绍。

由于\(e_i=Y_i-\widehat{Y}=(Y_i-\widehat{Y}_i)+(\widehat{Y_i}-\bar{Y})\)。我们一般称前者为随机误差((至少现阶段)无法消除)，而后者为系统误差(通过我们的模型可以消除)。通过简单的平方展开不难得到

\[\sum_{i=1}^{n}(Y_i-\widehat{Y})^2=\sum_{i=1}^{n}(Y_i-\widehat{Y}_i)^2+\sum_{i=1}^{n}(\widehat{Y_i}-\bar{Y})^2\] \[\mathrm{SSTO}:(n-1)=\mathrm{SSE}:(n-2) + \mathrm{SSR}:1\]

\(F\)检验

我们先使用\(F\)检验去检验\(H_0:\beta_1=0\)。可以证明[proof9] \[\mathbb{E}\{\mathrm{MSR}\}=\sigma^2+\beta_1^2\sum_{i=1}^{n}(X_i-\bar{X})^2\] \[\mathbb{E}\{\mathrm{MSE}\}=\sigma^2\]

故可取

\[F^\ast=\frac{\mathrm{MSR}}{\mathrm{MSE}}\sim F_{1,n-2}\]

当\(F^\ast>F_{1-\alpha,1,n-2}\)时拒绝原假设。

\(t\)检验

注意到下列等式关系就能发现\(F\)检验和\(t\)检验是等价的 \[\frac{\mathrm{MSR}}{\mathrm{MSE}}=F^\ast=(t^\ast)^2=\left(\frac{b_1}{s\{b_1\}}\right)^2\] \[F_{1-\alpha,1,n-2}=t_{1-\alpha/2,n-2}^2\]

事实上，该关系只在一元线性回归中成立。在多元线性回归中，\(t\)检验针对的是某个系数是否为0，而\(F\)检验针对的是所有系数是否为零。在接下来的章节中读者应当会有更深的体会。

残差

决定系数与相关系数

决定系数(coefficient of determination)为\(R^2=\mathrm{SSR}/\mathrm{SSTO}\)，衡量模型所能解释的方差比例。

相关系数(coefficient of correlation)为\(r=corr(\boldsymbol{X},\,\boldsymbol{Y})\)，衡量\(\boldsymbol{X},\,\boldsymbol{Y}\)的线性相关性。

二者有关系\(R^2=r^2\)

残差的性质

从最小二乘估计中我们得到

\[\sum_{i=1}^{n} e_i=\sum_{i=1}^{n}X_i e_i = 0\] \[\mathrm{MSE}=\sum_{i=1}^{n}e_i^2/(n-2)\]

从模型假设中我们得到

\[\mathbb{E}\{e_i\}=0\] \[\sigma^2\{e_i\}=\sigma^2\left(1-\left(\frac{1}{n}+\frac{(X_h-\bar{X})^2}{\sum_{i=1}^{n}(X_i-\bar{X})^2}\right)\right)=\sigma^2(1-h_i)\]

其中\(\sigma^2\{\widehat{Y}_i\}=\sigma^2h_i\)

标准化残差(standardized residuals)\(e_i^\ast=e_i/\sqrt{\mathrm{MSE}}\)，可被近似地视为服从标准正态分布。

而学生化残差(studentized residuals)\(r_i=e_i/\sqrt{\mathrm{MSE}(1-h_i)}\)，则严格服从\(t_{n-2}\)，是更精细的一个构造

模型诊断

可能出现的问题

  1. \(Y\)与\(X\)非线性关系
  2. \(\sigma^2\{\varepsilon_i\}\)并不能视为常数
  3. \(\varepsilon_i\)彼此并不独立
  4. 模型能够拟合但有离群值
  5. \(\varepsilon_i\)不服从正态分布
  6. 遗漏了某些重要的自变量

常用诊断方式

  1. 使用\(e_i\sim X_i\)诊断问题1,4
  2. 使用\(e_i\sim \widehat{Y}_i\)诊断问题2,4
  3. 使用\(e_i\sim e_{i+1}\)诊断问题3
  4. 使用qqplot等诊断问题5
  5. 使用\(e_i\sim X_{new}\)诊断问题6

一个简单的变种——过原点模型

计算方式和之前几乎无异，这里仅给出结果供参考。

\[Y_i=\beta_1 X_i+\varepsilon_i\]

\[b_1=\frac{\sum_{i=1}^{n}X_iY_i}{\sum_{i=1}^{n}X_i^2}\] \[\widehat{\sigma^2}=\mathrm{MSE}=\frac{\mathrm{SSE}}{n-1}\] \[\mathrm{SSE}=\sum_{i=1}^{n}Y_i^2-\frac{(\sum_{i=1}^{n}X_iY_i)^2}{\sum_{i=1}^{n}X_i^2}\] \[s^2\{b_1\}=\frac{\mathrm{SSE}}{\sum_{i=1}^{n}X_i^2}\] \[s^2\{\widehat{Y}_h\}=\frac{\mathrm{SSE}X_h^2}{\sum_{i=1}^{n}X_i^2}\] \[\frac{b_1-\beta_1}{s\{b_1\}}\sim t_{n-1}\] \[\frac{\widehat{Y}_h-\beta_1X_h}{s\{\widehat{Y}_h\}}\sim t_{n-1}\]

矩阵化描述——为多元线性回归铺路

如果读者对矩阵化描述熟练的话，这一节将是无比平凡的。否则，我建议读者将这里面的公式展开写写，并和上面的结果进行对比。

尽管本节可能十分平凡，但其内容的重要性丝毫不比前面少，因为这可以毫无障碍地拓展到多元的情形。

模型假设

对于 \(Y_i=\beta_0+\beta_1 X_i+\varepsilon_i\)，其矩阵形式为

\[ \begin{pmatrix} Y_1 \\ \vdots \\ Y_n \\ \end{pmatrix} = \begin{pmatrix} 1&X_1 \\ \vdots&\vdots \\ 1&X_n \\ \end{pmatrix} \begin{pmatrix} \beta_0\\ \beta_1 \\ \end{pmatrix} + \begin{pmatrix} \varepsilon_1 \\ \vdots \\ \varepsilon_n \\ \end{pmatrix} \]

我们通常用黑体表示向量、矩阵（甚至有教材对此用斜体和直体进一步区分，但本文并不特殊对待），即

\[\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}\]

模型假设的矩阵形式则为

\[\mathbb{B}\{\boldsymbol{Y}\}=\boldsymbol{X}\boldsymbol{\beta}\] \[\mathbb{E}\{\epsilon\}=\boldsymbol{0}\] \[\sigma^2\{\boldsymbol{Y}\}=\sigma^2\{\boldsymbol{\varepsilon}\}=\sigma^2\boldsymbol{I}\]

模型结果

\[Q=(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})^\intercal(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})\]

系数

\[\boldsymbol{b}=(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1}\boldsymbol{X}^\intercal\boldsymbol{Y}\]

\[ \begin{pmatrix} b_0 \\ b_1 \\ \end{pmatrix}= \begin{pmatrix} n&\sum_{i=0}^{n}X_i \\ \sum_{i=0}^{n}X_i&\sum_{i=0}^{n}X_i^2 \\ \end{pmatrix}^{-1} \begin{pmatrix} \sum_{i=0}^{n}Y_i\\ \sum_{i=0}^{n}X_iY_i \\ \end{pmatrix} \]

\[\mathbb{E}\{\boldsymbol{b}\}=(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1}\boldsymbol{X}^\intercal\boldsymbol{X}\boldsymbol{\beta}=\boldsymbol{\beta}\] \[\sigma^2{\boldsymbol{b}}=\sigma^2(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1}\]

拟合值

\[\boldsymbol{\widehat{\boldsymbol{Y}}}=\boldsymbol{X}(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1}\boldsymbol{X}^\intercal\boldsymbol{Y}\overset{\Delta}{=}\boldsymbol{H}\boldsymbol{Y}\] \[\mathbb{E}\{\widehat{\boldsymbol{Y}}\}=\boldsymbol{H}\boldsymbol{X}\boldsymbol{\beta}=\boldsymbol{X}\boldsymbol{\beta}\] \[\sigma^2{\widehat{\boldsymbol{Y}}}=\boldsymbol{H}(\sigma^2\boldsymbol{I})\boldsymbol{H}=\sigma^2\boldsymbol{H}^2=\sigma^2\boldsymbol{H}\]

残差

\[\boldsymbol{e}=\boldsymbol{Y}-\boldsymbol{\widehat{Y}}=(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\] \[\mathbb{E}\{\boldsymbol{e}\}=(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{X}\boldsymbol{\beta}=\boldsymbol{0}\] \[\sigma^2{\boldsymbol{e}}=\sigma^2(\boldsymbol{I}-\boldsymbol{H})\] \[\mathrm{SSE}=\boldsymbol{e}^\intercal\boldsymbol{e}=\boldsymbol{Y}^\intercal(\boldsymbol{I}-\boldsymbol{H})\boldsymbol{Y}\]

预测值

\[\boldsymbol{X}_h=(1,X_h)^\intercal\] \[\boldsymbol{\widehat{Y}}_h=\boldsymbol{X}_h^\intercal\boldsymbol{b}\] \[\mathbb{E}\{\boldsymbol{\widehat{Y}}_h\}=\boldsymbol{X}_h^\intercal\boldsymbol{\beta}\] \[\sigma^2\{\widehat{\boldsymbol{Y}}_h\}=\sigma^2\boldsymbol{X}_h^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})\boldsymbol{X}_h\]