回归分析09：自变量的选择(1)
Chapter 9：自变量的选择(1)
5.1 自变量选择的后果
5.1.1 全模型和选模型
第三章和第四章讨论的线性回归模型的参数估计和假设检验问题，都是基于模型设定是正确的这一前提成立，但是在处理实际问题时，我们不能确认真实的模型是什么样的，这就是模型选择问题。
这里我们只讨论线性模型，因此，模型选择问题就可以简略为自变量的选择问题。首先，我们来讨论自变量选择不当会造成什么后果。
假设初步确定一切可能对因变量 \(y\) 有影响的自变量共有 \(p\) 个，记为 \(x_1,x_2,\cdots,x_p\) ，相应的线性回归模型可以表示为矩阵形式
这里 \(X\) 为 \(n\times (p+1)\) 的列满秩设计矩阵，第一列元素全为 \(1\) ，我们将此模型称为全模型。
假设根据某些自变量选择的准则，剔除了全模型中的一些对因变量影响较小的自变量，不妨假设剔除了后面的 \(p-q\) 个自变量 \(x_{q+1},x_{q+2},\cdots,x_{p}\) ，记
可以得到一个新模型
其中 \(X_q\) 为 \(n\times(q+1)\) 的列满秩设计矩阵，\(\beta_q\) 为 \(q+1\) 维列向量，我们将此模型称为选模型。
在全模型假设下，回归系数 \(\beta\) 和 \(\sigma^2\) 的最小二乘估计为
在 \(x_0'=(x_{0q}',x_{0t}')\) 点上的预测为 \(\hat{y}_0=x_0'\hat\beta\) 。
在选模型假设下，回归系数 \(\beta_q\) 和 \(\sigma^2\) 的最小二乘估计为
在 \(x_0'=(x_{0q}',x_{0t}')\) 点上的预测为 \(\tilde{y}_{0q}=x_{0q}'\tilde\beta_q\) 。
如果对全模型的最小二乘估计 \(\hat\beta\) 作相应分块
则这里的 \(\hat\beta_q\) 和 \(\tilde\beta_q\) 并不相等，且具有完全不同的性质。
5.1.2 自变量选择对估计和预测的影响
这里我们需要引入均方误差矩阵的概念。若 \(\tilde\theta\) 是未知参数向量 \(\theta\) 的有偏估计，那么协方差阵就不能作为衡量估计精度的统计量，于是我们引入均方误差矩阵作为更合理的度量标准。
设 \(\theta\) 为一个列向量，\(\tilde\theta\) 为 \(\theta\) 的一个估计，定义 \(\tilde\theta\) 的均方误差矩阵为
注意，这里的均方误差矩阵和均方误差的定义有所区别：
容易证明，均方误差矩阵满足如下公式：
计算可得
\[\begin{aligned} {\rm MSEM}(\tilde\theta)&={\rm E}\left[(\tilde\theta-\theta)(\tilde\theta-\theta)'\right] \\ \\ &={\rm E}\left[\left(\tilde\theta-{\rm E}(\tilde\theta)+{\rm E}(\tilde\theta)-\theta\right)\left(\tilde\theta-{\rm E}(\tilde\theta)+{\rm E}(\tilde\theta)-\theta\right)'\right] \\ \\ ^*&={\rm E}\left[\left(\tilde\theta-{\rm E}(\tilde\theta)\right)\left(\tilde\theta-{\rm E}(\tilde\theta)\right)'+\left({\rm E}(\tilde\theta)-\theta\right)\left({\rm E}(\tilde\theta)-\theta\right)'\right] \\ \\ &={\rm Cov}(\tilde\theta)+\left[{\rm E}(\tilde\theta)-\theta\right]\left[{\rm E}(\tilde\theta)-\theta\right] \ . \end{aligned} \]其中带 \(^*\) 的等式中，省略了以下步骤：
\[{\rm E}\left[\left(\tilde\theta-{\rm E}(\tilde\theta)\right)\left({\rm E}(\tilde\theta)-\theta\right)'\right]={\rm E}\left[\left({\rm E}(\tilde\theta)-\theta\right)\left(\tilde\theta-{\rm E}(\tilde\theta)\right)'\right]=O \ . \]这里的 \(O\) 是全零方阵。
定理 5.1.1：假设全模型设定正确，则选模型对估计的影响有
(1) 设 \(G=\left(X_q'X_q\right)^{-1}X_q'X_t\) ，则有
所以除了 \(\beta_t=0\) 或者 \(X_q'X_t=0\) 外，\({\rm E}(\tilde\beta_q)\neq\beta_q\) ；
(2) \({\rm Cov}(\hat\beta_q)-{\rm Cov}(\tilde\beta_q)\) 为非负定矩阵；
(3) 当 \({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\) 为非负定矩阵时，\({\rm MSEM}(\hat\beta_q)-{\rm MSEM}(\tilde\beta_q)\) 为非负定矩阵；
(4) \({\rm E}\left(\tilde\sigma_q^2\right)\geq{\rm E}\left(\hat\sigma^2\right)=\sigma^2\) ，当且仅当 \(\beta_t=0\) 时等号成立。
(1) 由最小二乘估计的性质，显然 \({\rm E}(\hat\beta)=\beta\) ，下面考察 \(\tilde\beta_q\) 的均值：
\[\begin{aligned} {\rm E}(\tilde\beta_q)&=\left(X_q'X_q\right)^{-1}X_q'{\rm E}(Y) \\ \\ &=\left(X_q'X_q\right)^{-1}X_q'\begin{pmatrix} X_q & X_t \end{pmatrix}\begin{pmatrix} \beta_q \\ \beta_t \end{pmatrix} \\ \\ &=\begin{pmatrix} I_{q+1} & G \end{pmatrix}\begin{pmatrix} \beta_q \\ \beta_t \end{pmatrix} \\ \\ &=\beta_q+G\beta_t \ . \end{aligned} \]除了 \(\beta_t=0\) 或者 \(X_q'X_t=0\) 外，\({\rm E}(\tilde\beta_q)\neq\beta_q\) 由此显然成立。
(2) 回忆分块矩阵求逆公式，考虑 \(A\) 是非奇异对称矩阵的情况，
\[A=\begin{bmatrix} B & C \\ C'& D \end{bmatrix} \ , \]设 \(M=\left(D-C'B^{-1}C\right)^{-1}\) ，则有
\[A^{-1}=\begin{bmatrix} B & C \\ C' & D \end{bmatrix}^{-1}=\begin{bmatrix} B^{-1}+B^{-1}CMC'B^{-1} & -B^{-1}CM \\ -MC'B^{-1} & M \end{bmatrix} \ . \]将 \(X'X\) 用对应的分块形式表示：
\[X'X=\begin{pmatrix} X_q'\\ X_t' \end{pmatrix}\begin{pmatrix} X_q & X_t \\ \end{pmatrix}=\begin{bmatrix} X_q'X_q & X_q'X_t \\ X_t'X_q & X_t'X_t \end{bmatrix}\xlongequal{def}\begin{bmatrix} B & C \\ C' & D \end{bmatrix} \ . \]将 \(X'X\) 的逆矩阵记为
\[\left(X'X\right)^{-1}\xlongequal{def}\begin{bmatrix} B_1 & C_1 \\ C_1' & D_1 \end{bmatrix} \]因为
\[{\rm Cov}(\hat\beta)={\rm Cov}\begin{pmatrix} \hat\beta_q \\ \hat\beta_t \end{pmatrix}=\sigma^2\left(X'X\right)^{-1}=\sigma^2\begin{bmatrix} B_1 & C_1 \\ C_1' & D_1 \end{bmatrix} \ , \]所以
\[{\rm Cov}(\hat\beta_q)=\sigma^2B_1 \ . \]又因为
\[{\rm Cov}(\tilde\beta_q)=\sigma^2\left(X_q'X_q\right)^{-1}=\sigma^2B^{-1} \ . \]所以
\[{\rm Cov}(\hat\beta_q)-{\rm Cov}(\tilde\beta_q)=\sigma^2\left(B_1-B^{-1}\right)=\sigma^2B^{-1}CMC'B^{-1} \ , \]故 \({\rm Cov}(\hat\beta_q)-{\rm Cov}(\tilde\beta_q)\) 为非负定矩阵。
(3) 由均方误差矩阵的性质和前两条结论可知，
\[\begin{aligned} {\rm MSEM}(\tilde\beta_q)&={\rm Cov}(\tilde\beta_q)+\left[{\rm E}(\tilde\beta_q)-\beta_q\right]\left[{\rm E}(\tilde\beta_q)-\beta_q\right]' \\ \\ &=\sigma^2\left(X_q'X_q\right)^{-1}+G\beta_t\beta_t'G' \\ \\ &=\sigma^2B^{-1}+G\beta_t\beta_t'G' \ . \\ \\ {\rm MSEM}(\hat\beta_q)&={\rm Cov}(\hat\beta_q)=\sigma^2B_1 \ . \end{aligned} \]注意到 \(G=B^{-1}C\) ，所以当 \({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\) 为非负定矩阵时，
\[\begin{aligned} {\rm MSEM}(\hat\beta_q)-{\rm MSEM}(\tilde\beta_q)&=\sigma^2B_1-\sigma^2B^{-1}-G\beta_t\beta_t'G' \ . \\ \\ &=\sigma^2B^{-1}CMC'B^{-1}-B^{-1}C\beta_t\beta_t'C'B^{-1} \\ \\ &=B^{-1}C\left(\sigma^2M-\beta_t\beta_t'\right)C'B^{-1} \\ \\ &=B^{-1}C\left({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\right)C'B^{-1} \ . \end{aligned} \]故 \({\rm MSEM}(\hat\beta_q)-{\rm MSEM}(\tilde\beta_q)\) 为非负定矩阵。
(4) 已知 \({\rm E}(\hat\sigma^2)=\sigma^2\) ，于是
\[\begin{aligned} {\rm E}(\tilde\sigma_q^2)&=\frac{1}{n-q-1}{\rm E}\left\{Y'\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]Y\right\} \\ \\ &=\frac{1}{n-q-1}{\rm tr}\left\{\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]\cdot{\rm E}\left(YY'\right)\right\} \\ \\ &=\frac{1}{n-q-1}{\rm tr}\left\{\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]\cdot{\rm E}\left(\sigma^2I_n+X\beta\beta'X'\right)\right\} \\ \\ &=\frac{1}{n-q-1}\left\{(n-q-1)\sigma^2+\beta'X'\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]X\beta\right\} \\ \\ &=\sigma^2+\frac{1}{n-q-1}\beta'X'\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]X\beta \\ \\ &=\sigma^2+\frac{1}{n-q-1}\beta_t'X_t'\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]X_t\beta_t \\ \\ &=\sigma^2+\frac{1}{n-q-1}\beta_t'\left(X_t'X_t-X_t'X_q\left(X_q'X_q\right)^{-1}X_q'X_t\right)\beta_t \\ \\ &=\sigma^2+\frac{1}{n-q-1}\beta_t'\left(D-C'B^{-1}C\right)\beta_t \\ \\ &=\sigma^2+\frac{1}{n-q-1}\beta_t'M^{-1}\beta_t \\ \\ &\geq\sigma^2={\rm E}(\hat\sigma^2) \ . \end{aligned} \]当且仅当 \(\beta_t=0\) 时等号成立。
记全模型和选模型的预测偏差分别为
定理 5.1.2：假设全模型设定正确，则选模型对预测的影响有
(1) 设 \(G=\left(X_q'X_q\right)^{-1}X_q'X_t\) ，则有
一般情形下，\(\tilde{y}_{0q}\) 为有偏预测。
(2) \({\rm Var}(z_0)\geq{\rm Var}(z_{0q})\) ；
(3) 当 \({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\) 为非负定矩阵时，\({\rm MSE}(\hat{y}_0)-{\rm MSE}(\tilde{y}_{0q})\geq0\) 。
(1) 根据第四章点预测的性质，\({\rm E}(z_0)=0\) 显然成立。现考察 \(z_{0q}\) 的均值，由定理 5.1.1 (1) 可知
\[\begin{aligned} {\rm E}(z_{0q})&=x_0'\beta-x_{0q}'{\rm E}(\tilde\beta_q) \\ \\ &=x_0'\beta-x_{0q}'\left(\beta_q+G\beta_t\right) \\ \\ &=x_{0t}'\beta_t-x_{0q}'G\beta_t \ . \end{aligned} \](2) 根据第四章点预测的性质，容易看出
\[{\rm Var}(z_0)={\rm Var}\left(y_0-\hat{y}_0\right)=\sigma^2\left(1+x_0'\left(X'X\right)^{-1}x_0\right) \ . \]下面计算 \(z_{0q}\) 的方差，由于 \(e_0\) 和 \(e_1,e_2,\cdots,e_n\) 相互独立，所以 \(y_0\) 和 \(\tilde{y}_{0q}\) 相互独立，于是
\[\begin{aligned} {\rm Var}(z_{0q})&={\rm Var}\left(y_0-\tilde{y}_{0q}\right) \\ \\ &={\rm Var}\left(y_0\right)+{\rm Var}\left(\tilde{y}_{0q}\right) \\ \\ &=\sigma^2\left(1+x_{0q}'\left(X_q'X_q\right)^{-1}x_{0q}\right) \ . \end{aligned} \]注意到
\[\begin{aligned} &x_{0q}'\left(X_q'X_q\right)^{-1}x_{0q}=x_{0q}'B^{-1}x_{0q} \ , \\ \\ &\begin{aligned} x_{0}'\left(X'X\right)^{-1}x_{0}&=\begin{pmatrix} x_{0q}' & x_{0t}' \end{pmatrix}\begin{bmatrix} B_1 & C_1 \\ C_1' & D_1 \end{bmatrix}\begin{pmatrix} x_{0q} \\ x_{0t} \\ \end{pmatrix} \\ \\ &=x_{0q}'B_1x_{0q}+x_{0q}'C_1x_{0t}+x_{0t}'C_1'x_{0q}+x_{0t}'D_1x_{0t} \ . \end{aligned} \end{aligned} \]可以推得
\[\begin{aligned} &{\rm Var}(z_0)-{\rm Var}(z_{0q}) \\ \\ =\,&\sigma^2\left[x_0'\left(X'X\right)^{-1}x_0-x_{0q}'\left(X_q'X_q\right)^{-1}x_{0q}\right] \\ \\ =\,&\sigma^2\left[x_{0q}'\left(B_1-B^{-1}\right)x_{0q}+x_{0q}'C_1x_{0t}+x_{0t}'C_1'x_{0q}+x_{0t}'D_1x_{0t}\right] \\ \\ =\,&\sigma^2\left[x_{0q}'B^{-1}CMC'B^{-1}x_{0q}-x_{0q}'B^{-1}CMx_{0t}-x_{0t}'MC'B^{-1}x_{0q}+x_{0t}'Mx_{0t}\right] \\ \\ =\,&\sigma^2\left[x_{0q}'B^{-1}CM\left(C'B^{-1}x_{0q}-x_{0t}\right)-x_{0t}'M\left(C'B^{-1}x_{0q}-x_{0t}\right)\right] \\ \\ =\,&\sigma^2\left(x_{0q}'B^{-1}C-x_{0t}'\right)M\left(C'B^{-1}x_{0q}-x_{0t}\right) \\ \\ =\,&\sigma^2\left(C'B^{-1}x_{0q}-x_{0t}\right)'M\left(C'B^{-1}x_{0q}-x_{0t}\right)\geq0 \ . \end{aligned} \](3) 首先计算预测的均方误差
\[\begin{aligned} &{\rm MSE}\left(\hat{y}_0\right)={\rm E}\left(\hat{y}_0-y_0\right)^2={\rm E}\left(z_0^2\right)={\rm Var}(z_0) \ , \\ \\ &{\rm MSE}\left(\tilde{y}_{0q}\right)={\rm E}\left(\tilde{y}_{0q}-y_0\right)^2={\rm E}\left(z_{0q}^2\right)={\rm Var}(z_{0q})+\left[{\rm E}(z_{0q})\right]^2 \ . \end{aligned} \]由 (1) 的证明过程可得
\[\begin{aligned} \left[{\rm E}(z_{0q})\right]^2&=\left(x_{0t}'\beta_t-x_{0q}'G\beta_t\right)^2 \\ \\ &=\left(x_{0t}'-x_{0q}'G\right)\beta_t\beta_t'\left(x_{0t}'-x_{0q}'G\right)' \\ \\ &=\left(C'B^{-1}x_{0q}-x_{0t}\right)'\beta_t\beta_t'\left(C'B^{-1}x_{0q}-x_{0t}\right) \ . \end{aligned} \]所以当 \({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\) 为非负定矩阵时，根据 (2) 的证明过程可得
\[\begin{aligned} {\rm MSE}\left(\hat{y}_0\right)-{\rm MSE}\left(\tilde{y}_{0q}\right)&={\rm Var}(z_0)-{\rm Var}(z_{0q})-\left[{\rm E}(z_{0q})\right]^2 \\ \\ &=\left(C'B^{-1}x_{0q}-x_{0t}\right)'\left(\sigma^2M-\beta_t\beta_t'\right)\left(C'B^{-1}x_{0q}-x_{0t}\right) \\ \\ &=\left(C'B^{-1}x_{0q}-x_{0t}\right)'\left({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\right)\left(C'B^{-1}x_{0q}-x_{0t}\right) \ . \end{aligned} \]
定理 5.1.1 表明，即使全模型设定正确，剔除部分自变量也可以使剩余自变量的回归系数的最小二乘估计的方差减小，但此时的最小二乘估计一般是有偏的。若被剔除的自变量对因变量影响较小或难于掌握，则剔除这些自变量后可使得剩余自变量的回归系数的最小二乘估计的精度提高。
定理 5.1.2 表明，当全模型设定正确时，用选模型做预测，则预测一般是有偏的，但预测偏差的方差会有所减小。若被剔除的自变量对因变量影响较小或难于掌握，则剔除这些自变量后可使得剩余自变量的预测的精度提高。
这里，我们用非负定矩阵 \({\rm Cov}(\hat\beta_t)-\beta_t\beta_t'\) 刻画被剔除的自变量对因变量的影响，用均方误差或均方误差矩阵来刻画参数估计或预测的精度。
5.2 自变量选择的准则
5.2.0 残差平方和的局限性
首先给出以下论断：当自变量子集扩大时，残差平方和随之减少，如果按照“残差平方和越小越好”的原则选择自变量，则选入回归模型的自变量将越来越多，直到将所有自变量选入回归模型。因此，不能直接把“残差平方和越小越好”作为自变量选择的准则。
记选模型的残差平方和为 \({\rm RSS}_q\) ，则有
\[{\rm RSS}_q=Y'\left[I_n-X_q\left(X_q'X_q\right)^{-1}X_q'\right]Y \ . \]当在选模型中再增加一个自变量，不妨记为 \(x_{q+1}\) ，相应的设计矩阵为
\[X_{q+1}=\begin{pmatrix} X_q & x_{q+1} \end{pmatrix} \ , \quad x_{q+1}=\left(x_{1,q+1},x_{2,q+1},\cdots,x_{n,q+1}\right)' \ . \]相应的残差平方和为
\[{\rm RSS}_{q+1}=Y'\left[I_n-X_{q+1}\left(X_{q+1}'X_{q+1}\right)^{-1}X_{q+1}'\right]Y \ . \]由分块矩阵求逆公式可得
\[\left(X_{q+1}'X_{q+1}\right)^{-1}=\begin{bmatrix} X_q'X_q & X_q'x_{q+1} \\ x_{q+1}'X_q & x_{q+1}'x_{q+1} \end{bmatrix}^{-1}\xlongequal{def}\begin{bmatrix} \left(X_q'X_q\right)^{-1}+aba' & c \\ c' & b \end{bmatrix} \ , \]其中
\[\begin{aligned} &b=\left(x_{q+1}'x_{q+1}-x_{q+1}'X_q\left(X_q'X_q\right)^{-1}X_q'x_{q+1}\right)^{-1} \ , \\ \\ &c=-\left(X_q'X_q\right)^{-1}X_q'x_{q+1}b\xlongequal{def}-ab \ , \quad a\xlongequal{def}\left(X_q'X_q\right)^{-1}X_q'x_{q+1} \ . \end{aligned} \]于是
\[\begin{aligned} &X_{q+1}\left(X_{q+1}'X_{q+1}\right)^{-1}X_{q+1}' \\ \\ =\,&\begin{pmatrix} X_q & x_{q+1} \end{pmatrix}\begin{bmatrix} X_q'X_q & X_q'x_{q+1} \\ x_{q+1}'X_q & x_{q+1}'x_{q+1} \end{bmatrix}^{-1}\begin{pmatrix} X_q' \\ x_{q+1}' \end{pmatrix} \\ \\ =\,&\begin{pmatrix} X_q & x_{q+1} \end{pmatrix}\begin{bmatrix} \left(X_q'X_q\right)^{-1}+aba' & c \\ c' & b \end{bmatrix}\begin{pmatrix} X_q' \\ x_{q+1}' \end{pmatrix} \\ \\ =\,&X_q\left(X_q'X_q\right)^{-1}X_q'+X_qaba'X_q'+x_{q+1}c'X_q'+X_qcx_{q+1}'+x_{q+1}bx_{q+1}' \ . \end{aligned} \]注意到 \(b\) 是一非负常数，于是
\[\begin{aligned} &X_{q+1}\left(X_{q+1}'X_{q+1}\right)^{-1}X_{q+1}'-X_q\left(X_q'X_q\right)^{-1}X_q' \\ \\ =\,&X_qaba'X_q'+x_{q+1}c'X_q'+X_qcx_{q+1}'+x_{q+1}bx_{q+1}' \\ \\ =\,&b\left[X_qaa'X_q'-x_{q+1}a'X_q'-X_qax_{q+1}'+x_{q+1}x_{q+1}'\right] \\ \\ =\,&b\left[X_qa\left(a'X_q'-x_{q+1}'\right)-x_{q+1}\left(a'X_q'-x_{q+1}'\right)\right] \\ \\ =\,&b\left(X_qa-x_{q+1}\right)\left(X_qa-x_{q+1}\right)' \ . \end{aligned} \]即 \(X_{q+1}\left(X_{q+1}'X_{q+1}\right)^{-1}X_{q+1}'-X_q\left(X_q'X_q\right)^{-1}X_q'\) 是非负定矩阵，所以
\[{\rm RSS}_{q+1}-{\rm RSS}_q=Y'\left[X_q\left(X_q'X_q\right)^{-1}X_q'-X_{q+1}\left(X_{q+1}'X_{q+1}\right)^{-1}X_{q+1}'\right]Y\leq 0 \ . \]所以增加自变量可以使得残差平方和减少。
5.2.1 平均残差平方和准则
由于 \({\rm RSS}_q\) 随着 \(q\) 的增大而下降，为了防止选取过多的自变量，一个常见的做法是对 \({\rm RSS}_q\) 乘上一个随着 \(q\) 的增大而增大的函数，作为惩罚因子。定义平均残差平方和
我们按照 \({\rm RMS}_q\) 越小越好的准则选择自变量，并称其为平均残差平方和准则或 \({\rm RMS}_q\) 准则。
5.2.2 调整后的 \(R^2\) 准则
判定系数 \(R^2\) 度量了数据与模型的拟合程度，我们自然希望它越大越好。但根据定义
我们不能直接将 \(R^2_q\) 作为选择自变量的准则，否则会和 \({\rm RSS}_q\) 一样将所有的自变量选入模型。为了克服以上缺点，我们引入调整后的判定系数：
容易证明 \(\bar{R}_q^2\leq R_q^2\) ，且 \(\bar{R}_q^2\) 不一定随着自变量个数 \(q\) 的增大而增大。因此，我们按照 \(\bar{R}_q^2\) 越大越好的准则选择自变量，并称其为调整后的 \(R^2\) 准则。
