讀者提問：多元迴歸分析的變數選擇

多元迴歸分析 (multiple regression) 是統計學的重要突破，在「解釋」、「個體預測」、「趨勢預測」中都扮演著舉足輕重的地位。大鼻有位讀者詢問到「在多元迴歸模型中該如何進行變數選擇呢？」在此大鼻先簡單的說明傳統多元迴歸的模型，接著提供幾個變數選擇的方法。

多元迴歸模型與最小平方法(least square estimation)

在多元迴歸模型中，通常我們有數個解釋變數(explanatory variable)  $latex (1, X_1,X_2,\cdots,X_p)$  ，以及一個反應變數(response variable)  $latex Y$  ，在多元迴歸模型中，我們假設母體(population)的性質為：反應變數與解釋變數存在著線性關係，以條件期望值(conditional expectation)的形式表示可以寫為：

$latex \mathbb{E}(Y|X_1,\cdots,X_p) = \beta_0 \times 1+ \beta_1 X_1 + \cdots + \beta_p X_p$，其中  $latex (\beta_0,\cdots,\beta_p)$  為待估計參數。

實際上，我們會蒐集  $latex n$  個個體作為隨機樣本(random sample)，因此第 $i$ 個個體我們會得到一組解釋變數的實際值  $latex (1, X_{i1}, \cdots, X_{ip})$  與一個反應變數的實際值  $latex Y_i$  ，因此我們通常可以將樣本的迴歸模型寫為：

$latex Y_i = \beta_0 \times 1+ \beta_1 X_{i1} + \cdots + \beta_p X_{ip} + \epsilon_i$  ，其中  $latex i = 1,2,\cdots,n$。

我們可以將 $latex n$ 個個體反應變數的值寫作一個向量  $latex  \mathbf{Y}_{n\times1}=\left[\begin{matrix} ~Y_1~\\ ~Y_2~\\  \vdots\\  Y_n \end{matrix}\right]$  ，並將解釋變數寫作一個矩陣  $latex \mathbf{X}_{n\times (p+1)}=\left[\begin{matrix}
1& X_{11} & X_{12} & \cdots & X_{1p}\\
1& X_{21} & X_{22} & \cdots & X_{2p}\\
\vdots & \vdots & \vdots & \vdots & \vdots\\
1& X_{n1} & X_{n2} & \cdots & X_{np}\\
\end{matrix}\right] $  ，那麼就可以將多元迴歸模型以矩陣的形式表達：

$latex   \mathbf{Y} = \mathbf{X}\bold\beta + \bold\epsilon$  ，其中係數為  $latex \bold\beta_{(p+1)\times 1} = [\beta_0~~\beta_1 ~~ \cdots~~\beta_p]^T$   ，殘差項為 $latex  \bold\epsilon_{n\times1} = [\epsilon_1 ~~ \cdots~~\epsilon_n]^T$。

為什麼我們要寫成矩陣表達式呢？最主要的原因在於方便計算，另外我們也可以透過線性代數的理論來了解多元迴歸的投影性質，不過在此處不會探討這些問題。在我們寫出樣本的迴歸模型後，我們就需要用我們抽樣得到的解釋變數實際值  $latex \mathbf{X}$  與反應變數實際值  $latex   \mathbf{Y}$  來估計母體參數  $latex \bold\beta$  ，在傳統的多元迴歸模型，我們透過最小平方法 (least square estimation) 進行估計，目標是極小化我們估計的平方誤差和(sum of square error)，數學形式寫作：

$latex min_{\beta_0,\cdots,\beta_p} \sum_{i=1}^n\left[Y_i – \left( \beta_0 + \beta_1 X_{i1}+\cdots + \beta_pX_{ip}\right)\right]^2 $，或是寫作  $latex min_{\bold\beta}  \|\mathbf{Y}-\mathbf{X}\bold\beta\|^2_2$。

透過矩陣的微積分，我們可以很容易得出  $latex \bold\beta$  的最小平方估計式為  $latex \widehat{\bold\beta} _{ols} = (\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y}$，也就完成了整個多迴歸分析最重要的部分─參數估計。

為什麼需要進行變數選擇 (variable selection)？

大部分的時候，我們沒有辦法選出能夠完全解釋反映變數 $latex  \mathbf{Y}$  的解釋變數，有時候可能會選到不相關的變數，有時候可能會少放了重要的變數。如果我們少放了重要的變數，最小平方法得出的估計式將會是一個偏誤估計式 (biased estimator)，也就是說對於係數  $latex \beta_j$  而言，  $latex \mathbb{E}(\widehat{\beta_j}_{ols}) \neq \beta_j$  ，這樣的狀況我們稱作「省略變數的偏誤」(omitted variable bias)。

