回归分析第2讲——线性回归及其拟合

一、回归分析的概念

万事万物皆有联系。此消彼长是常态，同进同退不奇怪。如果我们想知道两个变量$X,Y$之间是否有关联，我们要怎么做呢？

很直接的想法就是测量多组$X$和$Y$的观测值，之后绘制散点图，看看随着$X$的变化，$Y$有没有发生变化。如果$X$增大，$Y$也增大，说明两个变量是有关联的，并且是正向的关联；如果$X$增大，$Y$却减小，说明两个变量是反向的关联。

再代数一点，假设在一次实验中，我们获得了两个变量$X,Y$的$n$组观测值$(x_i,y_i)$，其中$i=1,\cdots,n$。我们希望通过这$n$组观测值来推断两个变量的关系，也就是找到这样一个关系$f$，满足：

$$Y=f(\theta;X)$$

其中$\theta$是一个确定但未知的参数。

像这样的通过变量对应的大量观测值来推断关系$f$的过程，在统计上就称为回归(Regression)。

”回归“一词最早是由弗朗西斯·高尔顿提出的，用于描述他在优生学方面的发现。不过，现代的”回归“一词早已不等同于高尔顿所提出的”回归“的含义了。

照理来说，对于某一个确定的研究对象，这个关系$f$应当是一种未知但确定的规则。然而，研究对象的观测值是存在波动的，这可能是因为还有别的因素在影响观测值，或者就是一些不清楚、不知道的误差在作祟。因此，我们推断得到的$f$只能是近似的、拟真的。

为了表示不可控制的波动，在回归分析中，我们认为两个变量的基本关系应当形如：

$$Y=f(X_1,\cdots,X_p)+e$$

其中$Y$代表观测值，$X_p$代表第$p$个影响因素，$f(X_1,\cdots,X_p)$代表回归的拟合值，$e$代表观测值和拟合值之间的误差。

既然$f$是近似的，那么总是需要一种公认的标准来进行界定，否则百家争鸣，无法开展研究。那么如何设立标准呢？一般地，在回归分析领域，我们认为最优回归应当满足：

$$\min \sum_{i=1}^ne_i^2=\min \sum_{i=1}^n[y_i-\hat{y_i}]^2$$

其中$y_i$代表第$i$个样本的观测值，$\hat y_i$代表第$i$个样本的拟合值/估计值。

也就是说，我们需要尽可能使得拟合值与观测值接近，偏差最小。并且我们还希望控制$E(e)=0$，也就是误差的数学期望是0，即误差不成为影响观测值的主导因素，这样的误差又被称为随机误差。

另外，在保证拟合度的情况下，我们选择的曲线应当光滑；如果有经验曲线，应当考虑在内。

总的来说，寻找关系$f$的过程，就是寻找：

$$f(Y|X_1,\cdots,X_p)=E(Y|X_1,\cdots,X_p)$$

接下来的问题就是，如何寻找关系$f$？

简单来说，寻找关系$f$的方法无非两大类：

1. 假定关系的带参形式$f(X,\theta)$，然后通过观测值拟合参数$\theta$（传统的回归分析方法）
2. 无法假定关系的形式，采用学习方法实现端到端的映射（机器学习、深度学习、大模型）

回归分析专题主要探讨第一类方法，也就是假定关系的形式，再拟合参数；而第二类方法会放在其他专题当中单独学习。

既然我们要假定关系的形式，就可以把回归分为几类模型：

• 线性模型
  • 线性回归模型
  • 多项式回归模型
  • 方差分析模型
• 广义线性模型
  • 逻辑回归模型
  • 泊松回归、负二项回归模型
  • β-回归、γ-回归模型
• 拓展模型
  • 稳健回归模型
  • 正则化回归模型
  • 线性混合效应模型
  • 广义线性混合效应模型
  • 广义加性模型
• 生物信息学应用模型
  • Cox回归模型与生存分析
  • 全基因组关联分析（GWAS）

这一讲我们先从最简单的模型讲起：线性回归。

二、线性回归的原理

2.1 线性回归的概念

线性回归模型是回归分析中最简单、最基本的回归方法，但是其重要性不可小觑，目前很多回归方法都是基于线性回归来改造的。

笔者观点：统计学的研究应当服从奥卡姆剃刀原则，如果能用更简单的模型实现近似好的效果，就不应该使用复杂模型。不过，在实际数据中，最简单的这一种线性回归也不太好用了，但是学习基础原理可以帮我们更好地理解更复杂的模型。

