回归分析第5讲——线性回归的诊断与优化

前面我们介绍了线性回归如何建立，如何做统计推断，以及如何完成预测。这一节我们来学习如何判断拟合的线性回归模型是否足够好，如果不好应当如何优化。为了一般性，本讲所使用的模型都是带截距项的、进行过中心化处理的线性回归模型，即：

$$Y=\alpha E_n+\widetilde X_c\beta_I+e$$

另外，为了更好的进行说明，我们补充一点记号。我们已经知道：

$$\hat Y=X\hat\beta=X(X^TX)^{-1}X^TY$$

我们记$H\equiv X(X^TX)^{-1}X^T$，这个矩阵有个专有名词，称为帽子矩阵(Hat Matrix)。

帽子矩阵是一个投影矩阵，所以它对称幂等，且有：

$$rank(H)=tr(H)=p,HX=X,HE_n=E_n$$

为了后续叙述方便，我们记帽子矩阵的对角元是$h_{ii}$。

一、线性回归模型的诊断

1.1 决定系数

判断一个线性回归模型是否能够很好的解释观测数据，仅仅依靠回归系数是否显著是不够的。某个回归系数“显著”，只能说明该自变量发生一个单位变化时，观测值会发生显著的变化，但是观测值发生的变化中有多少比例是因为自变量发生变化而导致的？更严谨一点来说，我们想知道的是观测值的方差（波动情况）有多少可以被回归模型来解释。这就导出了决定系数(Coefficient of Determination)的概念。

首先，观测值的波动情况可以由总平方和(Total Sum of Squares, TSS)来描述，即：

$$TSS=\sum_{i=1}^n(y_i-\overline y)^2$$

其次，拟合值的波动情况可以由回归平方和(Regression Sum of Squares, RSS)来描述，即：

$$RSS=\sum_{i=1}^n(\hat y_i-\overline y)^2=\hat\beta_I^T\widetilde X_c^TY$$

于是，我们可以定义决定系数为：

$$R^2=\frac{RSS}{TSS}=\frac{\hat\beta_I^T\widetilde X_c^TY}{\sum_{i=1}^n(y_i-\overline y)^2}$$

从公式中可以看到，$R^2$越大，由回归带来的波动占总波动的比例就越大，也就说明了模型对样本的解释性越好，拟合效果越好。另外，由于$R^2$是个比例形式，那么$R^2$一定在0~1之间波动，所以$R^2$越接近1，拟合效果越好。

为什么回归平方和是这样构造的？

你会发现回归平方和本质是拟合值的波动情况，也就是求每个拟合值与观测值均值的差值。那么为什么可以用拟合值和观测值均值作差？这里涉及到线性回归的一个特点：带截距项的线性回归模型一定经过样本中心点$(\overline X,\overline Y)$。这是为什么？

根据正规方程：

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

移项可得：

$$X^T(Y-X\hat\beta)=0$$

如果存在截距项，那么设计矩阵的第一列一定是$E_n$，因此如果只看第一列，那么对应的方程是：

$$E_n^T(Y-X\hat\beta)=0$$

由于$Y-X\hat\beta$代表残差，所以根据方程，残差和为0，即：

$$\sum_{i=1}^n(y_i-\hat y_i)=0$$

同时除以$n$可得$\overline y=\overline{\hat y}$。又由于：

$$\overline y=\overline{\hat y}=\frac{1}{n}X^T\hat\beta=\overline x^T\hat\beta$$

所以模型经过样本中心的$(\overline X,\overline Y)$。因此拟合值对样本均值作差计算的波动就是回归拟合的波动。

你或许会好奇TSS中不能由RSS解释的部分是什么，也就是TSS-RSS的部分。通过推导，我们可以得到：

$$TSS-RSS=\sum_{i=1}^n(y_i-\hat y_i)^2\equiv ESS$$