另一種情況是，我們放了太多的變數，也就是說，實際上我們只有  $latex p$  個重要的解釋變數，但我們蒐集到了  $latex k>p$   個解釋變數  $latex (X_1,\cdots,X_k)$   ，蒐集到的解釋變數中存在不重要的變數。什麼樣的情況下我們會將一個變數  $latex X_i$  稱作「不重要的變數」呢？很簡單，就是變數  $latex X_i$   對應到的迴歸係數  $latex \beta_i$  為0時，也就是解釋變數  $latex X_i$   對於反應變數  $latex Y$  不會造成任何影響。放入太多不重要的變數會有什麼問題呢？從解釋的觀點來看，當我們放入太多不重要的變數，會造成我們對於其他重要變數係數估計的品質變差，詳言之，若  $latex \beta_j$  不等於0，且迴歸模型中存在不重要的解釋變數，此時  $latex Var(\widehat{\beta_j}_{ols})$  將會變大，也就是估計的誤差變大了。從預測的角度來看，由於多元迴歸的演算法本身在極小化平方誤差和，所以變數越多在樣本內的平方誤差和通常會越低，但可能會造成過度配適 (overfitting) 的問題。變數選擇對於迴歸模型的品質有著決定性的影響，以下將會介紹三種常用的變數選擇方法，分別是「子集合選取法」、「正規化」、以及「資訊準則法」。

「解釋」的變數選擇：子集合選取法 (subset selection)

「子集合選取法」的概念很簡單，就是在所有解釋變數  $latex X_1,\cdots,X_k$  ，找出這些變數中哪一個組合有最強的「解釋能力」。舉例而言，假設我們有  $latex k$  個解釋變數，則我們將會有  $latex 2^k$  種可能的組合(每一個變數都可以選/不選)，我們的演算法是在這  $latex 2^k$  種可能的解釋變數組合中，選出一組使得調整後  $latex R^2$  最大的變數組合，作為最終的多元迴歸模型，而這樣的方法，我們叫做最佳子集合選取法 (best subset selection)。

在上述的方法中，有一項很重要的議題必須在選取最佳的子集合前先行決定，就是決定衡量「解釋能力」的指標。常見的指標有：殘差平方和 (residual sum of squares，RSS)、調整後  $latex R^2$  (adjusted   $latex R^2$  )、Mallow’s Cp 等。在決定「解釋能力」指標後，我們可以試著用窮盡法去找出最佳的子集合，但這樣的方法最大的問題在於：假設我們現在有1000個解釋變數 (常見於醫療問題)  ，則我們將會有  $latex 2^{1000}$  種可能的子集合，對電腦的運算能力看來，很難運用最佳子集合選取法進行變數選擇，因此在找出「好的子集合」的方法中，也有許多其他不同的演算法，如 forward/backward regression、stagewise regression等不同的方法。

「個體預測」的變數選擇：正規化(regularization)

另外一種被大量使用在預測方法中的變數選擇方法叫做「正規化」(regularization)，或者叫做「懲罰化」(penalization) 。這個方法的概念很簡單，假設我們現在有解釋變數   $latex (X_1,\cdots,X_k)$  ，此時最小平方法會試著進行以下最佳化：

$latex min_{\beta_0,\cdots,\beta_k} \sum_{i=1}^n\left[Y_i – \left( \beta_0 + \beta_1 X_{i1}+\cdots + \beta_kX_{ik}\right)\right]^2 $，

此時，我們可以試著在本來試圖極小化的目標函數  $latex f(\beta_0, \cdots, \beta_p) = \sum_{i=1}^n\left[Y_i – \left( \beta_0 + \beta_1 X_{i1}+\cdots + \beta_kX_{ik}\right)\right]^2$  中加上「正規項」(regularizer)或是懲罰項 (penalty term)。正規化的目的是在於，當估計出的非0參數太多時，我們認為有可能會選擇到不重要的參數，因此我們在極小會平方誤差和時，也同時要考量到非0係數的大小或個數。

常用的正規化迴歸模型有兩種，一種稱為脊迴歸 (ridge regression，我們很少使用中文在稱呼此模型)，他的正規項是係數的平方和 $latex \sum_{j=0}^k \beta_j^2$，因此ridge regression的參數估計問題變為以下形式：

$latex min_{\beta_0,\cdots,\beta_k} \sum_{i=1}^n\left[Y_i – \left( \beta_0 + \beta_1 X_{i1}+\cdots + \beta_kX_{ik}\right)\right]^2  + \lambda \sum_{j=0}^k \beta_j^2$，其中  $latex \lambda$  為事先給定的常數。

此處的  $latex \lambda$  我們稱為「調整參數」(tuning parameter)，它的目的是用來控制「懲罰」的輕重程度，  $latex \lambda$  越大，代表我們對於係數大小的懲罰越重。通常我們會透過交叉驗證 (cross validation) 的方法來決定   $latex \gamma$  的值。Ridge regression有兩個特色，一是估計出來的參數  $latex \widehat{\beta_j}_{ridge}$  將會是有偏誤的估計式，二是 ridge regression 其實跟主成分分析 (principal component analysis)有很大的關係，由於背後的理論跟線性代數中的投影(projection)有關，因此我們在此暫不討論。

另外，ridge regression 有一個非常重要的「存在性定理」(existence theorem)，這也是為什麼許多與迴歸分析有關的機器學習演算法喜歡用ridge regression的原因，「存在性定理」是針對調整參數   $latex \lambda$  的定理，其敘述如下：