从几何上理解，线性回归就是认为观测值的分布是线性变化的。例如下图这组数据，描述了学生每周学习时间和学习成绩的关系，可以看到这两个变量之间大致是服从线性变化的，如果使用曲线来描述就不太能表达数据的特征。

从代数上理解，线性回归是指假定关系$f$的参数是线性的，满足$E(e)=0$，且$Var(e)=\sigma^2$越小越好。一般地，线性回归模型可以表示为：

$$y=\alpha+x_1\beta_1+\cdots+x_p\beta_p+e$$

其中$x_1,\cdots,x_p$是自变量，$\beta_1,\cdots,\beta_p$是回归系数，$\alpha$是截距项，$e$是随机误差，$y$是响应变量。

在这一小节我们不考虑截距项，即认为：

$$y=x_1\beta_1+\cdots+x_p\beta_p+e$$

如何理解”参数是线性的“？

参数是线性的，意味着$\beta_1,\cdots,\beta_p$不能出现高次项，例如不能出现$x_1\beta_1^2$，这样就不是线性回归了。

但是，形如$y=x_1\beta_1+x_2^2\beta_2+\beta_3\log x_3$这样的也可以称之为线性回归，虽然出现了非一次项，可以拟合曲线，但是参数仍然是线性的。这类模型又可以叫做多项式回归，求解方式与一般线性回归没有差别，只是自变量发生了变化。

如果使用多组观测值来拟合这个模型，就要进行一点变化。假设我们一共观测了$n$次（例如上图，$n$个学生），那么就可以写作：

$$y_i=x_{i1}\beta_1+\cdots+x_{ip}\beta_{p}+e_i,i=1,\cdots,n$$

如果写成矩阵形式就是：

$$\begin{pmatrix} y_1\\ \vdots\\ y_n \end{pmatrix}= \begin{pmatrix} x_{11}&\cdots&x_{1p}\\ \vdots&\ddots&\vdots\\ x_{n1}&\cdots&x_{np} \end{pmatrix} \begin{pmatrix} \beta_1\\ \vdots\\ \beta_p \end{pmatrix}+ \begin{pmatrix} e_1\\ \vdots\\ e_n \end{pmatrix}$$

即：

$$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}\ \ \ E(e)=0$$

我们称矩阵$X_{n\times p}$为设计矩阵(design matrix)，$X\beta$又被称为线性预测子。一般认为$rank(X)=p$，即列满秩。如果$rank(X)<p$，说明自变量存在冗余。

既然要拟合模型，那么我们首先要确定模型的参数量。在线性回归中，我们要拟合的参数有：

• $\beta_1,\cdots,\beta_p$，称为回归系数，用于描述每个自变量的贡献，估计值为$\hat\beta_1,\cdots,\hat\beta_p$，共$p$个
• $Cov(e_i,e_j)$，称为误差的协方差阵，用于描述误差的波动情况，估计值为$Cov(\hat e_i,\hat e_j)$，共$(n^2+n)/2$个

因此，线性回归模型的参数总量是$p+(n^2+n)/2$个。所以你可以看到，参数总量是很大的，尤其当观测次数$n$很大时，估计会变得非常困难。有没有简化方法呢？

当然有，我们会假定$n$次观测是独立等方差的，此时协方差阵可以简化为$Cov(e_i,e_j)=\sigma^2I_n$（$I_n$为单位矩阵），那么参数总量减少为$p+1$个，计算量大大降低。这个假定被称为高斯-马尔可夫假定(Gaussian-Markov, GM)，这也是为什么很多相关书籍会假定线性模型的每一次观测是独立等方差的。

2.2 线性回归的参数估计

对线性回归来说，最重要的参数就是回归系数，这是我们参数估计的重点。对回归系数的估计采用广为人知的最小二乘法(Ordinary Least Square, OLS)。最小二乘法的核心思想就是使得观测值与估计值的误差尽可能地小，从而完成对回归系数地估计。这与我们之前提到的”最优回归“不谋而合。

我们先从独立等方差的情形开始求解。对于线性模型$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}$，$Ee=0$，$Cov(e)=\sigma^2 I_n$来说，其最小二乘的估计函数为：

$$Q(\beta)= ||e||^2= ||Y-X\beta|| = (Y-X\beta)^T(Y-X\beta)$$

优化目标为：

$$\min_\beta Q(\beta)=\min_\beta\sum_{i=1}^ne_i^2=\min_\beta\sum_{i=1}^n(y_i-x_i^T\beta)^2$$

其中$x_i=(x_{i1},\cdots,x_{ip})^T$。