可以看到，TSS和RSS之间的差值就是残差平方和(Residual Sum of Squares, ESS)，它描述了观测值与拟合值之间的差距，说明观测值的波动无法完全由回归平方和解释，还存在的不可解释的差距就是残差平方和。

不同讲义对这三种平方和的记号可能会不同，有的讲义会把残差平方和写作RSS，对应“Residual”这个词，而把回归平方和重新命名为可解释的平方和(Explained Sum of Squares)，因此把回归平方和写成ESS。如果你发现其他地方有这种写法，请一定要对应好中英文的差异。

另外，本讲这样命名是因为从分子分母的构造出发，而不是从残差出发，并且笔者认为“回归平方和”这个词更加直观。不过你想用哪个记号并不重要，只要不影响定义就可以。

1.2 残差诊断

除了检验拟合程度，线性回归模型还有一个假设条件需要检验——GM假定。如果模型符合GM假定，那么模型的参数量就会下降；如果残差还服从正态分布，那么就有：

$$\hat e_i\sim N(0,\sigma^2(1-h_{ii}))$$

我们可以做标准化：

$$\frac{\hat e_i}{\sigma\sqrt{1-h_{ii}}}\sim N(0,1)$$

用样本的方差估计值代替$\sigma^2$，则：

$$r_i=\frac{\hat e_i}{\hat\sigma\sqrt{1-h_{ii}}}$$

$r_i$称为学生化内残差，是残差诊断中的重要指标。

虽然$r_i$与$\hat y_i$有轻微的相关性，但是可以近似认为$r_i\sim N(0,1)$且相互独立。因此，我们可以$\hat y_i$为横轴，$r_i$为纵轴绘制残差图，如果满足上述条件，那么图像应大致关于$r=0$对称，且分布在区间$0\pm1.96$中。

如果残差图的分布情况不理想，我们就需要分开进行检验，看看到底是哪一条不满足：

• 残差正态性诊断：检验残差是否服从正态分布，方法很多，比如Kolmogorov–Smirnov检验、Shapiro–Wilk检验等。
• 残差独立性检验：检验残差之间是否相互独立，一般用于时序数据，比如Durbin-Waston检验。
• 残差等方差检验：检验残差之间是否等方差，比如Breusch-Pagan检验(ncv检验)。

1.3 强影响点诊断

除了残差，观测值的变化也会显著影响模型的参数估计和推断，一个离群的点会极大影响拟合的结果，这些点我们统称为强影响点。下面是一张示意图，展示了强影响点对回归直线的影响：

因此，找出强影响点是回归诊断中的一个关键。找出强影响点有这么几种方法：

• 高杠杆点法
• 库克距离法
• 异常值检测

首先是高杠杆点法。所谓高杠杆点(High Leverage Point)，就是拟合直线必须经过的观测点$(x_i,y_i)$。什么意思呢？

当$h_{ii}\rightarrow1$时，不难想到$Cov(\hat e)\rightarrow0$，又由于$\hat e\rightarrow0$是模型的条件，那么此时$\hat y_i\rightarrow y_i$，也就是说拟合的回归直线几乎必须经过$(x_i,y_i)$，反过来说明这个观测点对回归结果的影响是很大的。那么什么样的点可以使得$h_{ii}$变大呢？由于$h_{ii}$：

$$\begin{aligned} h_{ii}&=\begin{pmatrix}1&x_i^T\end{pmatrix}\begin{pmatrix}X^TX\end{pmatrix}^{-1}\begin{pmatrix}1\\ x_i\end{pmatrix}\\ &=\frac{1}{n}+(x_i-\overline x)^T(\widetilde X_c^T\widetilde X_c)^{-1}(x_i-\overline x) \end{aligned}$$

所以一个点$x_i$距$\overline x$越远，$h_{ii}$就会越大，$\hat y_i$就会越接近$y_i$，那么这个点就越有将回归平面拉向自己的倾向。这样的点就是高杠杆点。

