【数理统计基础】 05 - 回归分析

合集 - 随机数学(11)

1.【初等概率论】 01 - 不确定中的确定性2017-02-072.【初等概率论】 02 - 条件概率和随机变量2017-02-073.【初等概率论】 03 - 常见分布和问题举例2017-02-074.【初等概率论】 04 - 数字特征2017-02-075.【初等概率论】 05 - 极限定理和正态分布2017-02-076.【数理统计基础】 01 - 关于预测的方法学2017-05-227.【数理统计基础】 02 - 统计量和三大分布2017-05-228.【数理统计基础】 03 - 参数估计2017-05-229.【数理统计基础】 04 - 假设检验2017-05-22

10.【数理统计基础】 05 - 回归分析2017-05-22

11.【数理统计基础】 06 - 相关分析和方差分析2017-05-24

收起

　　参数估计和假设检验是数理统计的两个基础问题，它们不光运用于常见的分布，还会出现在各种问题的讨论中。本篇开始研究另一大类问题，就是讨论多个随机变量之间的关系。现实生活中的数据杂乱无章，够挖掘出各种变量之间的关系非常有用，它可以预估变量的走势，能帮助分析状态的根源。关系分析的着手点可以有很多，我们从最简单直观的开始，逐步展开讨论。

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$独立的。在应用场景，一般给定回归函数一个含参表达式（比如后面的线性回归），这样的问题称为参数回归问题，否则叫非参数回归问题。

$$\begin{array}{}\text{(1)}& Y=f(X_1,X_2,\cdots ,X_p)+e,E(e)=0,0<D(e)<\mathrm{\infty }\end{array}$$

　　在微积分中我们知道，多元函数$f(x_1,\cdots ,x_p)$一般可以在任何点进行泰勒展开，其中最简单的就是线性展开。线性关系由于其形式简单，以及在局部能很好地逼近函数，在数学的各分支都被重点讨论。在回归分析中，这样的模型便称为线性回归，这里先从最简单的一元线性回归讨论起。

　　一元线性回归的模型是式（2）左，在提取出线性关系$b_0+b_{1}X$后，$Y$的剩余因素或随机性就都落在随机变量$e$上。所以从另一个角度看，回归分析是要找出随机变量的“确定”部分和随机部分，这种分解更能帮助分析随机现象。自然地，分析是基于$n$个样本点$(X_i,Y_i)$，其中$X_i$可能也是随机产生的，但在这个模型里一律看做定量。还要注意，这时每个$Y_i$是$e_i$的一个位移，它不再与$Y$同分布。

$$\begin{array}{}\text{(2)}& Y=b_0+b_{1}X+e;Y_i=b_0+b_1{X}_i+e_i=a_0+a_1(X_i-\overline{X})+e_i\end{array}$$

1.2 系数点估计