要求解这个优化问题，只需要求解$\frac{\partial Q(\beta)}{\partial \beta}=0$即可。最后解得：

$$X^TX\beta=X^TY$$

这个方程被称为正规方程(Normal Equation)，分量形式可以写作：

$$(\sum_{i=1}^nx_ix_i^T)\beta=\sum_{i=1}^nx_iy_i$$

求解过程

直接采用矩阵形式求解，首先可以证明

$$\begin{aligned}\frac{\partial C^T\beta}{\partial\beta}&=C\\ \frac{\partial\beta^TA\beta}{\partial\beta}&=(A+A^T)\beta\end{aligned}$$

于是可以求解：

$$\begin{aligned}\frac{\partial Q(\beta)}{\partial \beta}&=\frac{\partial}{\partial \beta}(Y^TY-2Y^TX\beta+\beta^TX^TX\beta)\\ &=-2X^TY+2X^TX\beta=0\end{aligned}$$

最后解得：$X^TX\beta=X^TY$，证毕。

得到正规方程后，我们还要考虑两个问题：

1. 正规方程是否有解？
2. 如果有解，有多少个解？

首先，正规方程一定有解，因为：

$$X^TY\in \mu(X^TX)=\left\{X^TX\beta|\beta\in\mathbb R^P\right\}$$

所以至少存在一个解为：

$$\hat\beta=(X^TX)^-X^TY$$

其次，有多少个解取决于设计矩阵的秩：

• 如果$rank(X)=p$，即列满秩，那么$\hat\beta$存在唯一解$\hat\beta=(X^TX)^-X^TY$；
• 如果$rank(X)<p$，说明存在变量冗余，$\hat\beta$是不可估计的，不存在对$\hat\beta$的无偏估计。

回归系数解决完以后，那么方差参数$\sigma^2$也就很好估计了。由于$\hat e=Y-X\hat\beta=(I_n-P_X)Y$，$E\hat e=0$，$Cov(\hat e)=\sigma^2(I_n-P_X)$，设$rank(X)=r$，则对方差参数的估计为：

$$\hat\sigma^2=\frac{\hat e^T\hat e}{n-r}=\frac{Y^T(I_n-P_X)Y}{n-r}$$

另外，如果误差还服从正态分布，即$e\sim N_r(0,\sigma^2I_n)$，那么这样的误差又被称为高斯噪声(Gaussian Noise)，还有新的结论：

1. $\hat \beta$不仅是$\beta$的最小二乘估计，还是极大似然估计，且$\hat\beta\sim N(\beta,\sigma^2(X^TX)^-)$；
2. $(n-r)\hat\sigma^2/n$是$\sigma^2$的极大似然估计，并且$(n-r)\hat\sigma^2/\sigma^2\sim \chi_{n-r}^2$；
3. $\hat\beta$与$\hat\sigma^2$相互独立。

三、有约束的线性回归

前面我们讨论了一般线性回归的参数估计，接下来我们加一点难度，探讨有约束的线性回归模型。

有约束，是指对回归系数$\beta$有约束条件$L\beta=d$，其中$d$为某常数。有约束的线性回归模型可以写作：

$$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}\ \ \ L_{q\times p}\beta=d$$

其中$q$表示约束条件的个数。

有约束线性回归在应用当中有很大用处，一方面可以限制回归系数的总和，另一方面可以在对应系数前设置0来控制参数个数，例如取$L=(0,1,0,\cdots,0)$就表示只取$\beta_2$参与回归拟合。

3.1 有约束线性回归的参数估计

为了不失一般性以及简洁，我们假定模型为：

$$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}\ \ \ L_{q\times p}\beta=0$$

且满足$rank(L_{q\times p})=q$，$p\gt q$。当然，$L\beta=0$一定是相容的条件才行。

如果$d≠0$，则取$\beta_0$使得$L\beta_0=d$，令$\widetilde Y=Y-X\beta$，$\widetilde\beta=\beta-\beta_0$，线性模型$\widetilde Y=X\widetilde\beta+e$的约束条件又变为$L\widetilde\beta=0$。所以这里直接假定$d=0$来讨论。

类似于最小二乘法，在约束条件下，回归系数$\hat\beta_L$的估计值为：

$$\hat\beta_L=\arg \min_{L\beta=0}||Y-X\beta||^2$$

直接对$\beta$求导是不够的，因为存在约束条件。根据优化理论，我们需要引入拉格朗日乘子$\lambda_{q\times1}$：

$$L(\beta,\lambda)= ||Y-X\beta||^2+2\lambda^TL\beta$$

