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\) 准則。