　　每一次试验相互独立，因此得到$n$个独立变量（式（2）右），其中$b_0,b_1$是待定系数。为了方便计算和讨论，一般还会把式（2）右的线性部分“中心化”。问题等价于讨论$a_0,a_1$的值，但要注意这里$\overline{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$的斜率，单点的斜率是${\displaystyle \frac{Y_i-\overline{Y}}{X_i-\overline{X}}}$，将分子分母同时乘以$X_i-\overline{X}$并相加，化简后便得到${\alpha }_1$，故它也是斜率的合理估计。

$$\begin{array}{}\text{(3)}& \sum _{i=1}^n(Y_i-b_0-b_1{X}_i{)}^2\end{array}$$

$$\begin{array}{}\text{(4)}& {\alpha }_0=\overline{Y};{\alpha }_1=\sum _{i=1}^n{\displaystyle \frac{X_i-\overline{X}}{S^2}}{Y}_i,S^2=\sum _{i=1}^n(X_i-\overline{X}{)}^2\end{array}$$

　　另外还要注意，式（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估计。

$$\begin{array}{}\text{(5)}& E(Y_i)=a_0+a_1(X_i-\overline{X})\Rightarrow E({\alpha }_0)=a_0,E({\alpha }_1)=a_1\end{array}$$

$$\begin{array}{}\text{(6)}& D(Y_i)={\sigma }^2\Rightarrow D({\alpha }_0)={\displaystyle \frac{{\sigma }^2}{n}},D({\alpha }_1)={\displaystyle \frac{{\sigma }^2}{S^2}}\end{array}$$

　　还有一点，把${\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$的无偏估计。

$$\begin{array}{}\text{(7)}& {\hat{\sigma }}^2={\displaystyle \frac{1}{n-2}}\sum _{i=1}^n{\delta }_i^2,{\delta }_i=Y_i-{\alpha }_0-{\alpha }_1(X_i-\overline{X})\end{array}$$

　　具体计算步骤参考教材（或自行证明），结果是得到式（8），这样就不难得到是（7）了。当然在实际计算时，可以直接展开得到式（9），然后利用现成的$X_i,Y_i,{\alpha }_j$来加速计算。而且从式（9）中还能得到更有用的结论，注意其中的后两项$n{\overline{Y}}^2=Z_1^2$和$S^2{\alpha }_1^2=Z_2^2$，$Z_1,Z_2$都是$Y_i$的线性函数，且系数向量是两个相互正交的标准化向量。

$$\begin{array}{}\text{(8)}& \sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n(e_i-\overline{e}{)}^2+{\displaystyle \frac{1}{S^2}}{(\sum _{i=1}^n(X_i-\overline{X})e_i)}^2\end{array}$$

$$\begin{array}{}\text{(9)}& \sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n{Y}_i^2-n{\overline{Y}}^2-S^2{\alpha }_1^2\end{array}$$

　　当$e$是正态分布时，$Y_i$也是正态分布，利用正交变换的性质，易知式（9）等于$Z_3^2+\cdots +Z_n^2$，其中$Z_j\sim N(0,{\sigma }^2)$，这便容易有式（10）的结论。关于残差，还有两点需要注意，式（8）如果很大或者残差体现出某些规律性，则说明线性模型不太合适，或还有重要因素没有被提取出来。

$$\begin{array}{}\text{(10)}& e\sim N(0,{\sigma }^2)\Rightarrow \sum _{i=1}^n{\displaystyle \frac{{\delta }_i^2}{{\sigma }^2}}\sim {\chi }_{n-2}^2\end{array}$$

1.4 区间估计

　　有了点估计，便可以做区间估计，为了能使用枢轴函数，这里还是假定$e$为正态分布。首先由公式（5）（6）可知${\alpha }_0,{\alpha }_1$满足式（11）的分布，当$\sigma$已知时，枢轴函数很容易得到。当$\sigma$未知时，由刚才的讨论知${\alpha }_0,{\alpha }_1$与${\hat{\sigma }}^2$是相互独立的，这样便能用$\hat{\sigma }$替代$\sigma$，得到式（12）的枢轴变量。

$$\begin{array}{}\text{(11)}& {\alpha }_0\sim N(a_0,{\displaystyle \frac{{\sigma }^2}{n}});{\alpha }_1\sim N(a_1,{\displaystyle \frac{{\sigma }^2}{S^2}})\end{array}$$

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\sqrt{n}({\alpha }_0-a_0)}{\hat{\sigma }}}\sim t_{n-2};{\displaystyle \frac{S({\alpha }_1-a_1)}{\hat{\sigma }}}\sim t_{n-2}\end{array}$$

　　线性回归的目的自然是为了进行预测，但在仅知道样本点且把$X_i$看成定量的情况下，其实是无法估计最初式（2）左中的$b_0,b_1$的。因此要注意，在用$y=a_0+a_1(x-\overline{X})$预估$Y$时，我们不光丢失了误差$e$，还丢失了$X$非连续得来的误差。前者通过合理建模来降低误差，后者则只能通过增加$X_i$的数量和密度来降低误差。

　　这一点容易通过估计值$y$的方差看出（式（13））。首先在样本数不变的情况下，$x$离样本中心$\overline{X}$越近方差越小，这个结论符合直觉，样本离预测点越近精度越高。另一方面，样本数越大方差也越小，这个很好理解。结合这两方面，当$n$足够大且$x$离样本中心足够近，估计的方差就可以任意小。

$$\begin{array}{}\text{(13)}& D(y)=({\displaystyle \frac{1}{n}}+{\displaystyle \frac{(x-\overline{X}{)}^2}{S^2}}){\sigma }^2\end{array}$$

2. 多元线性回归

2.1 系数估计

　　现实中的因变量可能受多个因素的影响，这些因素可能有主次之分，也可能是联合作用。无论如何，为了对因变量进行更加深入细致的分析，必须加入更多的自变量进行分析。另外同样的道理，多元函数在局部都可以用线性函数很好地近似，因此我们也可以建立式（14）中的模型和中心化样本表达式。为表达方便，本段下面就直接把$X_{ki}-{\overline{X}}_k$记作$X_{ki}$。

$$\begin{array}{}\text{(14)}& Y=b_0+\sum _{k=1}^p{b}_k{X}_k+e;Y_i=a_0+\sum _{k=1}^p{a}_k(X_{ki}-{\overline{X}}_k)+e_i\end{array}$$

　　多元模型的讨论内容和方法与一元的差别不大，但直接的讨论会很繁琐，必须借助于线性代数的工具，请注意前后对比。为讨论方便，首先规定式（15）的简写，并记$\gamma$的点估计为$\alpha$。然后定义列向量的期望$E(\alpha )=[E({\alpha }_i)]$，以及协方差$\text{Cov}(\alpha ,\beta )=[\text{Cov}({\alpha }_i,{\beta }_j)]$，且不难验证有式（16）成立。其实利用算子理论证明会很简单，但光凭形式化的假设，也不难完成证明，请独立尝试。