分别对$\beta,\lambda$求导得到：

$$\begin{aligned} \frac{\partial L(\beta,\lambda)}{\partial\beta}&=-2X^TY+2X^TX\beta+2L^T\lambda=0\\ \frac{\partial L(\beta,\lambda)}{\partial\lambda}&=2L\beta=0 \end{aligned}$$

最终得到：

$$\begin{pmatrix} X^TX&L^T\\ L&0 \end{pmatrix} \begin{pmatrix} \beta\\ \lambda \end{pmatrix}= \begin{pmatrix} X^TY\\0 \end{pmatrix}$$

可以证明，这个方程一定有解。

证明：

先给出两个引理，这两个引理的证明请看附录。

引理1：设$S=\{A_{n\times m}x|B_{k\times m}x=0,x\in\mathbb R^m\}$，那么$S$是线性子空间且$S$的维数是$dim S=rank\begin{pmatrix}A\\ B\end{pmatrix}-rank(B)$
引理2：设矩阵$V_{p\times p}\ge0$，$A_{p\times q}$，那么有：

• $\mu(A)\bigcap\mu(VA^{\bot})=\{0\}$
• $\mu(V\vdots A)=\mu(VA^{\bot}\vdots A)$

基于这两个引理，我们就可以开始证明了。根据$L\beta=0$先解出

$$\beta=(I_p-L^-L)Z,\forall Z_{p\times 1}$$

带入第一个方程有：

$$X^TX(I_p-L^-L)Z+L^T\lambda=X^TY$$

由于$I_p-L^-L=(L^T)^\bot$，根据引理2，有：

$$X^TY\in\mu(X^TX)\subset\mu(X^TX\vdots L^T)=\mu[X^TX(I_p-L^-L)\vdots L^T]$$

故方程有解。

既然方程一定有解，那么有多少个解呢？这就要使用可估的概念了。详情请跳转下一小节3.2，这里我们先把至少的那一个解尝试写出来。

由于正规方程非常复杂，如果要显式地写出$\hat\beta_L$的表达式，那我们只能设

$$\begin{pmatrix} X^TX&L^T\\ L&0 \end{pmatrix}^-= \begin{pmatrix} G_{11}&G_{12}\\ G_{21}&G_{22} \end{pmatrix}$$

因此$\hat\beta_L=G_{11}X^TY$，且$Var(\hat\beta_L)=\sigma^2G_{11}$。

当然，在一些条件下，$\hat\beta_L$有更简单的形式。仍假定$rank(L_{q\times p})=q$，那么：

• 若$\mu(L^T)\subset\mu(X^T)$，那么$\hat\beta_L=\hat\beta-(X^TX)^-L^T[L(X^TX)^-L^T]^{-1}L\hat\beta$，其中$\hat\beta=(X^TX)^-(X^TY)$。
• 若$rank(X_{n\times p})=p$，那么$\hat\beta_L=\hat\beta-(X^TX)^{-1}L^T[L(X^TX)^{-1}L^T]^{-1}L\hat\beta$，其中$\hat\beta=(X^TX)^-(X^TY)$。

最后，有约束线性回归的方差$\sigma_L^2$的估计值为：

$$\hat\sigma_L^2=\frac{||Y-X\hat\beta_L||^2}{n-s}=\frac{Y^T(I_n-P_M)Y}{n-s}$$

其中$M=X(L^T)^{\bot}$，$s=rank\begin{pmatrix}S\\ L\end{pmatrix}-rank(L)$。

同理，如果误差还服从正态分布，即$e\sim N_r(0,\sigma^2I_n)$，那么一定有：

1. $\hat\beta_L\sim N(\beta,\sigma^2G_{11})$；
2. $(n-s)\hat\sigma_L^2/\sigma^2\sim \chi_{n-s}^2$，其中$s=rank\begin{pmatrix}S\\ L\end{pmatrix}-rank(L)$；
3. $\hat\beta_L$与$\hat\sigma_L^2$相互独立。

3.2 可估与条件可估

为了探讨正规方程解的个数，我们首先要介绍可估函数和可估原理。

给定$c$，使得$c^T\beta$是关于$\beta$的线性函数。如果此时存在关于$Y$的线性函数$a^TY$使得对任意$\beta$都有：

$$Ea^TY=c^T\beta$$

则称$c^T\beta$是一个可估函数。

可估函数的意义就是把不可估的一组参数$\beta$转化为可估的形式。例如，两个物体的重量$\beta_1,\beta_2$未知，如果把它们同时放在天平上称$n$次，设第$i$次结果为$y_i$，那么可以构成模型：

