这节的内容是整个回归分析的核心，后面几节看似是不同的模型，其实都可以转化成 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$ 上，偏了或者弯了都有毛病。
结语
本文只是快速浏览了一遍线性回归的相关知识，很多地方写的不详细。复习时可以以此为线索，遇到不熟悉的地方再去查找其他资料。里面提到的没有证明的结论也最好自己亲自推导一下，很有意思。