一般认为，由于$h_{ii}$平均值为$\frac{1}{n}\sum h_{ii}=\frac{1}{n}tr(H)=p/n$，所以一般认为$h_{ii}>2p/n$的点是高杠杆点。

其次是库克距离(Cook’s Distance)。库克距离的原理是，根据删除一个点前后回归系数的变化来描述一个点的影响情况。要判断第$i$个点是不是强影响点，可以去掉第$i$个点之后再拟合，得到$\hat\beta_{(i)}$，构造库克距离统计量：

$$D_i=\frac{(\hat\beta-\hat\beta_{(i)})^TX^TX(\hat\beta-\hat\beta_{(i)})}{p\hat\sigma^2}$$

$D_i$的等值线是个椭圆，如果$D_i$值越大，说明这个点对回归系数和拟合值有实质性的影响，需要郑重考量。

另外，有如下定理可方便计算库克距离：

$$D_i=\frac{1}{p}\frac{h_{ii}}{1-h_{ii}}r_i^2$$

最后是异常值检测。异常值检测与库克距离有些类似，都是要去掉一个点，然后观察模型的变化。只不过异常值检测是基于拟合值（而不是回归系数）的变化来构造统计量的。如果要判断第$i$个点的观测值$y_{(i)}$是不是异常值，可以去掉第$i$个点之后再拟合，得到$\hat y_{(i)}$，之后构造统计量，检验拟合值和真实值是否具有显著差异，即$E(y_{(i)}-\hat y_{(i)})=0$。由于$y_{(i)}$没有参与拟合，所以与$\hat y_{(i)}$独立，于是：

$$t_i=\frac{y_{(i)}-\hat y_{(i)}}{\hat\sigma_{(i)}\sqrt{1+x_i^T(X_{(i)}^TX_{(i)})^{-1}x_i}}\sim t_{n-p-1}$$

这个被称为学生化外残差，之所以是“外”残差，是因为第$i$个点没参与运算。因此，我们可以遍历每一个点，取$t_i$最大或$p$值最小时的$y_{(i)}$作为异常值。

进一步地，有如下定理：

$$t_i=\frac{\hat e_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}$$

也可以方便计算学生化外残差。

1.4 共线性诊断

共线性(collinearity)，指的是设计矩阵$X$的列向量之间存在近似的相关性，这使得$rank(X_{n\times p})<p$。为什么说共线性不是一个很好的迹象呢？因为一旦出现共线性，那么$X^TX$就包含了特征值0，矩阵是奇异的，无法求逆，回归系数就没办法求解了。

共线性的出现是无法预知的。一般来说，经过精确设计的设计矩阵是不会出现共线性的，但是在实际数据中，由于自变量的选择和取值，$X^TX$的特征值可能会出现很小的数，这就近似于共线性，对回归系数估计的伤害很大。

要考察模型的共线性，有很多种指标：

• 方差膨胀因子
• ε-病态
• 条件数
• 回归系数均方误差

方差膨胀因子(Variance Inflation Factor, VIF)是一种量化共线性程度的指标，也是实践当中非常常用的指标。VIF的原理是，对于要考察的自变量$X_j$，我们以$X_j$为"伪因变量"，剩余自变量$X_i(i=1,\cdots,p,i\ne j)$为自变量进行回归，即构建模型：

$$X_j=\alpha+X_1\beta_1+\cdots+X_{j-1}\beta_{j-1}+X_{j+1}\beta_{j+1}+\cdots+X_p\beta_p$$

该模型可以计算出一个决定系数$R_j^2$，于是我们定义VIF的计算公式为：

$$VIF_j=\frac{1}{1-R_j^2},j=1,\cdots,p$$

如果有两个自变量之间存在共线性，那么上述回归模型的决定系数变大甚至接近1，这是因为这两个自变量之间能够互相线性表示，所以拟合程度很高。此时带入VIF计算时，最终的VIF值就会变得非常大。因此，我们就可以通过VIF值来判断共线性。