$$y_i=\beta_1+\beta_2+e_i,i=1,\cdots,n$$

假设该模型符合GM条件，令$c=(1,1)^T$，$\beta=(\beta_1,\beta_2)^T$，那么对于每一个单独的$\beta_1,\beta_2$来说都是不可估的，但是$c^T\beta=\beta_1+\beta_2$就是可估的，$y_1$就是一个无偏估计。

对于一个可估函数，下面三条判定定理都是等价的：

• $c^T\beta$可估；
• $X\beta_1=X\beta_2\Rightarrow c^T\beta_1=c^T\beta_2$；
• $c\in \mu(X^T)$

可估函数有如下性质：

1. 所有的$c^T\beta$可估$\Leftrightarrow rank(X)=p$；最多有$rank(X)$个可估函数。
2. 如果$c_1^T\beta,c_2^T\beta$都可估，那么线性组合$\lambda_1c_1^T\beta+\lambda_2c_2^T\beta$可估；如果$c_1,c_2$还是线性无关的，那么$c_1^T\beta,c_2^T\beta$也是线性无关的。
3. 如果$c^T\beta$可估，那么$c^T\hat\beta$是唯一的，与广义逆$(X^TX)^-$的取值无关。

接下来介绍可估原理。如果$c^T\beta$可估，那么对于$\forall b\in \mu(X)^{\bot}$，$c^T\beta$的所有无偏估计为$(a+b)^TY$，而其中方差最小的估计则被称为最优线性无偏估计(Best Linear Unbiased Estimator, BLUE)，或称为高斯-马尔可夫估计(GM估计)。进一步地，高斯-马尔可夫定理告诉我们，如果$c^T\beta$可估，那么$c^T\hat\beta$是唯一的BLUE。

（不想阅读可跳过）

证明：为什么$c^T\hat\beta$是唯一的BLUE？

任取$a^TY$是$c^T\beta$的一个无偏估计，那么有：

$$\begin{aligned}\forall\beta,c^T\beta&=Ea^TY=a^TX\beta\\ \Rightarrow c^T&=a^TX\end{aligned}$$

分别求方差：

$$\begin{aligned}Var(a^TY)&=a^TVar(Y)a=\sigma^2a^Ta\\ Var(c^T\hat\beta)&=Var[a^TX(X^TX)^-X^TY]\\ &=Var(a^TP_XY)\\ &=a^TP_XVar(Y)P_Xa\\ &=\sigma^2a^TP_Xa\end{aligned}$$

作差得：

$$Var(a^TY)-Var(c^T\beta)=\sigma^2a^T(I_n-P_X)a\ge0$$

当$a^T(I_n-P_X)a=0$时等号成立，此时$a=P_Xa$，于是有：

$$\begin{aligned}a^TY&=a^TP_XY\\ &=a^TX(X^TX)^-X^TY\\ &=c^T\hat\beta\end{aligned}$$

于是原命题得证。

所以，对于一般线性回归来说，如果$rank(X)<p$不满秩，就说明存在一些不可估的量，所以会有无穷多解；对于那些可估的量，最小二乘解是唯一的解。

那么拓展到有约束线性回归中，就可以衍生出条件可估的概念。条件可估是指存在$a$使得所有满足$L\beta=0$的$\beta$都有：

$$Ea^TY=c^T\beta$$

可以看到，条件可估添加了约束条件$L\beta=0$，这样使得原本的限制变得宽松，可估的函数数量增加，即满足：

$$c\in \mu(X^T\vdots L^T)$$

回忆：一般可估函数满足$c\in \mu(X^T)$

（不想阅读可跳过）

证明：为什么加入了约束条件反而使得可估函数数量增加？

根据条件可估的定义，一定有：

$$\forall\beta\in\{\beta|L\beta=0\},Ea^TY=c^T\beta$$

那么可以等价为：

$$\begin{aligned}&\Leftrightarrow a^TX\beta=c^T\beta\\ &\Leftrightarrow(X^Ta-c)^T\beta=0\\ &\Leftrightarrow(X^Ta-c)^T(I_p-L^-L)Z=0\\ &\Leftrightarrow(X^Ta-c)^T(I_p-L^-L)=0\\ &\Leftrightarrow X^Ta-c\in\mu(L^T)\\ &\Leftrightarrow\exists b,x^Ta-c=L^Tb\\ &\Leftrightarrow c=x^Ta-L^Tb\\ &\Leftrightarrow c\in\mu(X^T\vdots L^T)\end{aligned}$$