一定存在一個調整參數值  $latex \lambda_0$，使得 ridge regression得出的迴規參數其平方誤差比最小平方法得出來的迴規模型更小，也就是指，  $latex \mathbb{E}\|\widehat{\bold\beta}_{ridge} – \bold{\beta}\|^2 < \mathbb{E}\|\widehat{\bold\beta}_{ols} – \bold{\beta}\|^2$  。

嚴格來說，ridge regression 並沒有真正進行「變數選擇」，主要原因在於我們可能會估計出許多靠近0但非0的  $latex \widehat{\beta_j}_{ridge}$  ，因此有了另一個有劃時代的迴歸模型誕生─LASSO (‎least absolute shrinkage and selection operator)，當時這篇論文至經已經被引用了17,067次，非常的驚人！其實LASSO的想法非常簡單，就是將正規項改為  $latex  \sum_{j=0}^k |\beta_j|$，因此LASSO其實是在解以下的最佳化問題：

$latex min_{\beta_0,\cdots,\beta_k} \sum_{i=1}^n\left[Y_i – \left( \beta_0 + \beta_1 X_{i1}+\cdots + \beta_kX_{ik}\right)\right]^2  + \gamma  \sum_{j=0}^k |\beta_j|$，其中  $latex \gamma$  為事先給定的常數。

僅僅改動了正規項，LASSO卻能夠得到一項非常了不起的結果─稀疏表達 (sparse representation)。什麼是稀疏表達呢？前面有提到，ridge regression 會估計出許多靠近0但非0的  $latex \widehat{\beta_j}_{ridge}$  ，如此一來並沒有真正的解決變數選擇的問題，而LASSO的正規化可以強制使得不重要的解釋變數  $latex X_j$  其估計出的係數  $latex \widehat{\beta_j}_{lasso} = 0$，因此也達到了變數選擇的功能。LASSO在實務上最大的缺點就是「極小化」比較困難，此外LASSO在高維度的變數選擇(也就是變數個數比樣本數還多時)，會有一些不夠好的性質。因此，還有許多其他的分枝，如：結合ridge regression 與 LASSO 的 elastic net、將變數分群化的 group LASSO等，在此就不一一介紹。

「資訊角度」的變數選擇：Akaike資訊準則

迴歸分析中還有一種很常見的變數選擇方式，是從「資訊量」的角度出發，這種類型的變數選擇方法最常被應用在時間序列預測(time series forecasting)中，假設  $latex (X_1, X_2, X_3, \cdots, X_t)$  為一時間序列，我們常常建立自我迴歸模型 (autoregressive model，AR) 來預測  $latex X_t$  ，AR(p)模型的形式如下：

$latex  X_t = \beta_0 + \beta_1  X_{t-1} + \beta_2 X_{t-2} + \cdots + \beta_p X_{t-p} + \epsilon_t$

在AR的迴歸模型中，我們最常遇到的問題是，要運用落後幾期的資訊？也就是要決定  $latex p$  的值。這時我們時常會用到Akaike資訊準則(Akaike information criterion，AIC)，AIC的背後其實隱藏著一個非常深的理論，是關於「分配間距離」的Kullback–Leibler divergence (KL divergence)。KL divergence的概念其實不難，主要是想要衡量兩個不同的分配「差距多大」，運用在此處，就是想了解「估計出的迴歸模型」與「真實的迴歸模型」的差距，實際上我們需要用觀察到的樣本來估計KL divergence，因此AIC提供了一個大樣本底下的不偏估計式：

$latex AIC = -2 \ln \left(\mathcal{L}(\beta_0,\cdots,\beta_p)\right) + 2p$，其中  $latex \mathcal{L}(\beta_0,\cdots,\beta_p)$ 是最大化的概似函數。

透過極小化AIC，我們可以選出最適合的  $latex p$  ，如此一來也就幫我們在過去無數的落後其數中選擇出重要的前  $latex p$  期，進而達到變數選擇的功能。AIC的優點在於：極小化AIC在大樣本時等同於極小化一步預測(one-step prediction)的平均平方誤差(mean squared error, MSE)，同時也針對變數進行選擇，但最大的缺點就是對於模型有分配的假設，當分配假設錯誤時AIC便不是一個衡量KL divergence的好指標。

總結：多元迴歸的變數選擇

「變數選擇」可說是迴歸分析裡面最重要的議題之一，目前常見的手法是「子集合選擇」、「正規化」、「資訊準則」等三種方法。其實，每一種方法都是一個非常非常大的問題，在這裡只能帶過最簡單的觀念，最重要的是：了解不同情境與問題下該採用哪一種變數選擇的方法！日後若有時間，會在撰寫一篇以R進行多元迴歸分析與變數選擇的示範文。

※ 大鼻每兩周會針對「讀者提問」撰寫一篇文章，目前已經被預約了支持向量機(support vector machine)囉！有任何與資料分析/統計/機器學習相關的主題，歡迎來討論！