一般地，我们认为VIF在0_{5之间表示不存在共线性，5}10之间表示存在弱共线性，>10表示存在强共线性。

ε-病态。给定一个很小的数$\epsilon$，如果设计矩阵$X$存在一个特征值$\lambda_i<\epsilon$，那么我们称设计矩阵$X$是ε-病态的

条件数。定义条件数$cond(X)=\sqrt{\lambda_{max}(X^TX)/\lambda_{min}(X^TX)}$，其中$\lambda()$表示求特征值。条件数如果越大，说明设计矩阵越病态，近似共线性程度越高。

上面这两个指标都是从设计矩阵的特征值来判断共线性的，最后我们还可以直接从$\hat\beta$的角度来量化，通过对$\hat\beta$的好坏进行估计，侧面反映出共线性的情况。

回归系数的均方误差(Mean Squared Error, MSE)。回归系数的均方误差的计算公式为：

$$\begin{aligned} MSE(\hat\beta)&=tr\left[E(\hat\beta-\beta)^T(\hat\beta-\beta)\right]\\ &=tr\left[Var(\hat\beta-\beta)\right]\\ &=\sigma^2tr(X^TX)^{-1}\\ &=\sigma^2\sum_{i=1}^{p-1}\frac{1}{\lambda_i} \end{aligned}$$

如果$X^TX$有一个很小的特征值，那么$MSE(\hat\beta)$就会很大，也能反映共线性问题。

二、模型补救与优化

通过第1节的一系列诊断，我们已经可以找到拟合模型中的缺陷。在这一节我们将学习如何补救、优化效果不好的模型。

2.1 Box-Cox变换

我们首先介绍的是如何处理不满足GM假定的回归模型。虽然我们在之前学习过Aitken模型和奇异线性模型的解决办法，但是模型复杂、求解复杂，在实际运用中也没有那么实用。当然，在算力强劲的信息时代也不会有这些问题，但是我们可以有更简便的处理办法。

Box和Cox在1964年提出了一套系统的方法来处理不满足GM假定的模型，通过简单的变换使得模型的正态性和方差齐性都有明显的改善。Box-Cox变换是一个含参数的变换方法，针对的是响应变量的一种变化，即：

$$Y^{(\lambda)}=\begin{cases} \frac{Y^\lambda-1}{\lambda},&\lambda\ne0\\ \ln Y,&\lambda=0 \end{cases}$$

通过转换以后，将新的$Y^{(\lambda)}$作为响应变量用于拟合，模型的正态性和方差齐性会得到好转。

至于如何估计$\lambda$，可以假设$Y^{(\lambda)}\sim N(X\beta,\sigma^2I_n)$，之后通过极大似然估计求解得到：

$$\begin{aligned} \hat\lambda&=\arg\min_\lambda (Z^{(\lambda)})^T(I_n-P_X)Z^{(\lambda)}\\ z_i^{(\lambda)}&=\begin{cases} y_1^{(\lambda)}\left(\Pi_{i=1}^ny_i\right)^{\frac{1-\lambda}{n}},&\lambda\ne0\\ (\ln y_i)\left(\Pi_{i=1}^ny_i\right)^{\frac{1}{n}},&\lambda=0 \end{cases}\\ Z^{(\lambda)}&=(z_1^{(\lambda)},\cdots,z_n^{(\lambda)})^T \end{aligned}$$

2.2 主成分回归

对于出现共线性的模型，一种常用的方法是主成分分析(Principle Component Analysis, PCA)。这里我们不阐述主成分分析的具体内容，如果你不知道也能继续学习。