于是原命题得证。

同理一般可估函数，对于$\forall c^T\beta$，$c^T\hat\beta_L$是唯一的，不依赖广义逆的选取，且是所有$c^T\beta$中的唯一BLUE。

在小节的最后，我们讨论一种特殊情况。既然条件$L\beta=0$使得可估函数的范围扩大，那么当$L$取何值时，可以使得对$\forall c$，$c^T\beta$都是条件可估的？

要达到这个目的，我们要考虑两个条件：

1. 添加约束后投影空间不变，即$rank\begin{pmatrix}X\\ L\end{pmatrix}-rank(L)=rank(X)$。
2. 添加约束后列满秩，即$rank\begin{pmatrix}X\\ L\end{pmatrix}=p$且$\mu(X^T)\bigcap\mu(L^T)=\{0\}$。

满足两个条件的约束$L\beta=0$被称为边界条件(side condition)。

在边界条件下，$\hat\beta_L$是正规方程的唯一解：

$$\hat\beta_L=(X^TX+L^TL)^{-1}X^TY$$

四、带截距项的线性回归

前面我们假定的一般线性回归是不考虑截距项的。但是，为了获得更一般的结论，我们往往要考虑模型带有截距项，我们记作$\beta_0$。

4.1 带截距项线性回归的参数估计

对于带截距项的模型，其分量形式为：

$$y_i=\beta_0x_0+\beta_1x_1+\cdots+\beta_{p-1}x_{p-1}+e_i$$

其中$x_0=E_n=(1,\cdots,1)^T$。

如果写成矩阵形式，令$X=(E_n\vdots \widetilde{X}_{n\times(p-1)})$，$\widetilde{X}_{n\times(p-1)}=(x_1^T,\cdots,x_{p-1}^T)^T$，$\beta=(\beta_0,\cdots,\beta_{p-1})^T=(\beta_0,\beta_I)^T$，那么：

$$Y_{n\times1}=(E_n\vdots\widetilde X_{n\times(p-1)})\begin{pmatrix}\beta_0\\ \beta_I\end{pmatrix}_{p\times1}+e_{n\times 1}$$

扩展开来写就是：

$$\begin{pmatrix}y_1\\ \vdots\\ y_n\end{pmatrix}= \begin{pmatrix}1&x_{11}&\cdots&x_{1,p-1}\\ 1&x_{21}&\cdots&x_{2,p-1}\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_{n1}&\cdots&x_{n,p-1}\end{pmatrix} \begin{pmatrix}\beta_0\\ \beta_1\\ \vdots\\ \beta_{p-1}\end{pmatrix}+ \begin{pmatrix}e_1\\ \vdots\\ e_n\end{pmatrix}$$

这样写就完全变成了我们熟悉的$Y=X\beta+e$，前面的结论都可以照常搬过来使用。这里我们集中总结一下相关的结论，这里假定满足GM条件且列满秩：

• $\hat\beta=(X^TX)^{-1}X^TY$
• $\hat\sigma^2=||Y-X\hat\beta||^2/(n-p)$
• $E(\hat\beta)=\beta$，$Var(\hat\beta)=\sigma^2(X^TX)^{-1}$，$E\hat\sigma^2=\sigma^2$
• 对于任意$c^T\beta$，$c^T\hat\beta$是其唯一的BLUE

4.2 带截距项的实际应用处理

4.1小节展示了一般性的结论，但是在实际应用中，我们更关心的是回归系数$\beta_I$而非截距项$\beta_0$，所以一般来说我们会这样书写：

$$Y=E_n\beta_0+\widetilde X\beta_I+e$$

另外，为了计算方便，实际应用中我们还会进行中心化处理和标准化处理。

首先介绍中心化。中心化就是把样本点的实际值转化为到样本中心点的差值，这样做的好处就是很方便地计算出截距项，从而减少计算量。可以证明，中心化之后参数的估计值与不中心化是一致的。

设$\overline x_j=\sum_{i=1}^nx_{ij}/n$，$1\le j\le p-1$，那么模型改写为：

$$y_i=\alpha+\beta_1(x_{i1}-\overline x_1)+\cdots+\beta_{p-1}(x_{i,p-1}-\overline x_{p-1})+e_i$$

可以计算得到：$\alpha=\beta_0+\overline x^T\beta_I$。上述分量形式写成矩阵形式为：

$$Y=\alpha E_n+\widetilde X_c\beta_I+e,Ee=0,Cov(e)=\sigma^2I_n$$