$$\begin{array}{}\text{(15)}& \beta =\left[\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right],\gamma =\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_p\end{array}\right],A=\left[\begin{array}{ccc}1& \cdots & 1\\ X_{11}& \cdots & X_{1n}\\ ⋮& \ddots & ⋮\\ X_{p1}& \cdots & X_{pn}\end{array}\right]\end{array}$$

$$\begin{array}{}\text{(16)}& E(A\alpha )=AE(\alpha );\text{Cov}(A\alpha ,B\beta )=A\text{Cov}(\alpha ,\beta )B^T\end{array}$$

　　有了矩阵的定义，就可以直接利用线性代数中最小二乘的结论，得到（17）左的正则方程，以及式（17）的$\gamma$最小二乘解。式（18）推导出${\alpha }_i$是$a_i$的无偏估计，且都是$Y_i$的线性函数，继而还可以得到式（19）的协方差公式。注意到$A$中除第一列外，每行的和都是$0$，故$L$及$L^{-1}$都具有形式$\left[\begin{array}{cc}1& 0\\ 0& B_{p\times p}\end{array}\right]$。这说明${\alpha }_0=\overline{Y}$与其它${\alpha }_i$互不相关，这与一元的情况是一致的。

$$\begin{array}{}\text{(17)}& L\alpha =A\beta \Rightarrow \alpha =L^{-1}A\beta ,(L=A{A}^T)\end{array}$$

$$\begin{array}{}\text{(18)}& E(\alpha )=L^{-1}AE(\beta )=L^{-1}A{A}^T\gamma =\gamma \end{array}$$

$$\begin{array}{}\text{(19)}& \text{Cov}(\alpha ,\alpha )=L^{-1}A\text{Cov}(\beta ,\beta )A^T{L}^{-1}={\sigma }^2{L}^{-1}\end{array}$$

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$（主对角线之和，请参考线性代数）。

$$\begin{array}{}\text{(20)}& \text{Cov}(\hat{\alpha },\delta )=L^{-1}A\text{Cov}(\beta ,\beta )(I_n-A^T{L}^{-1}A)=0\end{array}$$

$$\begin{array}{}\text{(21)}& \text{Cov}(\delta ,\delta )=(I_n-B)\text{Cov}(\beta ,\beta )(I_n-B)={\sigma }^2(I_n-B),(B=A^T{L}^{-1}A)\end{array}$$

　　为了估计${\sigma }^2$，自然想到讨论残差平方和$\sum _{i=1}^n{\delta }_i^2$。式（22）计算了它的期望值，这样就可以用式（23）来无偏估计${\sigma }^2$。

$$\begin{array}{}\text{(22)}& E(\sum _{i=1}^n{\delta }_i^2)=\sum _{i=1}^{n}D({\delta }_i)=tr(\text{Cov}(\delta ,\delta ))={\sigma }^2(n-p-1)\end{array}$$

$$\begin{array}{}\text{(23)}& {\hat{\sigma }}^2={\displaystyle \frac{1}{n-p-1}}\sum _{i=1}^n{\delta }_i^2\end{array}$$

　　残差平方和${\delta }^T\delta$是一个半正定二次型，展开整理后有式（24）成立，它满足柯赫伦定理的条件。故假定$e$为正态分布的情况下，有式（25）左成立。另外由于${\alpha }_i$与${\delta }_j$互不相关，则${\alpha }_i$与$\hat{\sigma }$也不相关，正态分布下它们还是相互独立的，这就得到式（25）右的枢轴变量。

$$\begin{array}{}\text{(24)}& {\beta }^T\beta ={\beta }^{T}B\beta +{\delta }^T\delta \end{array}$$

$$\begin{array}{}\text{(25)}& \sum _{i=1}^n{\displaystyle \frac{{\delta }_i^2}{{\sigma }^2}}\sim {\chi }_{n-p-1};{\displaystyle \frac{{\alpha }_i-a_i}{\sqrt{L_{ii}^{-1}}\hat{\sigma }}}\sim t_{n-p+1}\end{array}$$

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=\overline{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 }^{T}B\beta$，而已知$B^2=B$，所以$R_2-R_1$也是一个半正定二次型。再次使用柯赫伦定理可有式（27）左成立，并且$R_2-R_1$与$R_1$互相独立，这等价于与${\hat{\sigma }}^2$互相独立，所以得到式（27）右的枢轴变量。注意到$R_2⩾R_1$，故检验否定的条件应当是$F>C$。

$$\begin{array}{}\text{(26)}& R_1={\delta }^T\delta ={\beta }^T(I_n-B)\beta ;R_2={\beta }_T\beta \end{array}$$

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{R_2-R_1}{{\sigma }^2}}\sim {\chi }_p^2;{\displaystyle \frac{R_2-R_1}{r{\hat{\sigma }}^2}}\sim F_{p,n-p-1}\end{array}$$