在介绍主成分回归之前，我们有必要介绍线性回归的典则形式(canonical form)。设$X^TX=\Phi\Lambda\Phi^T$，其中$\Lambda=diag(\lambda_1,\cdots,\lambda_{p-1})$，$\Phi=(\phi_1,\cdots,\phi_{p-1})$，于是线性回归模型可以改写为：

$$Y=\alpha E_n+Z\gamma+e,Z=X\Phi,\gamma=\Phi^T\beta$$

假设模型满足GM假定，那么有：

$$\hat\alpha=\overline y,\hat\gamma=\Lambda^{-1}Z^TY,Var(\hat\gamma)=\sigma^2\Lambda^{-1}$$

要完成主成分回归，首先要进行主成分分析。对于典则形式，我们选取$q$个特征值$\lambda_1,\cdots,\lambda_q$，再取这些特征值对应的$\phi_i$，那么$z_i=X\phi_i$称为第$i$个主成分，于是构建新模型：

$$Y=\alpha E_n+\widetilde Z\widetilde\gamma+e$$

其中$\widetilde Z=X(\phi_1,\cdots,\phi_q)$，$\widetilde\gamma=(\gamma_1,\cdots,\gamma_q)^T$，那么最小二乘解为：

$$\hat\alpha=\overline y,\hat{\widetilde\gamma}=(\widetilde Z^T\widetilde Z)^{-1}\widetilde Z^TY$$

在原模型存在共线性时，通过保留$q$个主成分使得新模型具有更小的均方误差。当然，这种估计是有偏的。

还有一个问题是，要选择哪$q$个特征值呢？一般地，要满足

$$\frac{\sum_{i=1}^q\lambda_i}{\sum_{i=1}^{p-1}\lambda_i}\ge c$$

$c$为常数，可以取80%、85%或更大。

主成分回归也等价于下面这个优化问题的解。设$X$已经中心化，那么第一个主成分$z_1=X\phi_1$，其中：

$$\phi_1=\arg\max_{||\phi||=1}\phi^TX^TX\phi$$

第$i$个主成分$z_i=X\phi_i$，其中：

$$\phi_i=\arg\max_{||\phi||=1,\phi_l^TX^TX\phi=0}\phi^TX^TX\phi,l=1,\cdots,i-1$$

选择合适的$q\le p-1$个主成分，最后得到回归：

$$\hat{y}^{PCR}=\overline y+\sum_{i=1}^q\hat y_iz_i$$

2.3 偏最小二乘回归

2.2小节讲述的主成分回归只考虑了自变量，没有考虑响应变量。偏最小二乘法(Partial Least Squares, PLS)则利用响应变量来分解，但是基本原理是一致的。

偏最小二乘回归等价于下面这个优化问题的解。设$X,Y$均已经中心化，那么第一个PLS方向$z_1=X\phi_1$，其中：

$$\phi_1=\arg\max_{||\phi||=1}\phi^TX^TYY^TX\phi$$

第$i$个PLS方向$z_i=X\phi_i$，其中：

$$\phi_i=\arg\max_{||\phi||=1,\phi_l^TX^TX\phi=0}\phi^TX^TYY^TX\phi,l=1,\cdots,i-1$$

选择合适的$q\le p-1$个PLS方向，最后得到回归：

$$\hat{y}^{PLS}=\overline y+\sum_{i=1}^q\hat y_iz_i$$

2.4 岭估计

对于共线性问题，我们还有更现代的方法。为了能提高对病态设计矩阵的耐受性，Hoerl和Kennard在1970年提出，可以使用

$$\hat\beta(k)=(X^TX+kI_{p-1})^{-1}X^TY$$

作为$\beta$的估计值，这种估计被称为岭估计(Ridge Estimate)。

从直观上来看，当$X^TX$出现接近于0的特征值时，$X^TX+kI_{p-1}$可以使得其接近0的程度得到改善，从而提高对病态矩阵的耐受性。

岭估计又被称为收缩方法，这是因为它相当于为原始的最小二乘估计$\hat\beta$赋予了一个收缩因子：