其中$\widetilde X_c=(I_n-\frac{E_nE_n^T}{n})\widetilde X$，$\widetilde X_c^TE_n=0$。

于是我们可以得到正规方程：

$$\begin{aligned} \begin{pmatrix}E_n^T\\ \widetilde X_c^T\end{pmatrix} \begin{pmatrix}E_n&\widetilde X_c\end{pmatrix} \begin{pmatrix}\alpha\\ \beta_I\end{pmatrix}&= \begin{pmatrix}E_n^T\\ \widetilde X_c^T\end{pmatrix}Y\\ \Rightarrow\begin{pmatrix}n&0\\ 0&\widetilde X_c^T\widetilde X_c\end{pmatrix} \begin{pmatrix}\alpha\\ \beta_I\end{pmatrix}&= \begin{pmatrix}n\overline y\\ \widetilde X_c^TY\end{pmatrix} \end{aligned}$$

于是解得：

$$\begin{aligned} \hat\alpha&=\overline y\\ \hat\beta_I&=(\widetilde X_c^T\widetilde X_c)^{-1}\widetilde X_c^TY \end{aligned}$$

进一步地：

$$Cov\begin{pmatrix}\hat\alpha\\ \hat\beta_I\end{pmatrix}=\sigma^2 \begin{pmatrix}\frac{1}{n}&0\\ 0&(\widetilde X_c^T\widetilde X_c)^{-1}\end{pmatrix}$$

通过中心化我们知道，截距项就可以通过观测值的均值来估计，方便快捷。

另一种是标准化。标准化是在中心化的基础上再除以标准差，其目的是消除量纲，这一点在分析协变量相关性上有重要应用。

令$s_j^2=\sum_{i=1}^n(x_{ij}-\overline x_j)^2$，$1\le j\le p-1$，$z_{ij}=\frac{x_{ij}-\overline x_j}{s_j}$，$Z=(z_{ij})_{n\times(p-1)}$，那么一定有$E_n^TZ=0$。

我们定义相关系数$R=(r_{ij})_{(p-1)\times(p-1)}=Z^TZ$，于是：

$$r_{ij}=\frac{\sum_{k=1}^n(x_{ki}-\overline x_i)(x_{kj}-\overline x_j)}{s_is_j}$$

相关系数可以分析协变量之间的关系，后续判断共线性问题中有重要应用。

五、特殊线性回归的处理

5.1 非独立观测与Aitken模型

前面的内容我们都在独立等方差的GM条件下探讨线性回归的参数估计，如果协方差阵具有一定结构应该怎么处理呢？

这里我们探讨观测不独立的情形，此时模型结构为

$$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}\ \ E(e)=0\ \ Cov(e)=\sigma^2\Sigma$$

其中$\Sigma>0$且已知。

这个模型被称为Aitken模型。Aitken模型仍然假定模型有公共的方差参数$\sigma^2$，但是观测之间不一定独立。

模型看似复杂，但是仍然可以转化为我们熟悉的模型。我们可以做如下变换：

$$\widetilde Y=\Sigma^{-\frac{1}{2}}Y,\widetilde X=\Sigma^{-\frac{1}{2}}X,\widetilde e=\Sigma^{-\frac{1}{2}}e$$

变换之后你会发现，$Cov(\widetilde e)=\sigma^2I_n$，就又满足GM假定了，因此上面两个小节的参数估计方法就可以套用在这个新模型上。这样我们就得到了关于回归系数$\beta$的广义解：

$$\beta^*=(X^T\Sigma^{-1}X)^-X^T\Sigma^{-1}Y$$

特别地，如果$\Sigma$是个对角阵，则说明观测之间是独立的，但是每次观测的权重不同，那我们称此时的$\beta^*$是加权最小二乘解。

5.2 奇异线性回归的处理

事实上，5.1小节中的Aitken模型还有一点缺陷，就是没有考虑$|\Sigma|=0$的情况（当然，$\Sigma$肯定不是负定的）。对于模型$Y_{n\times 1}=X_{n\times p}\beta_{p\times 1}+e_{n\times 1}$，$Ee=0$，$Cov(e)=\sigma^2\Sigma$，其中$\Sigma$已知但$|\Sigma|=0$，我们称之为奇异线性模型。

对于奇异线性模型，由于$|\Sigma|=0$，导致最小二乘法

$$Q(\beta)=(Y-X\beta)^T\Sigma^-(Y-X\beta)$$

是没有意义的。那么如何解决呢？

