【数理统计基础】 05 - 回归分析
参数估计和假设检验是数理统计的两个基础问题，它们不光运用于常见的分布，还会出现在各种问题的讨论中。本篇开始研究另一大类问题，就是讨论多个随机变量之间的关系。现实生活中的数据杂乱无章，够挖掘出各种变量之间的关系非常有用，它可以预估变量的走势，能帮助分析状态的根源。关系分析的着手点可以有很多，我们从最简单直观的开始，逐步展开讨论。
1. 一元线性回归
1.1 回归分析
如果把每个量都当做随机变量，问题的讨论会比较困难，或者得到的结论会比较受限。一个明智做法就是只把待考察的量\(Y\)看做随机变量，而把其它量\(X_i\)看成是自主选定的。即使都看成变量，也是把\(Y\)看成因变量，而把\(X_i\)看成自变量。该模型同样是研究某个随机变量的情况，不同之处在于更加关注变量与各因素的函数关系，希望能找到影响随机变量的主要因素并给出表达式。
如式（1）所示，选定要关注的因素\(X_i\)，并假定它们以函数\(f(X_1,\cdots,X_p)\)形式影响变量\(Y\)，其它的因素统一放到随机变量\(e\)中。其中函数\(f\)称为\(Y\)对\(X_i\)的回归函数或回归方程，\(e\)则是随机误差。由于已经提取出主要因素，这里假定\(e\)的均值为\(0\)，并且它是与\(f\)独立的。在应用场景，一般给定回归函数一个含参表达式（比如后面的线性回归），这样的问题称为参数回归问题，否则叫非参数回归问题。
\[Y=f(X_1,X_2,\cdots,X_p)+e,\;\;E(e)=0,\;0<D(e)<\infty\tag{1}\]
在微积分中我们知道，多元函数\(f(x_1,\cdots,x_p)\)一般可以在任何点进行泰勒展开，其中最简单的就是线性展开。线性关系由于其形式简单，以及在局部能很好地逼近函数，在数学的各分支都被重点讨论。在回归分析中，这样的模型便称为线性回归，这里先从最简单的一元线性回归讨论起。
一元线性回归的模型是式（2）左，在提取出线性关系\(b_0+b_1X\)后，\(Y\)的剩余因素或随机性就都落在随机变量\(e\)上。所以从另一个角度看，回归分析是要找出随机变量的“确定”部分和随机部分，这种分解更能帮助分析随机现象。自然地，分析是基于\(n\)个样本点\((X_i,Y_i)\)，其中\(X_i\)可能也是随机产生的，但在这个模型里一律看做定量。还要注意，这时每个\(Y_i\)是\(e_i\)的一个位移，它不再与\(Y\)同分布。
\[Y=b_0+b_1X+e;\;\;Y_i=b_0+b_1X_i+e_i=a_0+a_1(X_i-\bar{X})+e_i\tag{2}\]
1.2 系数点估计
每一次试验相互独立，因此得到\(n\)个独立变量（式（2）右），其中\(b_0,b_1\)是待定系数。为了方便计算和讨论，一般还会把式（2）右的线性部分“中心化”。问题等价于讨论\(a_0,a_1\)的值，但要注意这里\(\bar{X}\)是依赖于具体样本的。在得知样本点\((X_i,Y_i)\)的情况下，如何确定系数比较合理？在式（2）中，我们把\(Y_i\)看做是有误差\(e_i\)的变量，因此让误差的平方和达到最小是一个比较好的模型。
式（3）取最小值时\(a_0,a_1\)便是合理的参数估计，利用偏导为零容易算得式（4）中对\(a_0,a_1\)的估计，这个结论非常得益于刚才的中心化。求解的方法其实就是最小二乘法，这在后面再展开讨论。\(a_0\)表示\(Y\)的中心，估计值\(\alpha_0\)十分合理。\(a_1\)应当是\(Y\)关于\(X\)的斜率，单点的斜率是\(\dfrac{Y_i-\bar{Y}}{X_i-\bar{X}}\)，将分子分母同时乘以\(X_i-\bar{X}\)并相加，化简后便得到\(\alpha_1\)，故它也是斜率的合理估计。
\[\sum\limits_{i=1}^n(Y_i-b_0-b_1X_i)^2\tag{3}\]
\[\alpha_0=\bar{Y};\;\;\alpha_1=\sum_{i=1}^n\dfrac{X_i-\bar{X}}{S^2}Y_i,\;S^2=\sum_{i=1}^n(X_i-\bar{X})^2\tag{4}\]
另外还要注意，式（4）中\(X_i\)是定值，而\(Y_i\)独立随机变量，\(\alpha_0,\alpha_1\)都是\(Y_i\)的线性函数，这对于下面的讨论很重要。估计合理的另一个基本要求应当是误差估计、即统计量（随机变量）\(\alpha_0,\alpha_1\)的期望值应当就是\(a_0,a_1\)，利用式（5）左很容易验证结论成立（式（5）右）。以下令\(e\)的方差为\(\sigma^2\)，利用\(Y_i\)的无关性也容易有式（6）。其中\(D(\alpha_1)\)分母中的\(S^2\)有直观的含义，当\(X_i\)比较分散时，得到的斜率估计越准确。另外还可以证明，\(\alpha_0,\alpha_1\)是\(a_0,a_1\)的MVU估计。
\[E(Y_i)=a_0+a_1(X_i-\bar{X})\;\Rightarrow\;E(\alpha_0)=a_0,\;E(\alpha_1)=a_1\tag{5}\]
\[D(Y_i)=\sigma^2\;\Rightarrow\;D(\alpha_0)=\dfrac{\sigma^2}{n},\;D(\alpha_1)=\dfrac{\sigma^2}{S^2}\tag{6}\]
还有一点，把\(\alpha_0,\alpha_1\)看成\(Y_i\)的线性函数，观察两者的“系数向量”，发现它们内积为\(0\)。从向量的角度它们就是直交的，经验证\(\alpha_0,\alpha_1\)也的确是（线性）不相关的，这个结论非常重要，也显示了前面中心化的意义。另外，当\(e\)是正态分布时，\(\alpha_0,\alpha_1\)也都是正态分布，故可知它们独立。
1.3 误差估计
对于模型（2）来说，目前还有\(e\)的方差\(\sigma^2\)没有讨论，在有了系数估计（4）后，现在来估计误差的方差。随着\(X\)的变化，\(Y\)的中心也跟着变化，其误差的方差自然也要以具体的中心为准。在样本点\((X_i,Y_i)\)处，误差\(\delta_i\)（式（7））也称为残差，它们的平均平方和理应作为方差的估计。但由于\(\alpha_0,\alpha_1\)的估计中消耗了两个自由度，故可验证式（7）才是\(\sigma^2\)的无偏估计。
\[\hat{\sigma}^2=\dfrac{1}{n-2}\sum_{i=1}^n\delta_i^2,\;\;\delta_i=Y_i-\alpha_0-\alpha_1(X_i-\bar{X})\tag{7}\]
具体计算步骤参考教材（或自行证明），结果是得到式（8），这样就不难得到是（7）了。当然在实际计算时，可以直接展开得到式（9），然后利用现成的\(X_i,Y_i,\alpha_j\)来加速计算。而且从式（9）中还能得到更有用的结论，注意其中的后两项\(n\bar{Y}^2=Z_1^2\)和\(S^2\alpha_1^2=Z_2^2\)，\(Z_1,Z_2\)都是\(Y_i\)的线性函数，且系数向量是两个相互正交的标准化向量。
\[\sum_{i=1}^n\delta_i^2=\sum_{i=1}^n(e_i-\bar{e})^2+\dfrac{1}{S^2}\left(\sum_{i=1}^n(X_i-\bar{X})e_i\right)^2\tag{8}\]
\[\sum_{i=1}^n\delta_i^2=\sum_{i=1}^nY_i^2-n\bar{Y}^2-S^2\alpha_1^2\tag{9}\]
当\(e\)是正态分布时，\(Y_i\)也是正态分布，利用正交变换的性质，易知式（9）等于\(Z_3^2+\cdots+Z_n^2\)，其中\(Z_j\sim N(0,\sigma^2)\)，这便容易有式（10）的结论。关于残差，还有两点需要注意，式（8）如果很大或者残差体现出某些规律性，则说明线性模型不太合适，或还有重要因素没有被提取出来。
\[e\sim N(0,\sigma^2)\;\Rightarrow\;\sum_{i=1}^n\dfrac{\delta_i^2}{\sigma^2}\sim\chi_{n-2}^2\tag{10}\]
1.4 区间估计
有了点估计，便可以做区间估计，为了能使用枢轴函数，这里还是假定\(e\)为正态分布。首先由公式（5）（6）可知\(\alpha_0,\alpha_1\)满足式（11）的分布，当\(\sigma\)已知时，枢轴函数很容易得到。当\(\sigma\)未知时，由刚才的讨论知\(\alpha_0,\alpha_1\)与\(\hat{\sigma}^2\)是相互独立的，这样便能用\(\hat{\sigma}\)替代\(\sigma\)，得到式（12）的枢轴变量。
\[\alpha_0\sim N(a_0,\dfrac{\sigma^2}{n});\;\;\alpha_1\sim N(a_1,\dfrac{\sigma^2}{S^2})\tag{11}\]
\[\dfrac{\sqrt{n}(\alpha_0-a_0)}{\hat{\sigma}}\sim t_{n-2};\;\;\dfrac{S(\alpha_1-a_1)}{\hat{\sigma}}\sim t_{n-2}\tag{12}\]
线性回归的目的自然是为了进行预测，但在仅知道样本点且把\(X_i\)看成定量的情况下，其实是无法估计最初式（2）左中的\(b_0,b_1\)的。因此要注意，在用\(y=a_0+a_1(x-\bar{X})\)预估\(Y\)时，我们不光丢失了误差\(e\)，还丢失了\(X\)非连续得来的误差。前者通过合理建模来降低误差，后者则只能通过增加\(X_i\)的数量和密度来降低误差。
这一点容易通过估计值\(y\)的方差看出（式（13））。首先在样本数不变的情况下，\(x\)离样本中心\(\bar{X}\)越近方差越小，这个结论符合直觉，样本离预测点越近精度越高。另一方面，样本数越大方差也越小，这个很好理解。结合这两方面，当\(n\)足够大且\(x\)离样本中心足够近，估计的方差就可以任意小。
\[D(y)=\left(\dfrac{1}{n}+\dfrac{(x-\bar{X})^2}{S^2}\right)\sigma^2\tag{13}\]
2. 多元线性回归
2.1 系数估计
现实中的因变量可能受多个因素的影响，这些因素可能有主次之分，也可能是联合作用。无论如何，为了对因变量进行更加深入细致的分析，必须加入更多的自变量进行分析。另外同样的道理，多元函数在局部都可以用线性函数很好地近似，因此我们也可以建立式（14）中的模型和中心化样本表达式。为表达方便，本段下面就直接把\(X_{ki}-\bar{X}_k\)记作\(X_{ki}\)。
\[Y=b_0+\sum_{k=1}^p b_kX_k+e;\;\;Y_i=a_0+\sum_{k=1}^p a_k(X_{ki}-\bar{X}_k)+e_i\tag{14}\]
多元模型的讨论内容和方法与一元的差别不大，但直接的讨论会很繁琐，必须借助于线性代数的工具，请注意前后对比。为讨论方便，首先规定式（15）的简写，并记\(\gamma\)的点估计为\(\alpha\)。然后定义列向量的期望\(E(\alpha)=[E(\alpha_i)]\)，以及协方差\(\text{Cov}(\alpha,\beta)=[\text{Cov}(\alpha_i,\beta_j)]\)，且不难验证有式（16）成立。其实利用算子理论证明会很简单，但光凭形式化的假设，也不难完成证明，请独立尝试。
\[\beta=\begin{bmatrix}Y_1\\Y_2\\\vdots\\Y_n\end{bmatrix},\;\gamma=\begin{bmatrix}a_0\\a_1\\\vdots\\a_p\end{bmatrix},\;A=\begin{bmatrix}1&\cdots&1\\X_{11}&\cdots&X_{1n}\\\vdots&\ddots&\vdots\\X_{p1}&\cdots&X_{pn}\end{bmatrix}\tag{15}\]
\[E(A\alpha)=AE(\alpha);\;\;\text{Cov}(A\alpha,B\beta)=A\text{Cov}(\alpha,\beta)B^T\tag{16}\]
有了矩阵的定义，就可以直接利用线性代数中最小二乘的结论，得到（17）左的正则方程，以及式（17）的\(\gamma\)最小二乘解。式（18）推导出\(\alpha_i\)是\(a_i\)的无偏估计，且都是\(Y_i\)的线性函数，继而还可以得到式（19）的协方差公式。注意到\(A\)中除第一列外，每行的和都是\(0\)，故\(L\)及\(L^{-1}\)都具有形式\(\begin{bmatrix}1&0\\0&B_{p\times p}\end{bmatrix}\)。这说明\(\alpha_0=\bar{Y}\)与其它\(\alpha_i\)互不相关，这与一元的情况是一致的。
\[L\alpha=A\beta\;\Rightarrow\;\alpha=L^{-1}A\beta,\;\;(L=AA^T)\tag{17}\]
\[E(\alpha)=L^{-1}AE(\beta)=L^{-1}AA^T\gamma=\gamma\tag{18}\]
\[\text{Cov}(\alpha,\alpha)=L^{-1}A\text{Cov}(\beta,\beta)A^TL^{-1}=\sigma^2L^{-1}\tag{19}\]
2.2 误差估计
现在来分析误差\(e\)，首先记残差向量\(\delta=\beta-A^T\alpha\)，容易证明\(E(\delta)=0\)，而且根据式（20）的推导可知\(\alpha_i\)与\(\delta_j\)互不相关。另外可以算得式（21）的协方差，其中\(B\)是一个秩为\(p+\)的非负定方阵。根据\(B^2=B\)可以证明，\(B\)的\(p+1\)个特征值都是\(1\)，从而它的迹\(tr(B)=p+1\)（主对角线之和，请参考线性代数）。
\[\text{Cov}(\hat{\alpha},\delta)=L^{-1}A\text{Cov}(\beta,\beta)(I_n-A^TL^{-1}A)=0\tag{20}\]
\[\text{Cov}(\delta,\delta)=(I_n-B)\text{Cov}(\beta,\beta)(I_n-B)=\sigma^2(I_n-B),\;\;(B=A^TL^{-1}A)\tag{21}\]
为了估计\(\sigma^2\)，自然想到讨论残差平方和\(\sum\limits_{i=1}^n\delta_i^2\)。式（22）计算了它的期望值，这样就可以用式（23）来无偏估计\(\sigma^2\)。
\[E(\sum_{i=1}^n\delta_i^2)=\sum_{i=1}^nD(\delta_i)=tr(\text{Cov}(\delta,\delta))=\sigma^2(n-p-1)\tag{22}\]
\[\hat{\sigma}^2=\dfrac{1}{n-p-1}\sum_{i=1}^n\delta_i^2\tag{23}\]
残差平方和\(\delta^T\delta\)是一个半正定二次型，展开整理后有式（24）成立，它满足柯赫伦定理的条件。故假定\(e\)为正态分布的情况下，有式（25）左成立。另外由于\(\alpha_i\)与\(\delta_j\)互不相关，则\(\alpha_i\)与\(\hat{\sigma}\)也不相关，正态分布下它们还是相互独立的，这就得到式（25）右的枢轴变量。
\[\beta^T\beta=\beta^TB\beta+\delta^T\delta\tag{24}\]
\[\sum_{i=1}^n\dfrac{\delta_i^2}{\sigma^2}\sim\chi_{n-p-1};\;\;\dfrac{\alpha_i-a_i}{\sqrt{L_{ii}^{-1}}\hat{\sigma}}\sim t_{n-p+1}\tag{25}\]
2.3 假设检验
线性回归的假设往往是针对线性系数\(a_k\)的，如果仅是对单个系数的检验，直接利用式（25）的枢轴变量即可。实际应用中最常用的假设是\(a_k=0\)，它说明因素\(X_k\)对\(Y\)其实是不相关的，这对检验变量相关性很有用（但更偏重\(X_k\)对\(Y\)的影响）。观察式（17），你会发现\(\alpha_k\)并不只与\(X_k,Y\)有关，它与上面的一次模型得到的结论不一样。可以这样解释：更多因素的加入使得模型更加精确。
但是不是因素越多越好？如果加入的是真正影响\(Y\)的元素，对模型自然是有益的，否则多加入的元素只能增加随机性，从而对结论精度造成影响。样本不足的情况下，以上模型容易把无效元素估计成“假”的关系，从而影响真实因素的作用。但逐个地检验无效元素，有时效果并不好，因为元素之间的复杂关系和随机性会使得检验出现较大偏差。
检验较多无关参数时，最好能将它们捆绑操作，当选定好要检验的无关参数后，甚至可以将将模型中的其它参数去除，以简化讨论，也就是说假设条件变成\(a_1=\cdots=a_p=0\)。但这个多变量的假设很难建立之前单变量的枢轴变量，我们需要另外找一个变量作为度量的对象。在鉴别“有效、无效”元素的问题中，注意“有效”的元素的典型特征，就是使得残差平方和变小，或者说使得\(\hat{\sigma}^2\)尽量小。这便是我们要找的“值”，具体来说，就是要度量\(\hat{\sigma}^2\)和\(\sigma^2\)之间的差别。
但由于\(\sigma^2\)未知，必须找统计量来替代它，在假设条件下，自然是用\(S_Y^2\)。当直接用\(\hat{\sigma}^2\)和\(S_Y^2\)难产生好的枢轴变量，原因主要是系数的影响，这时我们自然想到直接比较残差平方和。为此记\(R_1=\delta^T\delta\)，并记假设条件下的残差平方和为\(R_2\)。为了讨论方便，这里把式（14）稍作修改，就是先作出估计\(\alpha_0=\bar{Y}\)，然后用\(Y_i\)取代\(Y_i-\alpha_0\)重新建模，随之\(\gamma,A\)中的第一列也去除。
但新的模型仍然能得到式（17）的估计式，以及残差向量\(\delta=\beta-A^T\alpha\)。这个模型下\(R_1,R_2\)如式（26）所示，不难发现\(R_2-R_1=\beta^TB\beta\)，而已知\(B^2=B\)，所以\(R_2-R_1\)也是一个半正定二次型。再次使用柯赫伦定理可有式（27）左成立，并且\(R_2-R_1\)与\(R_1\)互相独立，这等价于与\(\hat{\sigma}^2\)互相独立，所以得到式（27）右的枢轴变量。注意到\(R_2\geqslant R_1\)，故检验否定的条件应当是\(F>C\)。
\[R_1=\delta^T\delta=\beta^T(I_n-B)\beta;\;\;R_2=\beta_T\beta\tag{26}\]
\[\dfrac{R_2-R_1}{\sigma^2}\sim\chi_p^2;\;\;\dfrac{R_2-R_1}{r\hat{\sigma}^2}\sim F_{p,n-p-1}\tag{27}\]