$$\hat\beta(k)=(X^TX+kI_{p-1})^{-1}X^TY=(X^TX+kI_{p-1})^{-1}X^TX\hat\beta$$

接下来的问题是，添加$k$是否能够让回归系数的估计变好？可以证明，一定能找到这样的$k$，在实现收缩的情况下还能降低回归系数的均方误差。

证明：一定能找到这样的$k$，使得$MSE(\hat\beta(k))<MSE(\hat\beta)$。

采用典则形式来证明，那我们令

$$\begin{aligned}\hat\gamma(k)&=(Z^TZ+kI_{p-1})^{-1}Z^TY\\ &=(\Lambda+kI_{p-1})^{-1}Z^TY\end{aligned}$$

于是

$$\hat\beta(k)=(X^TX+kI_{p-1})^{-1}X^TY=\Phi\hat\gamma(k)$$

那么：

$$MSE[\hat\gamma(k)]=MSE[\hat\beta(k)]$$

对$k$求导得到：

$$\begin{aligned}\frac{\partial}{\partial k}MSE[\hat\gamma(k)]&=\frac{\partial}{\partial k}\left\{\sum_{i=1}^{p-1}\frac{k^2\gamma_i^2}{(\lambda_i+k)^2}+\sigma^2\sum_{i=1}^{p-1}\frac{\lambda_i}{(\lambda_i+k)^2}\right\}\\ &=2\sum_{i=1}^{p-1}\frac{\lambda_i(k\gamma_i^2-\sigma^2)}{(\lambda_i+k)^3}\end{aligned}$$

可以发现，当$k=0$时，$\hat\gamma(k)=\hat\gamma$，且偏导数小于0，所以一定存在$k^*$，当$k\in[0,k^*)$时，MSE单调递减，此时一定有：

$$MSE[\hat\gamma(k)]<MSE[\hat\gamma(0)]=MSE(\hat\gamma)$$

因此原命题得证。

那么如何选取这个$k$呢？Hoerl和Kennard建议取：

$$\hat k=\frac{\hat\sigma^2}{\max_{1\le i\le p-1}\hat\gamma_i^2}$$

三、基于子集的变量选择方法

前面我们介绍了如何诊断回归模型，如何补救和优化，但我们忽略了对自变量的选择问题。在实际应用中，由于我们要探讨的因素很多，所以导致自变量的数量往往巨大，所有自变量都参与运算不太现实，也没有必要，毕竟我们感兴趣的是真正显著的、有作用的变量。因此，有必要对自变量进行筛选，从而实现降维。

3.1 所有子集选取

所有子集选取应当是最直接的一种想法，也就是说，我只要把所有的自变量组合（如果自变量个数为$p$个，那么一共$2^p$种组合）都拟合一遍，再相互比较就可以了。这是没问题的，只要总的自变量个数不要太多，通过暴力拟合完全没问题。

那么现在的一个问题时，如何比较两个拟合模型呢？你当然可以使用我们之前学过的$F$检验，不过这里我们再额外提供4个更实用的准则供你使用：

• Adjusted $R^2$准则
• Mallow’s $C_p$准则
• AIC准则
• BIC准则

1. Adjusted $R^2$准则

相信你已经了解了决定系数，但是$R^2$很依赖于自变量的个数，因为不同的数据量带来的拟合效果差异很大，因此我们需要进行调整，即：

$$Adj\ R^2=1-\frac{n-1}{n-p}(1-R^2)$$

$p$是你当前模型的自变量的个数。在遍历所有的变量组合后，选取$Adj\ R^2$最大的模型即可。

2. Mallow’s $C_p$准则

假定全模型有$s$个变量，即$Y=X_{n\times s}\beta+e$，给定只有$p$个参数的子集模型，得到经验方程$\hat Y=X_p\hat\beta_p$，用$\hat Y$预测$EY$的均方误差为：