C.R.Rao认为，可以定义$T=\Sigma+XUX^T$，且$rank(T)=rank(\Sigma\vdots X)$，此时有：

$$Q(\beta)=(Y-X\beta)^TT^-(Y-X\beta)$$

进而可以求解。

再进一步，如果协方差阵含有参数，即$Cov(e)=\sigma^2\Sigma(\theta)$，其中$\Sigma(\theta)>0$且含未知参数$\theta$，那么处理起来就会更加麻烦。此时可以先对$\theta$进行估计得到$\hat\theta$，从而得到$\hat\Sigma(\hat\theta)$，最后获得回归系数的估计：

$$\hat\beta(\hat\theta)=\left[ X^T[\hat\Sigma(\hat\theta)]^{-1}X \right]X^T\left[ \hat\Sigma(\hat\theta) \right]^{-1}Y$$

附录

3.1节引理1：设$S=\{A_{n\times m}x|B_{k\times m}x=0,x\in\mathbb R^m\}$，那么$S$是线性子空间且$S$的维数是$dim S=rank\begin{pmatrix}A\\ B\end{pmatrix}-rank(B)$。证明：$S$是非空的，且

$$\begin{aligned} \forall Ax_1,Ax_2\in S&\Rightarrow Ax_1+Ax_2=A(x_1+x_2)\in S\\ \forall Ax\in S,c\in \mathbb R&\Rightarrow A(cx)\in S \end{aligned}$$

满足加法封闭和乘法封闭，所以属于线性子空间。

设$V$是$Bx=0$的解空间，定义线性变换$\varphi:V\rightarrow S$，则：

$$\begin{aligned} dim S &= dim V - ker\ \varphi\\ &=dim\{x|Bx=0,x\in\mathbb R^m\}-dim\{x|\begin{pmatrix}A\\ B\end{pmatrix}x=0,x\in\mathbb R^m\}\\ &=(m-rank(B))-(m-rank\begin{pmatrix}A\\ B\end{pmatrix})\\ &=rank\begin{pmatrix}A\\ B\end{pmatrix}-rank(B) \end{aligned}$$

原命题得证。

3.1节引理2：设矩阵$V_{p\times p}\ge0$，$A_{p\times q}$，那么有：

1. $\mu(A)\bigcap\mu(VA^{\bot})=\{0\}$
2. $\mu(V\vdots A)=\mu(VA^{\bot}\vdots A)$

证明：

（1）取$\forall x\in \mu(A)\bigcap\mu(VA^{\bot})$，则$\exists u,x=Au$；$\exists w,x=VA^{\bot}w$。由于$A^TA^{\bot}=0$，那么有

$$\begin{aligned} &(A^{\bot})^TA=0\\ &\Rightarrow (A^{\bot})^TAu=0\\ &\Rightarrow (A^{\bot})^Tx=0\\ &\Rightarrow (A^{\bot})^TVA^{\bot}w=0\\ &\Rightarrow w^T(A^{\bot})^TVA^{\bot}w=0 \end{aligned}$$

由于$V\ge0$，那么$\exists B,V=B^TB$，于是推出：

$$\begin{aligned} &w^T(A^{\bot})^TVA^{\bot}w=0\\ &\Rightarrow w^T(A^{\bot})^TB^TBA^{\bot}w=0\\ &\Rightarrow BA^{\bot}w=0\\ &\Rightarrow B^TBA^{\bot}w=0\\ &\Rightarrow VA^{\bot}w=0\\ &\Rightarrow x=0 \end{aligned}$$

原命题得证。

（2）不难验证，$\mu(VA^{\bot}\vdots A)\subset \mu(V\vdots A)$，下证$dim\ \mu(VA^{\bot}\vdots A)=dim\ \mu(V\vdots A)$

由于

$$\begin{aligned} dim\ \mu(VA^{\bot}\vdots A)&=rank(VA^{\bot}\vdots A)\\ &=rank(VA^{\bot})+rank(A)\\ dim\ \mu(V\vdots A)&=rank(V\vdots A) \end{aligned}$$

设$S=\{Vx|A^Tx=0,x\in\mathbb R^p\}=\mu(VA^{\bot})$，那么有：

$$\begin{aligned} rank(V\vdots A)-rank(A)&=rank\begin{pmatrix}V\\A^{\bot}\end{pmatrix}-rank(A)\\ &=dim\ S \end{aligned}$$

所以：

$$rank(V\vdots A)=rank(VA^{\bot})+rank(A)$$

于是$dim\ \mu(VA^{\bot}\vdots A)=dim\ \mu(V\vdots A)$，原命题得证。