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

%e8%a1%a8%e6%a0%bc123

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

 

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

 

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

 

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

 

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

 

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

 

我們可以將 n 個個體反應變數的值寫作一個向量   \mathbf{Y}_{n\times1}=\left[\begin{matrix} ~Y_1~\\ ~Y_2~\\  \vdots\\  Y_n \end{matrix}\right]  ,並將解釋變數寫作一個矩陣  \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]   ,那麼就可以將多元迴歸模型以矩陣的形式表達:

 

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

 

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

 

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 ,或是寫作  min_{\bold\beta}  \|\mathbf{Y}-\mathbf{X}\bold\beta\|^2_2

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

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

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

 

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

 

此時,我們可以試著在本來試圖極小化的目標函數  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,我們很少使用中文在稱呼此模型),他的正規項是係數的平方和 \sum_{j=0}^k \beta_j^2,因此ridge regression的參數估計問題變為以下形式:

 

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,其中  \lambda  為事先給定的常數。

 

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

 

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

 

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

 

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

 

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|,其中  \gamma  為事先給定的常數。

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

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

About David Huang

目前在台灣大學就讀統計碩士學位學程。我的研究領域是特徵表達與降維分析、序列決策模型、以及財務時間序列,我喜歡用商業的觀點切入大數據與資料科學!

8 Comments

  1. > 一定存在一個調整參數值 lambda_0,使得 ridge regression得出的迴規模型其平方誤差和比最小平方法得出來的迴規模型更小。

    這個寫法讓我頭痛了一陣子,頭痛的點在於直覺上 ridge 在對參數逞罰之下,配適誤差應會上升。看了一些參考資料後有些想法,但不確定是否正確。我的理解是:ridge regression得出的「模型 MSE」仍然較大,但「 MSE( hat{beta} )」較小。ridge 透過犧牲參數 unbias 取得 較小的 risk 。所以應該是

    > 一定存在一個調整參數值 lambda_0,使得 ridge regression得出的迴規模型其參數的平方誤差和比最小平方法得出來的更小。

    這樣?

    參考資料 http://stats.stackexchange.com/questions/122936/under-exactly-what-conditions-is-ridge-regression-able-to-provide-an-improvement

      1. 這個說法我覺得有一點點不夠精準,通常我們定義「最好的估計式」衡量的方法是平均平方誤差(mean squared error MSE),也就是 mathbb{E}left(|widehat{boldbeta} - bold{beta}|^2right),我們可以將此式拆解成以下形式:
        |mathbb{E}left(widehat{boldbeta} - bold{beta}right)|^2 + Var((widehat{boldbeta})=Bias^2 + Variance
        ridge regression在MSE的表現上比較好,所以我們會說他是比較好的估計式,但你的想法沒錯,ridge regression在估計時犧牲了一點點bias,但估計時的variance卻能夠大幅下降。

        至於資料略有改變對於參數的影響,這件事情跟評估迴歸分析估計式的好壞不太一樣,評估估計式的好壞是從理論觀點來看我們得出的估計式是否具有不偏性、有效性(跟估計式的variance有關)、跟一致性等良好性質,而「資料略有改變對於參數的影響」會跟最佳化或是把解釋變數展開過度時比較有關。

        一般而言,不論是ols或是ridge得出的參數在模型一樣的情況下,就算資料略有改變,估計值也不會差太多,原因是 hat{boldsymbol beta}_{ridge} = (mathbf X^T mathbf X + lambda mathbf I)^{-1} mathbf X^top mathbf y,因此只要 mathbf X^T mathbf X 以及 mathbf X不要改變太多,參數並不會差太多。

  2. 跪著看完,感謝David大神無私的分享。再請教David大神PLS regression 和 PCA regression常用嗎?若實務上有被使用,可以再與LASSO做個比較嗎?
    或是自變數和應變數都是 time series 的時候做 regression是什麼概念呢?(r arima 函數中 xreg的參數是如何估計的)

    1. 在這裡要先把PLSR拆開來討論,因為PLSR主要是用來處理多變量反應變數的情景,他可以說是由典型相關分析(canonical correlation analysis)做比較。

      PCA regression常常被使用在兩種情境,一是解釋變數 mathbf{X} 之間有高度的共線性(collinearity) 時,為了能夠讓估計出的迴歸係數其變異數不要太高,或是為了讓 mathbf{X} 不要變成不可逆矩陣,因此我們將原本的 mathbf{X} 轉換成正交主成分 (orthogonal principal component),轉換後的變數雖然沒有增加或減少,但共線性的問題卻被解決了。另外一個使用情境也是解釋變數 p 很大時,我們會希望將解釋變數進行降維度的處理 (dimension reduction),這時 PCA 算是一種降維度方法。我自己認為,PCA regression的重點是在於特徵轉換 (feature transformation),而 LASSO 的重點在於變數選擇 (variable selection),兩者在實務上嘛… 我覺得PCA regression 跟 ridge regression 算是最常看見的,LASSO則是在 p > n時非常常用,但現在更常用LASSO作者發明的 elastic net。

      關於時間序列中迴歸係數的估計,其實得先回到隨機過程 (時間序列算是離散化的隨機過程),一般來說 AR(p) 模型我們比較常透過了解 AR(p) 過程的理論性質後,推導出用來估計迴歸參數的 Euler-Walker Equation 進行估計,跟一般least square的估計方式不大一樣喔!

發表迴響

你的電子郵件位址並不會被公開。 必要欄位標記為 *