$$\begin{aligned} MSE_p&=E(\hat Y-EY)^T(\hat Y-EY)\\ &=(X_p\beta_p-EY)^T(X_p\beta_p-EY)+p\sigma^2 \end{aligned}$$

我们定义：

$$J_p=MSE_p/\sigma^2=\frac{(X_p\beta_p-EY)^T(X_p\beta_p-EY)}{\sigma^2}+p$$

虽然我们只需要找最小的$J_p$对应的回归子集就可以，但是$J_p$不是统计量，所以需要找到$J_p$的一个估计。我们可以构造：

$$C_p=\frac{(Y-\hat Y)^T(Y-\hat Y)}{\hat\sigma^2}-n+2p$$

其中$\hat\sigma^2=\frac{Y-X\hat\beta_s}{n-s}$，那么$C_p$就是一个合理的估计，虽然是有偏的但是可用。

因此，Mallow’s $C_p$准则就是找所有子集中$C_p$值最小的子集模型。

3. AIC准则

AIC准则是一种信息准则，不仅用于线性回归，在其他方面也有广泛运用。对于一组样本$Y$，含有参数向量$\theta_{p\times1}$，其极大似然函数估计为$\hat\theta$。设似然函数$l_p(Y|\theta)$，那么定义AIC值为：

$$AIC_p=-2\log l_p(Y|\hat\theta)+2p$$

应当取AIC值最小的那个模型。

例如，对于线性回归来说，可以计算得到其AIC值为：

$$AIC_p=n\log||Y-X_p^*\beta_p^*||^2+2p$$

AIC准则的一个缺点就是对$p$的惩罚不足，适合小样本。

4. BIC准则

BIC准则在AIC准则的基础上增加了对$p$的约束，通过引入与样本量$n$相关的量$C(n)$来完成。BIC中一般取$C(n)=\log n$，所以：

$$BIC_p=-2\log l_p(Y|\hat\theta)+2p\log n$$

3.2 贪婪算法/逐步回归

所有子集选取必然会碰到一个问题，就是变量太多，计算量太大。因此，可以选择一些贪婪算法来降低计算复杂度，这里我们介绍逐步回归。逐步回归(Step-wise Regression)是一种贪婪算法，其原理是，从初始模型开始不断添加/删去变量，比较前后模型的差异，直到没有显著性的变化，从而找到局部最优的子集模型。

逐步回归有三种方法：

1. 向前选择(Forward Selection)
2. 向后削去(Backward Elimination)
3. 交错逐步(Hybrid Stepwise)

向前选择，就是从只有截距项开始，依次向模型中添加一个自变量。设已经有$k$个自变量，新添加一个自变量后$\hat\beta$变成了$\widetilde\beta$，如果：

$$F=\frac{ESS(\hat\beta)-ESS(\widetilde\beta)}{ESS(\widetilde\beta)/(n-k-1)}>F_{1,n-k-1}(\alpha)$$

那么就在满足这个条件的所有$F$值中取最大的$F$值对应的自变量，加入到模型当中。如果所有新添加的自变量的$F$值都不显著，则停止。

向后削去，就是逆向的过程，从全模型开始，依次向模型中删去自变量。设已经有$k$个自变量，删除一个自变量后$\hat\beta$变成了$\widetilde\beta$，如果：

$$F=\frac{ESS(\hat\beta)-ESS(\widetilde\beta)}{ESS(\widetilde\beta)/(n-k-1)}<F_{1,n-k-1}(\alpha)$$

那么就在满足这个条件的所有$F$值中取最小的$F$值对应的自变量，从模型当中剔除。如果所有剔除的自变量的$F$值都显著，则停止。

交错逐步，也就是轮流运用向前选择和向后削去，直到停止。

逐步回归是非常常用的方法，在变量数目不太多时很好用。另外，模型比较也可以使用3.1小节的准则来进行。