LASSO與其衍生方法之特性比較 - 政大學術集成

全文

(1)國立政治大學統計學系研究所碩士學位論文. LASSO 與其衍生方法之特性比較政治. 大. 立 of LASSO and Its Derivative Methods Property Comparison ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 指導教授：蔡政安博士薛慧敏博士研究生：黃昭勳撰. 中華民國一 ○ 六年六月.

(2) 摘要本論文比較了幾種估計線性模型係數的方法，包括 LASSO、Elastic Net、 LAD-LASSO、EBLASSO 和 EBENet。有別於普通最小平方法，這些方法在估計模型係數的同時，能夠達到變數篩選，也就是刪除不重要的解釋變數，只將重要的變數保留在模型中。在現今大數據的時代，資料量有著愈來愈龐大的趨勢，其中不乏上百個甚至上千個解釋變數的資料，對於這樣的資料，變數篩選就顯得更加重要。本文主要目的為評估各種估計模型係數方法的特性與優劣，當中包含了兩種模擬研究與兩筆實際資料應用。由模擬的分析結果來看，每種估計方法都有. 政治大. 不同的特性，沒有一種方法使用在所有資料都是最好的。. 立. ‧. ‧ 國. 學. 關鍵詞：Elastic Net、LASSO、懲罰函數、迴歸、變數篩選. n. er. io. sit. y. Nat. al. Ch. engchi. I. i n U. v.

(3) Abstract In this study, we compare several methods for estimating coefficients of linear models, including LASSO, Elastic Net, LAD-LASSO, EBLASSO and EBENet. These methods are different from Ordinary Least Square (OLS) because they allow estimation of coefficients and variable selection simultaneously. In other words, these methods eliminate non-important predictors and only important predictors remain in the model. In the age of big data, quantity of data has become larger and larger. A datum with hundreds of or thousands of predictors is also common. For this type of. 政治大 to compare properties of different 立 variable selection methods as well as to find which. data, variable selection is apparently more essential. The primary goal of this article is. ‧ 國. 學. method best fits a large number of data. Two simulation scenarios and two real data applications are included in this study. By analyzing results from the simulation study,. ‧. we can find that every method enjoys different characteristics, and no standard. al. er. io. sit. y. Nat. method can handle all kinds of data.. n. Keywords: Elastic Net; LASSO; penalty function; regression; variable selection. Ch. engchi. II. i n U. v.

(4) 目次第一章研究背景......................................................................................................1 第二章研究方法......................................................................................................4 第三章模擬研究......................................................................................................9 第一節前言......................................................................................................9 第二節模擬過程..............................................................................................9 第三節模擬結果與討論..................................................................................13 第四章變數分群模擬研究......................................................................................19. 政治大第二節模擬過程..............................................................................................19 立第一節前言......................................................................................................19. ‧ 國. 學. 第三節模擬結果與討論..................................................................................20 第五章實際資料應用..............................................................................................25. ‧. 第一節攝護腺癌 (Prostate Cancer) 研究應用 ..............................................25. sit. y. Nat. 第二節白血病 (Leukemia) 研究應用 ...........................................................28. al. er. io. 第六章結論..............................................................................................................39. v. n. 參考文獻....................................................................................................................40. Ch. engchi. III. i n U.

(5) 表次表 3-1. 估計係數評估指標 ...................................................................................12. 表 3-2. 正確配適、過度配適、低度配適定義 ...................................................13. 表 3-3. εi ～N(0, 1)模擬結果 .................................................................................14. 表 3-4. εi ～dbexp(0, 1/√2)模擬結果 ....................................................................15. 表 3-5. εi ～Cauchy(0, 1)模擬結果 .......................................................................16. 表 4-1. x 的相關係數矩陣 ....................................................................................20. 表 5-1. 攝護腺癌資料：各方法估計係數 ...........................................................26. 表 5-2. 攝護腺癌資料：各方法平均預測均方差，括弧內為標準差 ...............27. 表 5-3. 政治大白血病資料：各方法平均基因篩選個數 ...............................................29 立白血病資料：各方法基因篩選與次數 ...................................................30. 表 5-5. 白血病資料：各方法平均預測正確率 ...................................................31. 表 5-6. 基因選擇一：LASSO (篩選次數至少 25 次) .........................................32. 表 5-7. 白血病資料：各主成分資訊 (I) .............................................................33. 表 5-8. 基因選擇二：Elastic Net (篩選次數至少 25 次) ....................................34. 表 5-9. 白血病資料：各主成分資訊 (II) ............................................................35. 表 5-10. 基因選擇三：Elastic Net (篩選次數 50 次) ............................................36. 表 5-11. 白血病資料：各主成分資訊 (III) ...........................................................37. 表 5-12. 白血病資料：各基因選擇個數與模型平均預測正確率 .......................38. ‧. ‧ 國. 學. 表 5-4. n. er. io. sit. y. Nat. al. Ch. engchi. IV. i n U. v.

(6) 圖次圖 1-1. p=2 時 LASSO (左) 和脊迴歸 (右) 係數估計示意圖 ..........................2. 圖 3-1. N(0, 1), dbexp(0, 1/√2), Cauchy(0, 1)機率分佈圖 ..........................10. 圖 4-1. LASSO 係數路徑圖..................................................................................20. 圖 4-2. Elastic Net 係數路徑圖.............................................................................21. 圖 4-3. LAD-LASSO 係數路徑圖 ........................................................................21. 圖 4-4. EBLASSO 係數路徑圖 ............................................................................22. 圖 4-5. EBENet 係數路徑圖 .................................................................................22. 圖 4-6. L1-Norm 與λ關係圖 (LASSO) ................................................................23. 圖 4-7. 政治大 LASSO 係數路徑圖立 (橫軸為λ) ...............................................................23 白血病資料：雙標圖 (I) .........................................................................33. 圖 5-2. 白血病資料：雙標圖 (II) ........................................................................35. 圖 5-3. 白血病資料：雙標圖 (III).......................................................................37. ‧. ‧ 國. 學. 圖 5-1. n. er. io. sit. y. Nat. al. Ch. engchi. V. i n U. v.

(7) 第一章研究背景考慮一筆 n 個樣本、p 個解釋變數的資料：(xi1 , xi2 , … , xip , yi ), i = 1, 2, … , n。其中xi1 , xi2 , … , xip 為解釋變數，yi 為反應變數，目標為建立一線性模型 (線性迴歸、羅吉斯迴歸 (logistic regression))。以線性迴歸模型 Yi = β0 + β1 Xi1 + β2 Xi2 + ̂ = (β̂0 , β̂1 , … , β̂p )，普通最小平方估計 ⋯ + βp Xip + εi 為例，欲得到估計係數 𝛃 量 (Ordinary Least Square (OLS) Estimator) 為一種常見的估計迴歸模型係數方法，其做法是最小化殘差平方和，求得估計係數，即 ̂ = arg minβ ∑ni=1(yi − β0 − ∑p βj xij )2。 𝛃 j=1. 政治大. (1). 然而，OLS 建立的模型，經常無法令資料分析者感到滿意，當中可以分成兩個. 立. 部分來討論：第一是預測準確度，OLS 估計量偏誤 (bias) 小但變異 (variance) 大，. ‧ 國. 學. 導致預測準確度差。第二是解釋力，OLS 無法篩選變數，若模型中解釋變數多，分析者難以判斷哪些是真正顯著的變數，也就不容易有很好的解釋力 (Tibshirani,. ‧. 1996)。. y. Nat. sit. 為了改善 OLS 估計係數變異大的問題，過去曾提出脊迴歸 (ridge regression). al. n. p. er. io. (Hoerl and Kennard, 1988)，其延用了 OLS 最小化殘差平方和的想法，並將係數. i n U. v. 的平方總和 (∑j=1 β2j , 即 L2-Norm) 限制在某常數以內，即. Ch. engchi. ̂ = arg minβ ∑ni=1(yi − β0 − ∑p βj xij )2 s.t. ∑p β2j ≤ t。 𝛃 j=1 j=1. (2). p. 加入限制式∑j=1 β2j ≤ t後，估計係數的偏誤會增加，但能夠縮小係數，降低其變異。換句話說，藉由此限制式，犧牲一些偏誤，以得到變異較小的估計係數，增加預測準確度。然而，脊迴歸無法使任何估計係數為 0，因此無法篩選變數，增加模型解釋力。有鑑於脊迴歸無法達到變數篩選，Tibshirani (1996) 提出了 Least Absolute Shrinkage and Selection Operator (LASSO)，其延續了脊迴歸的想法，只是將限制 p. 式改為係數絕對值總和 (∑j=1 |βj |, 即 L1-Norm) 小於或等於某常數，即 1.

(8) ̂ = arg minβ ∑ni=1(yi − β0 − ∑p βj xij )2 s.t. ∑p |βj | ≤ t。 𝛃 j=1 j=1. (3). p. 有別於脊迴歸，LASSO 的限制式∑j=1 |βj | ≤ t不但能縮小估計係數，減少其變異，增加預測準確度，還能使部分估計係數為 0，亦即刪除不顯著的變數，達到變數篩選，增加模型解釋力。其中 t 的大小可以控制係數縮減程度，令β̂oj 為 OLS 估 p 計係數，則當 t 滿足0 < t < ∑j=1 |β̂oj |有係數縮減效果，且 t 愈小則估計係數縮減 p 程度愈大，篩選的變數也就愈少。若t ≥ ∑j=1 |β̂oj |，則 LASSO 估計係數解與 OLS. 相同。. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. 圖 1-1. Ch. engchi. i n U. v. p=2 時 LASSO (左) 和脊迴歸 (右) 係數估計示意圖 (Tibshirani, 1996). 圖 1-1 為 LASSO (左) 和脊迴歸 (右) 在 p=2 時的係數估計示意圖，其中橫軸和縱軸分別為β1 和β2 。黑色區域為限制式範圍，分別是 LASSO (左) 的 ̂ 𝐨 為 OLS 估計係數解，此點是殘差 |β1 | + |β2 | ≤ t和脊迴歸 (右) 的β12 + β22 ≤ t。𝛃 平方和最小的點，其周圍橢圓形為殘差平方和的等高線，由內往外逐漸擴大。 LASSO 或脊迴歸在求解估計係數時，欲最小化殘差平方和，且同時必須滿足限 ̂ 𝐨 有最小的殘差平方和，但未落在限制式範圍內。為了同制條件。從圖形來看，𝛃 2.

(9) 時達成最小化殘差平方和與滿足限制條件，等高線由內往外第一個接觸到黑色區域的點，即為 LASSO (左) 和脊迴歸 (右) 的估計係數解。由於 LASSO 的限制式範圍有「角」，等高線容易交會在角上，因而形成 0 係數。反觀脊迴歸，限制式範圍是平滑的，沒有角能夠交會，因此不容易形成 0 係數。運用圖 1-1 二維的幾何概念，可以推廣至多維度的解釋變數，當變數多時，LASSO 可以達到變數篩選，脊迴歸則無法。 p. p. LASSO 的限制式∑j=1 |βj | ≤ t可以改成λ ∑j=1 |βj |的形式，並與殘差平方和相加，同時最小化求解 (Gill et al., 1981)，亦即原先的式子 (3) 可以改寫為. 政治大其中包含了兩個部分：損失函數立 (loss function) ∑ (y − β. ̂ = arg minβ {∑ni=1(yi − β0 − ∑p βj xij )2 + λ ∑p |βj |}, λ > 0。 𝛃 j=1 j=1 n i=1. i. (4). p. 0. − ∑j=1 βj xij )2 和懲罰. p. ‧ 國. 學. 函數 (penalty function) λ ∑j=1 |βj |。關於這兩種函數，第二章將有更完整的說明。 LASSO 由 Tibshirani (1996) 所提出，日後統計學家又以此方法為基礎，發. ‧. 展出幾種衍生的方法。第二章針對 LASSO 與幾種其衍生的估計線性模型係數方. sit. y. Nat. 法做介紹。第三章模擬當反應變數服從常態分配 (normal distribution) 或存在厚. al. er. io. 尾偏差 (heavy-tailed error) 時，對各方法的影響，並考慮在不同的樣本數 (n) 與. v. n. 解釋變數個數 (p) 情況下做比較。第四章做變數分群模擬，並根據係數路徑圖. Ch. engchi. i n U. (coefficient path) 來觀察各方法是否存在變數分群效應 (grouping effect)。第五章將兩筆實際資料「攝護腺癌 (Prostate Cancer) 」臨床研究資料及「白血病 (Leukemia)」微陣列資料 (microarray data) 套入各方法估計模型係數，並比較各方法的好壞。第六章為總結。. 3.

(10) 第二章研究方法以下將介紹 LASSO 與幾種其衍生的估計線性模型係數方法。. LASSO (Tibshirani, 1996) 此方法的係數估計式子為 ̂ = arg minβ {∑ni=1(yi − β0 − ∑p βj xij )2 + λ ∑p |βj |}, λ > 0。 𝛃 j=1 j=1 p. p. 其中包含兩個部分：損失函數∑ni=1(yi − β0 − ∑j=1 βj xij )2 和懲罰函數λ ∑j=1 |βj |。損失函數的部分延用了 OLS 的概念，即最小化殘差平方和，盡可能使誤差達到. 政治大 p. p. 最小。懲罰函數則如同第一章提到的，λ ∑j=1 |βj |等價於限制式∑j=1 |βj | ≤ t，能. 立. 縮小估計係數，減少其變異，且根據圖 1-1，這樣的限制條件能使部分估計係數. ‧ 國. 學. 為 0，亦即刪除不顯著的變數，達到變數篩選。當中的λ為調諧參數 (tuning parameter)，其值大小會影響係數估計結果，λ愈大則係數絕對值會壓縮得愈小，. ‧. 篩選的變數也就愈少。此方法的運算流程如下：. y. Nat. sit. 1. 給定λ (假設為λm )。其最佳值可經由交叉驗證 (cross validation) 來決定。. n. al. er. io. 2. 最小化係數估計式子，即. i n U. v. ̂ = arg minβ {∑ni=1(yi − β0 − ∑p βj xij )2 + λm ∑p |βj |}。 𝛃 j=1 j=1. Ch. engchi. ̂ = (β̂0 , β̂1 , … , β̂p )。若β̂j ≠ 0，則第 j 個解釋變數被篩選，反 3. 得到估計係數 𝛃 之則被刪除。. Elastic Net (Zou and Hastie, 2005) 此方法的係數估計式子為 ̂ = arg minβ {∑ni=1(yi − β0 − ∑p βj xij )2 + λ (α ∑p |βj | + (1 − α) ∑p β2j )} 𝛃 j=1 j=1 j=1 , λ > 0, 0 < α < 1。. 4.

(11) p 其中包含兩個部分：損失函數 ∑ni=1(yi − β0 − ∑j=1 βj xij )2 和懲罰函數 p. p. λ(α ∑j=1 |βj | + (1 − α) ∑j=1 β2j )。雖然 LASSO 改善了脊迴歸的缺陷，可以達成變數篩選，但當資料的解釋變數多於樣本 (p > n) 時，LASSO 最多只能選出 n 個變數 (Efron et al., 2004)。為了解決此使用限制，Zou 和 Hastie (2005) 提出了 Elastic Net，其損失函數的部分同樣延用了 OLS 的概念。懲罰函數則包含係數絕 p. 對值總和 (∑j=1 |βj |, 即 L1-Norm) 與係數平方總和 (∑pj=1 β2j , 即 L2-Norm)，相當於 LASSO 與脊迴歸懲罰函數的結合，α與1 − α分別為兩者的權重。相較於 LASSO， Elastic Net 能選出較多的變數，且當p > n時，Elastic Net 能篩選多於 n 個變數。. 政治大中的一個，而 Elastic Net 則有篩選一整群的趨勢，且能夠達到變數分群。在本文立當資料中存在幾個高度相關且顯著的解釋變數時，LASSO 傾向於只選出一群當. 中，α皆設定為 0.5，也就是給予 LASSO 與脊迴歸的懲罰函數相等的權重。λ為. 1. 給定λ (假設為λm )。其最佳值可經由交叉驗證來決定。. sit. y. Nat. 2. 最小化係數估計式子，即. ‧. ‧ 國. 學. 調諧參數。此方法的運算流程如下：. io. al. er. ̂ = arg minβ {∑ni=1(yi − β0 − ∑p βj xij )2 + λm (0.5 ∑p |βj | + 0.5 ∑p β2j )}。 𝛃 j=1 j=1 j=1. v. n. ̂ = (β̂0 , β̂1 , … , β̂p )。若β̂j ≠ 0，則第 j 個解釋變數被篩選，反 3. 得到估計係數 𝛃 之則被刪除。. Ch. engchi. i n U. LAD-LASSO (Gao and Huang, 2010) 此方法的係數估計式子為 ̂ = arg minβ {∑ni=1 |yi − β0 − ∑p βj xij | + λ ∑p |βj |}, λ > 0。 𝛃 j=1 j=1 p. p. 其中包含兩個部分：損失函數∑ni=1 |yi − β0 − ∑j=1 βj xij | 和懲罰函數λ ∑j=1 |βj | 。 LAD-LASSO 延續了 LASSO 的做法，只是將損失函數的部分由 LASSO 的最小殘差平方形式改為最小絕對離差 (Least Absolute Deviation, LAD)。懲罰函數則同樣為 LASSO 的 L1-Norm 型式。當反應變數存在離群值 (outliers) 或厚尾偏差時， 5.

(12) LAD-LASSO 的係數估計結果較為穩健 (robust)。λ為調諧參數。此方法的運算流程如下： 1.. 給定λ (假設為λm )。其最佳值可經由交叉驗證來決定。. 2.. 最小化係數估計式子，即 ̂ = arg minβ {∑ni=1 |yi − β0 − ∑p βj xij | + λm ∑p |βj |}。 𝛃 j=1 j=1. 3.. ̂ = (β̂0 , β̂1 , … , β̂p )。若β̂j ≠ 0，則第 j 個解釋變數被篩選，反得到估計係數 𝛃 之則被刪除。. 政治大此方法全名為 Empirical 立 Bayesian LASSO，其係數需滿足階層事前分配. EBLASSO (Cai et al., 2011). ‧ 國. 學. (hierarchical prior distribution). βj ～N(0, σ2j )，σ2j ～exp(λ)，λ～gamma(a, b), j = 1, 2, … , p。. ‧. 前面介紹的三種方法皆是利用懲罰函數來壓縮估計係數和篩選變數，現在將 p. sit. y. Nat. LASSO 懲罰函數λ ∑j=1 |βj | (L1-Norm) 的想法延伸到貝氏，此函數對估計係數的. v ni. n. al. er. io. 懲罰機制可以視為係數βj 服從階層事前分配βj ～N(0, σ2j )，σ2j ～exp(λ)。給定λ後， λ. 可以得到 βj 的邊際分配 (marginal distribution) 為 p(βj ) = √2 exp(−√2λ|βj |)，當. Ch. U i e h n c g 中的√2λ|βj |即為 L -Norm 懲罰函數的型式。同樣地，λ的大小會影響估計係數的 1. 壓縮程度，為了避免選擇不適當的λ，將λ再加入一個階層的分配：λ～gamma(a, b)，其中 a 和 b 為超參數 (hyperparameters)。EBLASSO 以𝐮 ̂ (𝛃′ 事後平均數 (posterior mean)) 作為估計係數，其中𝛃′ 為所有係數中非 0 的部分。此方法的運算流程如下： 1.. 給定 a、b (假設為am 、bm )。其最佳值可經由交叉驗證來決定。. 2.. 目標模型：𝐘 = μ + 𝐗𝛃 + 𝐞，其中𝐞～N(0, σ20 𝐈)。. 3.. 設αj = 1/σ2j , j = 1, 2, … , p，並令 𝛂 = (α1 , α2 , … , αp )。. 4.. μ, σ0 , 𝛂為需要被估計的參數。令 𝛉 = (μ, σ0 , 𝛂)。 6.

(13) 5.. 𝛉的 log 邊際事後分配 (marginal posterior distribution)： 1. p. L(𝛉) = − 2 [log|𝐂| + (𝐘 − μ)T 𝐂 −1 (𝐘 − μ)] − ∑j=1(am + 1)log. 1+bm αj bm αj. + constant. p. ，其中 𝐂 = σ20 𝐈 + ∑j=1 αj−1 𝐱 j 𝐱 jT。 6.. ̂ = arg max𝛉 L(𝛉) = (μ̂, σ ̂) ，其中最大化 L(𝛉) ，得到 𝛉 估計值，即 𝛉 ̂0 , 𝛂 ̂ = (α 𝛂 ̂1 , α ̂2 , … , α ̂p )。. 7.. 若α ̂j = ∞，由αj = 1/σ2j ，可得σ ̂2j = 0，再由βj ～N(0, σ2j )，可得β̂j = 0。簡言之，若α ̂j = ∞，則β̂j = 0 (第 j 個解釋變數被刪除)，藉此達到變數篩選。. 8.. 將模型簡化，考慮所有 0 係數被排除後的模型：𝐘 = μ + 𝐗 ′ 𝛃′ + 𝐞，其中𝛃′ 為. 政治大 𝛃 事後分配 (posterior 立distribution) 服從多元常態分配 (multivariate normal 所有係數中非 0 的部分，𝐗 ′ 為這些係數對應的變數。. 9.. ′. ‧ 國. T. ′ T ′ −1 ̂ = (diag(𝛂′ ) + σ−2 ∑ ，期望值 0 𝐗 𝐗 ). 學. distribution) ，共變異數. 𝐮 ̂=. ′ ̂ ′ ̂ 中非無限大的元素，μ = 𝟏𝐓 𝐘/n。 σ−2 0 ∑𝐗 (𝐘 − μ)，其中 𝛂 為𝛂. ‧. 10. 以𝐮 ̂ (𝛃′ 事後平均數) 作為𝛃′ 的估計係數。. sit. y. Nat. al. er. io. EBENet (Huang et al., 2015). v. n. 此方法全名為 Empirical Bayesian Elastic Net，其係數需滿足階層事前分配. Ch. engchi. i n U. −1 βj ～N(0, (λ1 + σ ̃−2 ̃2j ～generalized gamma(λ1 , λ2 ) j ) )，σ. , j = 1, 2, … , p。當λ1 > 0, λ2 > 0，且 c = √λ1 λ2 /π exp(−λ2 /λ1 ) 時，f(σ ̃2j )可以改寫成平移的伽 1. 瑪分配 (gamma distribution). f(σ ̃2j ). =. λ2 2. ̃2j 1 (σ Γ( ) 2. 1. +. 1 −2 ) exp(−λ2 (σ ̃2j λ1. 1. + λ ))，其中λ1 1. 和λ2 為超參數。給定λ1 和λ2 後，可以得到 p(βj ) ∝ exp(−(√2λ2 |βj | + 中的√2λ2 |βj | +. λ1 2. λ1 2. β2j ))，當. β2j 即為 Elastic Net 懲罰函數的型式 (L1-Norm+L2-Norm)。當. ̃2j ～exp(λ2 )，相當於 EBLASSO 的事前分配 λ1 = 0, λ2 > 0時，βj ～N(0, ̃ σ2j )，σ (prior distribution)，為此方法的一個特例。EBENet 以𝐮 ̂ (𝛃′ 事後平均數) 作為估 7.

(14) 計係數，其中𝛃′ 為所有係數中非 0 的部分。如同 Elastic Net 相較於 LASSO，EBENet 相較於 EBLASSO 能選出較多的變數。當資料中存在幾個高度相關且顯著的解釋變數時，EBENet 有篩選一整群的趨勢，且能夠達到分群。此方法的運算流程如下： 1.. 給定λ1 、λ2 (假設為λ1m 、λ2m )。其最佳值可經由交叉驗證來決定。. 2.. 目標模型：𝐘 = μ + 𝐗𝛃 + 𝐞，其中𝐞～N(0, σ20 𝐈)。. 3.. ̃ = (α 設α ̃j = 1/σ ̃2j , j = 1, 2, … , p，並令 𝛂 ̃1 , α ̃2 , … , α ̃p )。. 4.. ̃ 為需要被估計的參數。令 𝛉 = (μ, σ0 , 𝛂 ̃ )。 μ, σ0 , 𝛂. 5.. 𝛉的 log 邊際事後分配： 1. 立. 政治大 1. λ1m. 2. ̃j α. p L(𝛉) = − [log|𝐂| + (𝐘 − μ)T 𝐂 −1 (𝐘 − μ)] − ∑j=1 log ( 2. p. + 1) − ∑j=1. λ2m ̃j α. + constant. p. ‧ 國. 6.. 學. ̃j )−1 𝐱 j 𝐱 jT。，其中 𝐂 = σ20 𝐈 + ∑j=1(λ1m + α. ̂) ，其中 ̂ = arg max𝛉 L(𝛉) = (μ̂, σ ̃ 最大化 L(𝛉) ，得到 𝛉 估計值，即 𝛉 ̂0 , 𝛂. ‧. ̂ = (α ̂1 , α ̂2 , … , α ̂p )。 ̃ 𝛂 ̃ ̃ ̃. −1 ̂ ̂ 若α ̃j = ∞，由α ̃2j ，可得σ ̃−2 ̃−2 ̃j = 1/σ j = ∞，再由βj ～N(0, (λ1 + σ j ) )，可得. sit. y. Nat. 7.. 8.. al. n. 變數篩選。. er. io. ̂j = ∞，則β̂j = 0 (第 j 個解釋變數被刪除)，藉此達到 β̂j = 0。簡言之，若α ̃. Ch. engchi. i n U. v. 將模型簡化，考慮所有 0 係數被排除後的模型：𝐘 = μ + 𝐗 ′ 𝛃′ + 𝐞，其中𝛃′ 僅包含非 0 係數，𝐗 ′ 為這些係數對應的變數。. 9.. ′ T ′ −1 ̂ = (diag(𝛂′ ) + σ−2 𝛃′ 事後分配服從多元常態分配，共變異數 ∑ 0 𝐗 𝐗 ) ，期 T. ̂ 𝐗 ′ (𝐘 − μ)，其中 𝛂′ 為𝛂 ̂ 中非無限大的元素，μ = 𝟏𝐓 𝐘/n。望值 𝐮 ̂ = σ0−2 ∑ 10. 以𝐮 ̂ (𝛃′ 事後平均數) 作為𝛃′ 的估計係數。. [註]：σ ̃2j ～generalized gamma(λ1 , λ2 )，σ ̃2j 之機率密度函數為 f(σ ̃2j ) = c(λ1 σ ̃2j + 1)−1/2 exp(−λ2 σ ̃2j ), j = 1, 2, … , p (c 為一常數)。. 8.

(15) 第三章模擬研究第一節前言在第二章介紹了幾種估計線性模型係數方法，究竟這些方法的估計表現如何？當反應變數服從常態分配或存在厚尾偏差時，各方法的係數估計表現是好是壞？另外，在不同的樣本數 (n) 與解釋變數個數 (p) 之下，各方法的表現又是如何？為了研究這些議題，設計了以下的模擬情境，並藉由估計係數評估指標來評估各方法的估計表現，其中可以分成「變數篩選表現」和「估計係數準確度」來討論。. 治政第二節模擬過程大立學. ‧ 國. 模擬假設Xij 與Yi 滿足關係式 p. Yi = ∑j=1 Xij βj + εi , i = 1, 2, … , n。. io. sit. y. s = 1, 2, … , p。. 2, j = 1, 2, … , 5. er. 真實係數. Nat. Corr(Xr , Xs ) = 0.5|r−s| , r = 1, 2, … , p,. ‧. 其中Xij 利用 N(0, 1)隨機產生，且兩兩變數間之相關係數滿足. n. ， a l βj = {0, o. w. v i n C h ，希望藉由模擬來看各方法能否將即 p 個解釋變數中有 5 個為顯著變數 5 個顯著 engchi U 變數估計為非 0，並將 p-5 個不顯著的估計為 0。為了觀察當反應變數服從常態分配或存在厚尾偏差時，各方法的係數估計表現，考慮三種隨機誤差項 (εi ) 的分配 (Gao and Huang, 2010)： I.. εi ～N(0, 1)：反應變數無厚尾偏差。. II.. εi ～dbexp(0, 1/√2)：反應變數存在輕微厚尾偏差。. III. εi ～Cauchy(0, 1)：反應變數存在嚴重厚尾偏差。以 N(0, 1)為基準，若分配的機率分佈兩側比 N(0, 1)厚，則存在厚尾現象。圖 3-1 為此三種分配之機率分佈圖，可以發現 dbexp(0, 1/√2)的兩側略比 N(0, 1)厚，存 9.

(16) 在輕微厚尾現象；Cauchy(0, 1)的兩側比 N(0, 1)厚許多，存在嚴重厚尾現象。另外，欲觀察在不同的 n 與 p 之下，各方法的估計表現，考慮五種 n 與 p 的組合： (n, p) = (100, 200), (200, 200), (500, 200), (100, 10), (100, 300)。. 政治大. 圖 3-1. 學. N(0, 1), dbexp(0, 1/√2), Cauchy(0, 1)機率分佈圖. ‧. ‧ 國. 立. 將三種εi 分配分別搭配五種 n 與 p 組合來模擬，看看不同情況下各方法估計. y. Nat. sit. 線性迴歸模型係數的表現。模擬的流程如下：. 2.. 隨機產生一筆數據：(xi1 , xi2 , … , xip , yi ), i = 1, 2, … , n。. 3.. 將數據隨機分成兩類 (各 n/2 個)：. n. al. er. 給定εi 分配、n 與 p。. io. 1.. 4.. Ch. a.. 訓練資料 (training data). b.. 測試資料 (testing data). engchi. i n U. v. 各方法使用訓練資料，分別採用 3 折 (3-fold) 交叉驗證，找出使預測均方差最小之調諧參數/超參數，作為最佳參數。. 5.. ̂ = (β̂1 , β̂2 , … , β̂p )。各方法給定最佳參數，使用測試資料求得估計係數 𝛃. 6.. 重複步驟 [3]～[5] 50 次，各方法皆得到 50 組估計係數 ̂ (𝟏) = (β̂1(1) , β̂(1) ̂ (1) ̂ (𝟐) = (β̂1(2) , β̂(2) ̂ (2) ̂ (𝟓𝟎) = (β̂1(50) , β̂(50) ̂ (50) 𝛃 2 , … , βp ), 𝛃 2 , … , βp ), … , 𝛃 2 , … , βp ). 7.. 各方法分別取出 50 組估計係數中每一係數的中位數 10.

(17) ̂ ∗ = (Median(β̂1(1) , β̂1(2) , … , β̂1(50) ), Median(β̂(1) ̂ (2) ̂ (50) ̂ (1) ̂ (2) ̂ (50) 𝛃 2 , β2 , … , β2 ), … , Median(βp , βp , … , βp )). 作為此組數據之最終估計係數。 8.. 模擬 100 次：重複步驟 [2]～[7] 100 次，各方法皆得到 100 組估計係數 ̂ ∗{1} , 𝛃 ̂ ∗{2} , … , 𝛃 ̂ ∗{100} 。 𝛃 依據各方法之 100 次模型估計係數計算評估指標，並作為評估各方法表現之依據 (表 3-3～3-5)。各指標的說明如表 3-1 和 3-2。. 立. 政治大. 學 ‧. ‧ 國 io. sit. y. Nat. n. al. er. 9.. Ch. engchi. 11. i n U. v.

(18) 表 3-1. 估計係數評估指標. 說明. TN：非 0 係數總個數. 100 次模擬中，平均每次非 0 係數總個數，即平均每次 TN. (number of total nonzero coefficients). 篩選變數總個數。(0 ≤ TN ≤ p). CN：正確非 0 係數個數. 100 次模擬中，平均每次正確非 0 係數個數，即平均每 CN. (number of correct nonzero coefficients). 次篩選正確變數個數。已知 p 個解釋變數中有 5 個為. +. 顯著變數，故 CN 理想值為 5。(0 ≤ CN ≤ 5). IN. IN：錯誤非 0 係數個數. 100 次模擬中，平均每次錯誤非 0 係數個數，即平均每. 政治大立估計係數”有”5 個正確非 0 係數且”無”錯誤非 0 係數的. (number of incorrect nonzero coefficients). CFR：正確配適比例. 表. (0% ≤ CFR ≤ 100%). 3-2. y. sit. er. 模擬次數佔總模擬次數 (100) 的比例。OFR 愈低愈好。. a l (0% ≤ OFR ≤ 100%) i v n 0 係數且”無”錯誤非 0 係數 Ch 估計係數”不足”5 個正確非 U engchi. n. (under fitting rate). ‧. ‧ 國. 學. io UFR：低度配適比例. 的模擬次數佔總模擬次數 (100) 的比例。UFR 愈低愈好。(0% ≤ UFR ≤ 100%). MSE：均方差 (mean squared error of coefficients). 詳. 數篩選表現」的主要準則，其值愈高愈好。. 估計係數”有”5 個正確非 0 係數且”有”錯誤非 0 係數的. Nat. (over fitting rate). 次篩選錯誤變數個數。IN 理想值為 0。(0 ≤ IN ≤ p − 5). 模擬次數佔總模擬次數 (100) 的比例。CFR 為評估「變見. (correct fitting rate). OFR：過度配適比例. 註. =. 評估指標. p. 2. 100 次估係數均方差的平均，即E[∑j=1(β̂j − βj ) ]。MSE 為評估「估計係數準確度」的主要準則，其值愈低愈好。. 12.

(19) 表 3-2. 正確配適、過度配適、低度配適定義. (β̂1 , β̂2 , … , β̂5 ) 皆為非 0 係數有 0 係數 (β̂6 , β̂7 , … , β̂p ) 有非 0 係數. 過度配適. 皆為 0 係數. 正確配適. 低度配適. 第三節模擬結果與討論表 3-3～3-5 為三種εi 分配分別搭配五種 n 與 p 組合的模擬結果。接下來，針. 政治大. 對「變數篩選表現」和「估計係數準確度」來探討，兩個面向分別以 CFR (correct. 立. fitting rate) 和 MSE (mean squared error of coefficients) 作為評估各方法優劣的主. ‧. ‧ 國. 學. io. sit. y. Nat. n. al. er. 要主要準則。. Ch. engchi. 13. i n U. v.

(20) 表 3-3. 立. εi ～N(0, 1)模擬結果. 政治大. ‧. ‧ 國. 學. io. sit. y. Nat. n. er. 14. al. Ch. engchi. i n U. v.

(21) 表 3-4. 立. εi ～dbexp(0, 1/√2)模擬結果. 政治大. ‧. ‧ 國. 學. io. sit. y. Nat. n. er. 15. al. Ch. engchi. i n U. v.

(22) 表 3-5. 立. εi ～Cauchy(0, 1)模擬結果. 政治大. ‧. ‧ 國. 學. io. sit. y. Nat. n. er. 16. al. Ch. engchi. i n U. v.

(23) 首先評估「變數篩選表現」。當εi ～N(0, 1)時 (表 3-3)，先考慮 n 固定為 100 的情況，不論 p 為多少，LASSO 的 CFR 皆最高，且為 100%，亦即 100 次模擬中，100 次皆正確選出 5 個顯著變數，且正確刪除非顯著變數，是所有方法中變數篩選表現最好的。LAD-LASSO 的 CFR 皆最低，是所有方法中篩選表現最差的。另外，可以看到僅有 EBLASSO 存在低度配適 (under-fitting) 的情形 (UFR > 0%且CN < 5)，也就是 EBLASSO 有選到較少顯著變數的傾向。接著考慮 p 固定為 200 的情況，當樣本數小 (n=100) 時，LASSO 的 CFR 最高， LAD-LASSO 最低；當樣本數大 (n=500) 時，EBLASSO 的 CFR 最高，EBENet. 政治大換句話說，當 n=500 時，EBENet 相較於其它方法，明顯有誤選非顯著變數的傾立. 最低，且此時 EBENet 的 OFR 為 100%，IN 為 7，兩項數據皆遠高於其它方法。. 向，且平均每次選到 7 個。顯然地，樣本數對各方法的變數篩選表現有很大的影. ‧ 國. 學. 響。. ‧. 當εi ～dbexp(0, 1/√2)時 (表 3-4)，先觀察 n 固定為 100 的情況，不論 p 為多. y. Nat. 少，LASSO 的 CFR 皆最高，是所有方法中變數篩選表現最好的。接著觀察 p 固. er. io. sit. 定為 200 的情況，當 n=100 時，LASSO 的 CFR 最高；當 n=500 時，EBLASSO 最高。比較表 3-3 與 3-4，可以發現兩張表格的結果十分相似。由於 dbexp(0, 1/√2). al. n. v i n 的厚尾現象並不明顯，因此結果與 C h N(0, 1)接近是可以預期的。 engchi U. 當εi ～Cauchy(0, 1)時 (表 3-5)，考慮 p 固定為 200 的情況，可以發現除了. EBLASSO，其餘方法的 CFR 皆是隨著 n 增加而減少，亦即樣本數愈大，愈不容易正確找出顯著變數。不論 n 和 p 為哪種組合，LAD-LASSO 的 CFR 皆比其它方法高出許多，也就是無論樣本數和變數個數為多少，LAD-LASSO 的變數篩選表現皆比其它方法優異許多。另外，可以發現 EBLASSO 皆存在嚴重低度配適的情形 (UFR 接近 100%、CN < 1且 IN 接近 0)，換句話說，此方法不太容易選到變數。由於 Cauchy(0, 1)存在明顯的厚尾現象，因此表 3-5 的結果與表 3-3 和 3-4 有相當大的差異。. 17.

(24) 接著評估「估計係數準確度」。當εi ～N(0, 1)時 (表 3-3)，先考慮 n 固定為 100 的情況，當 p=200 或 300 時，LAD-LASSO 的 MSE 最小，亦即 LAD-LASSO 的估計係數誤差最小，準確度最高。EBLASSO 的 MSE 最大，也就是 EBLASSO 的估計係數誤差最大，準確度最差。當 p=10 時，則是 EBENet 的 MSE 最小， EBLASSO 最大。接著考慮 p 固定為 200 的情況，當樣本數小 (n=100) 時， LAD-LASSO 的 MSE 最小，EBLASSO 最大；當樣本數大 (n=500) 時，EBLASSO 的 MSE 最小，EBENet 最大。由此可見，樣本數小時，EBLASSO 的估計係數準確度最差，但樣本數大時，EBLASSO 的準確度卻又最好。顯然地，當解釋變數. 政治大當ε ～dbexp(0, 1/√2)時 (表 3-4)，先觀察 n 固定為 100 的情況，不論 p 為多立. 個數固定時，樣本數對各方法的估計係數準確度有很大的影響。 i. 少，LAD-LASSO 的 MSE 皆最小，是所有方法中估計係數準確度最高的。接著. ‧ 國. 學. 觀察 p 固定為 200 的情況，EBLASSO 在 n=100 時 MSE 最大，在 n=500 時卻又. ‧. 最小，如同εi ～N(0, 1)的結果。由於 dbexp(0, 1/√2)的厚尾現象並不明顯，因此結. y. Nat. 果與 N(0, 1)相似是可以預期的。. er. io. sit. 當εi ～Cauchy(0, 1)時 (表 3-5)，不論 n 和 p 為哪種組合，LAD-LASSO 的 MSE 皆遠小於其它方法，亦即無論樣本數和變數個數為多少，LAD-LASSO 都. al. n. v i n 是所有方法中估計係數準確度最高的，且遠勝過其它方法。 Ch engchi U. 縱觀來看，當εi ～N(0, 1)或εi ～dbexp(0, 1/√2)，亦即反應變數無厚尾偏差或. 僅有輕微厚尾偏差時，沒有任何一個方法在不同的 n 與 p 組合都有最好的估計表現。然而，當εi ～Cauchy(0, 1)，也就是反應變數存在嚴重厚尾偏差時，不論 n 和 p 為哪種組合，LAD-LASSO 的「變數篩選表現」與「估計係數準確度」皆比其它方法優越許多。這樣的結果充分驗證了 LAD-LASSO 的穩健性 (robustness)，當反應變數存在明顯的厚尾偏差時，LAD-LASSO 不但變數篩選正確性高，係數也能估計得準確，其它方法則沒有很好的估計表現。. 18.

(25) 第四章變數分群模擬研究第一節前言在第二章提到當資料中存在幾個高度相關且顯著的解釋變數時，Elastic Net 和 EBENet 有篩選一整群的趨勢，且能夠達到變數分群。為了驗證這樣的特性，設計了以下的模擬情境，並藉由係數路徑圖，觀察變數間存在高度相關時，各方法變數篩選的特性。. 第二節模擬過程. 模擬假設. 政治大 𝐗 = (X , X , … , X ) 與 Y 滿足以下關係 (Zou and Hastie, 2005) 立 1. 2. 9. 學. ‧ 國. Y～N(Z1 + 0.5Z2 + 0.1Z3 , 1). X3 = Z1 + ε3 ,. X4 = Z2 + ε4 ,. X5 = −Z2 + ε5 ,. X6 = Z2 + ε6 ,. X7 = Z3 + ε7 ,. X8 = −Z3 + ε8 ,. X9 = Z3 + ε9 。. io. al. εj iid N(0, 1/16), j = 1, 2, … , 9。樣本數 n=100。此～. er. 其中Zi iid U(0, 20), i = 1, 2, 3, ～. sit. y. ‧. X2 = −Z1 + ε2 ,. Nat. X1 = Z1 + ε1 ,. v. n. 數據經由Z1 、Z2 、Z3 三個均勻分配 (uniform distribution) 生成 X 與 Y 的關係，. Ch. engchi. i n U. 其中X1 ～X9 分別為一個均勻分配加上一個隨機誤差項。藉由三個均勻分配，將 X 分成三個分群：(X1 , X 2 , X3 )、(X4 , X5 , X6)、(X7 , X 8 , X9 )，而 Y 服從常態分配，其期望值為三個均勻分配的線性組合。隨機產生一組數據(x, y)，觀察 x 的相關係數矩陣 (表 4-1)，可以發現同一群變數中，兩兩之間的相關係數都非常接近 1 或-1，也就是彼此間呈現高度正相關或負相關。不同群的變數，兩兩間相關係數則趨近於 0。接下來，以這組數據為樣本來繪製各方法的線性迴歸模型係數路徑圖，並觀察X1 ～X9 的估計係數走向。由於模擬設定不同分群對 Y 期望值影響的規模 (scale) 不同，故可以藉由係數路徑之間的重合或對稱關係，看出同一群的變數。 19.

(26) 表 4-1. X1 X1 X2. X2. X3. 1 -0.998 -0.998. x 的相關係數矩陣. X4. 0.998. X5. X6. 0.018 -0.021. 1 -0.998 -0.033. X7. X8. 0.031 -0.058. 0.037 -0.046. X9. 0.063 -0.057. 0.054 -0.059. 0.053. X3. 0.998 -0.998. 1. 0.026 -0.029. 0.038 -0.052. 0.057 -0.051. X4. 0.018 -0.033. 0.026. 1 -0.998. 0.998 -0.133. 0.134 -0.149. X5 X6 X7 X8 X9. -0.021. 0.037 -0.029 -0.998. 0.031 -0.046 -0.058. 0.038. 0.998 -0.998. 0.054 -0.052 -0.133. 0.063 -0.059 -0.057. 1 -0.998. 0.057. 0.127 -0.128. 1 -0.122. 0.127 -0.122. 0.134 -0.128. 0.053 -0.051 -0.149. 0.124 -0.139. 1 -0.997. 0.124 -0.997. 0.143 -0.139. 0.998. 1 -0.998. 0.998 -0.998. 政治大. 0.143. 1. 立模擬結果與討論第三節. ‧ 國. 學. 圖 4-1 ～ 4-5 為各方法的係數路徑圖，其中縱軸為 X1 ～X9 的估計係數. ‧. ̂ |, 即 L1-Norm) 為其橫軸的對應值。 (β̂1, β̂2 , … , β̂9 )，9 個值的絕對值總和 (∑9j=1 |β j. sit. y. Nat. 黑色垂直虛線為最佳調諧參數/超參數 (經由交叉驗證決定) 對應的係數解絕對. io. n. al. er. 值總和。可以發現係數路徑皆有由右至左逐漸收斂的趨勢。. Ch. 圖 4-1. engchi. i n U. LASSO 係數路徑圖 20. v.

(27) 圖 4-2. 立. 政治大 Elastic Net 係數路徑圖. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. 圖 4-3. engchi. i n U. LAD-LASSO 係數路徑圖. 21. v.

(28) 圖 4-4. 立. 政治大 EBLASSO 係數路徑圖. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. 圖 4-5. engchi. i n U. EBENet 係數路徑圖. 22. v.

(29) 圖 4-6. 立. 政治大. L1-Norm 與λ關係圖 (LASSO). ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. 圖 4-7. Ch. engchi. i n U. v. LASSO 係數路徑圖 (橫軸為λ). 23.

(30) 已知(X1 , X2 , X3 )、(X4 , X5 , X6)和(X7 , X8 , X9 )為三個分群。先看 Elastic Net 的係數路徑圖 (圖 4-2)，從圖中可以發現X1 和X3 幾乎是兩條重合的線，且和X2 對稱於橫軸，這表示不論λ如何變動，X1 和X3 的估計係數都非常接近，X2 則與X1 、X3 差個負號。(X4 , X5 , X6 )和(X7 , X8 , X9 )亦有相似的情形。簡言之，Elastic Net 可以從係數路徑之間的重合或對稱關係，清楚看出哪些變數屬於同一群。接著看到 EBENet 的係數路徑圖 (圖 4-5)，(X4 , X5 , X6 )和(X7 , X8 , X9 )這兩個分群都和 Elastic Net 的情況相似，係數路徑之間具有重合或對稱的關係。(X1 , X2 , X3 ) 則較為不同，可以看到X3 的估計係數值特別顯著，但X1 和X2 都僅在最右端出現. 政治大 10：5：1，規模不相同，而(X , X , X )是三個分群當中規模最大的，猜測是這樣立. 一小段，由右至左很快就縮減到 0。由於三個分群對 Y 期望值的影響為 1. 2. 3. 的原因造成與另外兩群的差異。. ‧ 國. 學. 再觀察 LASSO (圖 4-1)、LAD-LASSO (圖 4-3) 和 EBLASSO (圖 4-4) 三張. ‧. 係數路徑圖，可以發現三種方法皆無路徑之間的重合或對稱關係，也就沒辦法藉. y. Nat. 由係數路徑圖來看出哪些變數屬於同一群。LASSO 和 EBLASSO 都有同一群變. 僅在右端出現一小段，由右至左很快就縮減到 0。. n. al. Ch. engchi. 24. er. io. sit. 數中只選出一個的傾向。以 LASSO 為例，X2 的估計係數值特別顯著，X1 和X3 都. i n U. v.

(31) 第五章實際資料應用第一節攝護腺癌 (Prostate Cancer) 研究應用資料背景這筆資料使用於攝護腺癌的研究 (Stamey et al., 1989)，主要探討攝護腺特異抗原 (Prostate Specific Antigen, PSA) 與腫瘤體積、攝護腺重量等變數之間的關聯性。PSA 是由攝護腺上皮組織所分泌的特殊蛋白，為攝護腺癌的良好生物指標。. 政治大 log(prostate specific antigen, 攝護腺特異抗原) (lpsa)。解釋變數有 8 個，分別為立本資料的樣本為 97 位即將接受攝護腺切除手術的男性患者。反應變數為. log(cancer volume, 腫瘤體積) (lcavol)；log(prostate weight, 攝護腺重量) (lweight)；. ‧ 國. 學. age, 年齡 (age)；log(benign prostatic hyperplasia amount, 良性攝護腺增生量). ‧. (lbph)；seminal vesicle invasion, 儲精囊侵犯 (svi)；log(capsular penetration, 囊膜. y. Nat. 渗透) (lcp)；Gleason score, 格利森評分 (gleason) 和 percentage Gleason scores 4. er. io. sit. or 5, 格利森評分 4 或 5 的比例 (pgg45)。. 接下來，使用各方法估計線性迴歸模型係數，並比較好壞。資料操作的流程. n. al. 如下：. Ch. engchi. i n U. v. 1. 將 97 個樣本隨機分成兩類 (個數固定)： a. 訓練資料 (67 個) b. 測試資料 (30 個). 2. 各方法使用訓練資料尋找最佳調諧參數/超參數 (採用 3 折交叉驗證)，並求得 ̂ = (β̂0 , β̂1 , … , β̂8 )。模型估計係數 𝛃 3. 各方法使用測試資料代入模型做預測，並計算預測均方差： ̂i − Yi )2 , i = 1, 2, … , 30 (Y ̂i ：預測值，Yi ：真實值)。 MSE = E(Y 4. 重複步驟 [1]～[3] 50 次，各方法皆得到 50 組估計係數 25.

(32) ̂ (𝟏) = (β̂1(1) , β̂(1) ̂ (1) ̂ (𝟐) = (β̂1(2) , β̂(2) ̂ (2) ̂ (𝟓𝟎) = (β̂1(50) , β̂(50) ̂ (50) 𝛃 2 , … , β8 ), 𝛃 2 , … , β8 ), … , 𝛃 2 , … , β8 ). 與 50 個預測均方差MSE (1) , MSE (2) , … , MSE (50) 。 5. 各方法分別取出 50 組估計係數中每一係數的中位數 ̂ ∗ = (Median(β̂1(1) , β̂1(2) , … , β̂1(50) ), Median(β̂(1) ̂ (2) ̂ (50) ̂ (1) ̂ (2) ̂ (50) 𝛃 2 , β2 , … , β2 ), … , Median(β8 , β8 , … , β8 )). 作為最終估計係數，成為最終的模型。觀察各方法篩選的變數 (係數不為 0 者) 與係數值 (表 5-1)。 6. 計算各方法之 50 個預測均方差的平均： Avg. MSE = Mean(MSE (1) , MSE (2) , … , MSE (50) ) ，並以此作為比較各方法優劣的主要準則 (表 5-2)。. 立. ‧ 國. 學. 結果與討論. 政治大. 由表 5-1，可以發現 Elastic Net 的估計係數中有 6 個為非 0 (不含截距項)，. ‧. 亦即選到 6 個變數，是所有方法中最多的，而 EBLASSO 只選到 1 個變數則是最. sit. y. Nat. 少的。觀察各個變數的選取結果，可以發現 lcavol (log 腫瘤體積) 是所有方法都. al. er. io. 有選到的，可能是影響 PSA 最重要的變數。lcp (log 囊膜渗透)、gleason (格利森. n. 評分) 則完全沒有被任何方法選中。. Ch. 表 5-1. 方法＼變數. (intercept). lcavol. engchi. i n U. v. 攝護腺癌資料：各方法估計係數. lweight. age. lbph. svi. lcp. gleason. pgg45. LASSO. 0.400. 0.525. 0.390. 0. 0.060. 0.625. 0. 0. 0.002. Elastic Net. 0.333. 0.503. 0.396. -0.004. 0.046. 0.592. 0. 0. 0.002. LAD-LASSO. 1.865. 0.255. 0. 0. 0. 0. 0. 0. 0.010. EBLASSO. 1.632. 0.638. 0. 0. 0. 0. 0. 0. 0. EBENet. 1.572. 0.574. 0. 0. 0.106. 0.601. 0. 0. 0. 26.

(33) 表 5-2. 攝護腺癌資料：各方法平均預測均方差，括弧內為標準差. 方法. Avg. MSE. LASSO. 0.614 (0.152). Elastic Net. 0.612 (0.151). LAD-LASSO 0.926 (0.294) EBLASSO. 0.825 (0.231). EBENet. 0.597 (0.146). 政治大資料的方法。其估計係數中有三個為非 0 (表 5-1)，亦即選到 3 個變數，分別是立根據表 5-2，EBENet 為所有方法中平均預測均方差最小的，是最適合此筆. ‧ 國. 學. lcavol (log 腫瘤體積)、lbph (log 良性攝護腺增生量)、svi (儲精囊侵犯)，這三個變數可能是影響 PSA 的重要變數。. ‧. PSA 雖然是攝護腺癌的良好生物指標，但 PSA 上升並不一定代表癌症。其. sit. y. Nat. 它如攝護腺良性增生、感染、近期內有做過攝護腺超音波/切片或長時間騎乘腳. io. er. 踏車等攝護腺相關活動都可能造成 PSA 上升 (黃書彬, 2014)。. al. 本資料的解釋變數皆是與攝護腺癌息息相關的醫學檢驗數據 (年齡除外)。. n. v i n Ch 以估計誤差最小的 EBENet 篩選的變數為例，三個非 e n g c h i U 0 係數皆為正值，也就是對. PSA 有正向的影響。腫瘤體積 (cancer volume) 愈大，則 PSA 愈大；良性攝護腺增生量 (benign prostatic hyperplasia amount) 愈大，則 PSA 也愈大。腫瘤體積增大和良性攝護腺增生確實都會造成 PSA 上升，如上一段提到的，PSA 是攝護腺癌的指標，但其值上升不一定是癌症。儲精囊侵犯 (seminal vesicle invasion) 為腫瘤擴散至攝護腺包膜外，侵犯到儲精囊，是攝護腺癌晚期的症狀。此症狀也確實會造成 PSA 上升。而格利森評分 (Gleason score) 雖然沒有任何方法篩選此變數，但它卻是一種診斷攝護腺癌的重要指標。由於資料中的格利森評分大多是 6 或 7，彼此間差異小，可能因此造成這個變數未被任何方法選取。 27.

(34) 第二節白血病 (Leukemia) 研究應用資料背景這筆資料使用於白血病的研究 (Golub et al., 1999)，為一微陣列資料，主要探討白血病類型與各基因之間的關聯性。微陣列資料通常有上千個基因和小於 100 個樣本，搜集自微陣列實驗。微陣列實驗能同時大量地、全面性地偵測上萬個基因表現量，透過基因晶片，可在短時間內找出可能受疾病影響的基因，作為早期診斷的生物指標 (biomarker) (蔡政安, 2009)。本資料的樣本為 72 位白血病患者。白血病分為兩類，分別是急性淋巴細胞. 政治大. 白血病 (Acute Lymphoblastic Leukemia, ALL) 與急性骨髓性白血病 (Acute. 立. Myeloid Leukemia, AML)。72 位患者中有 47 位為 ALL、25 位為 AML。反應變. ‧ 國. 學. 數為白血病類型 (1=ALL, 0=AML)，解釋變數為 7129 個基因。由於反應變數為二分 (dichotomous) 的型態，目標模型為羅吉斯迴歸。. ‧. 接下來，使用各方法估計羅吉斯迴歸模型係數，並比較好壞。資料操作的流. sit. y. Nat. 程如下：. n. al. er. io. 1. 將 72 個樣本隨機分成兩類 (ALL 和 AML 的個數固定)：. i n U. a. 訓練資料 (48 個)：31 個 ALL、17 個 AML. Ch. engchi. b. 測試資料 (24 個)：16 個 ALL、8 個 AML. v. 2. 各方法使用訓練資料尋找最佳調諧參數/超參數 (採用 3 折交叉驗證)，並求得 ̂ = (β̂0 , β̂1 , … , β̂7129 )。統計各方法篩選的基因 (係數不為 0 模型估計係數 𝛃 者) 與篩選個數n′ 。 3. 各方法使用測試資料代入模型做類型預測，預測機率大於或等於 0.5 者為 ALL，小於 0.5 者為 AML。 4. 計算各方法預測正確率： CR = (測試資料預測正確數/測試資料總個數(24)) × 100%。 ̂ (𝟏) , 𝛃 ̂ (𝟐) , … , 𝛃 ̂ (𝟓𝟎)、50 5. 重複步驟 [1]～[4] 50 次，各方法皆得到 50 組估計係數𝛃 28.

(35) 個篩選個數n′. (1). , n′. (2). , … , n′. (50). 與 50 個預測正確率CR(1) , CR(2) , … , CR(50) 。. 6. 統計各方法平均每次篩選個數： Avg. n′ = Mean(n′. (1). , n′. (2). , … , n′. (50). )，並做比較 (表 5-3)。. 7. 統計各方法所有篩選次數至少一半 (25 次) 的基因，並做比較 (表 5-4)。 8. 計算各方法之 50 個預測正確率的平均： Avg. CR = Mean(CR(1) , CR(2) , … , CR(50) )，並以此作為比較各方法優劣的主要準則 (表 5-5)。 [註]：由於使用上的限制，此資料忽略 LAD-LASSO 這個方法。. 立. 結果與討論. 政治大. ‧ 國. 學. 根據表 5-3，Elastic Net 平均每次篩選 83.68 個基因，是所有方法中最多的，. Nat. 白血病資料：各方法平均基因篩選個數. io. n. al. Avg. 𝐧′. 方法. C LASSO h e n g c19.28 hi Elastic Net. 83.68. EBLASSO. 2.08. EBENet. 8.68. 29. er. 表 5-3. sit. y. ‧. 而 EBLASSO 只篩選 2.08 個則是最少的。. i n U. v.

(36) 表 5-4. 白血病資料：各方法基因篩選與次數. LASSO (8 個) Y07604_at. M19507_at. X95735_at. M23197_at. X17042_at. M13690_s_at. M31994_at. M84526_at. 44. 43. 40. 34. 29. 28. 25. 25. Elastic Net (46 個) M19507_at. M23197_at. M27891_at. M84526_at. X95735_at. Y07604_at. M63138_at. U46499_at. 50. 50. 50. 50. 50. 50. 47. 46. U82759_at. X17042_at. M92287_at. M96326_rna1_at. X59417_at. M13690_s_at. X85116_rna1_s_at. M84371_rna1_s_at. 46. 46. 45. 45. 45. 44. M16038_at. M22960_at. D49950_at. M34344_at. L09209_s_at. 42. 42. 39. 39. M11147_at. Z29067_at. 38. 37. X51521_at. HG2562-HT2658_s_at. M26708_s_at. 33. 33. 33. M20902_at. D42043_at. M65214_s_at. 28. 27. 27. 立. 政 45治大45. M31994_at. 41. 41. 39. 39. HG1612-HT1612_at. M62762_at. M55150_at. Z24727_at. M31303_rna1_at. Y08612_at. 36. 36. 35. 35. 34. 34. M31211_s_at. M11722_at. U05259_rna1_at. U59632_s_at. M28170_at. 33. 30. 30. 29. 29. X63097_at. Ch. S50223_at. 26. 25. n. EBLASSO (0 個) 無 EBENet (1 個) X95735_at. 34 [註]：各基因名稱對應的數字為篩選次數 (至少 25 次)。. 30. y. sit. io. engchi. er. Nat. al. ‧. ‧ 國. J05243_at. 學. U50136_rna1_at. iv X16546_at n U 25.

(37) 將表 5-3 與 5-4 相比對，Elastic Net 平均每次能選到約 84 個基因，篩選次數至少一半的有 46 個。Elastic Net 篩選的基因中，M19507_at、M23197_at、 M27891_at、M84526_at、X95735_at 和 Y07604_at 這 6 個基因的篩選次數達 50 次，亦即平均每次選到的 84 個基因中，這 6 個是每次都能選到的，其餘的基因則每次篩選不盡相同，但共有 46 個篩選次數達到一半 (25 次)。EBLASSO 平均每次能選到約 2 個基因，但每次篩選皆不太一致，所以沒有任何基因篩選次數至少一半。由表 5-4，可以觀察到 LASSO 篩選次數至少一半的 8 個基因皆包含在 Elastic. 政治大也包含在 LASSO 和 Elastic Net 當中。由此可見，各方法篩選的基因有重複的現立. Net 的 46 個基因裡。另外，EBENet 只有 X95735_at 篩選次數達到一半，此基因. 學. 表 5-5. 白血病資料：各方法平均預測正確率. LASSO. 93.7%. n. al. Elastic Net. 96.8%. C EBLASSO h e n g c89.1% hi. sit er. io. Avg. CR. y. Nat. 方法. ‧. ‧ 國. 象。. i n U. v. 91.0%. EBENet. 根據表 5-5，Elastic Net 為所有方法中平均預測正確率最高的，是最適合此筆資料的方法，其平均正確率為 96.8%，亦即每次的 24 個測試資料中，平均有 23.2 個分類正確。此方法有 46 個基因篩選次數至少一半，這些基因可能是影響白血病類型的重要變數。LASSO 的預測正確率為 93.7%，僅次於 Elastic Net。為了驗證各方法篩選的基因 (表 5-4) 是否為好的變數選擇，接下來將選取幾種基因選擇，並透過主成分分析 (Principal Component Analysis, PCA) 來繪製 31.

(38) 雙標圖 (Biplot)，觀察樣本的分佈狀況及各基因的影響情形。進行 PCA 前需先將原始資料數值標準化。由於 EBLASSO 和 EBENet 分別只選出 0 個和 1 個基因，無法進行 PCA，故忽略這兩種方法的結果。. 表 5-6. 基因選擇一：LASSO (篩選次數至少 25 次). Y07604_at. M19507_at. X95735_at. M23197_at. X17042_at. M13690_s_at. M31994_at. M84526_at. 44. 43. 40. 34. 29. 28. 25. 25. 政治大縱軸分別為第一主成分 (PC1) 與第二主成分 (PC2)，其中 PC1 能解釋 47.6%的立. 利用基因選擇一 (表 5-6) 的 8 個基因所繪製的雙標圖如圖 5-1 所示，橫軸與. 變異，PC2 能解釋 14.7%的變異，共 62.3% (表 5-7)。從圖中可以發現 ALL 群聚. ‧ 國. 學. 的情形相當明顯，AML 則較為分散。另外，兩種類型的分佈能夠以PC1 = 0作為. ‧. 分界線，此線以左的樣本為 ALL，以右為 AML，代表只需透過 PC1 便能區別兩. sit. y. Nat. 種類型的白血病。. io. er. 觀察樣本的分佈與各基因的指向，可以發現 ALL 皆分佈在這 8 個基因相反指向的位置，代表這些基因對 ALL 有反向的影響。一個樣本中的這 8 個基因表. al. n. v i n Ch 現量值 (標準化後) 整體來看若愈小，愈有可能歸類在 ALL。 engchi U. 32.

(39) 政治大. 立. n. al. PC1. 白血病資料：各主成分資訊 (I). PC2. Ch. PC3. i e n0.980 g c h0.830. 1.083. PC5. v. PC6. PC7. PC8. 0.651. 0.637. 0.583. 0.443. i n U. PC4. sit. io. 表 5-7. er. Nat. y. ‧. ‧ 國. 白血病資料：雙標圖 (I). 學. 圖 5-1. Standard deviation. 1.952. Proportion of Variance. 0.476. 0.147. 0.120. 0.086. 0.053. 0.051. 0.042. 0.025. Cumulative Proportion. 0.476. 0.623. 0.743. 0.829. 0.882. 0.933. 0.975. 1.000. 33.

(40) 表 5-8. 基因選擇二：Elastic Net (篩選次數至少 25 次). M19507_at. M23197_at. M27891_at. M84526_at. X95735_at. Y07604_at. M63138_at. U46499_at. 50. 50. 50. 50. 50. 50. 47. 46. U82759_at. X17042_at. M92287_at. M96326_rna1_at. X59417_at. M13690_s_at. X85116_rna1_s_at. M84371_rna1_s_at. 46. 46. 45. 45. 45. 45. 45. 44. M16038_at. M22960_at. D49950_at. U50136_rna1_at. J05243_at. M31994_at. M34344_at. L09209_s_at. 42. 42. 41. 41. 39. 39. 39. 39. M11147_at. Z29067_at. HG1612-HT1612_at. M62762_at. M55150_at. Z24727_at. M31303_rna1_at. Y08612_at. 38. 37. 36. 36. 34. 34. X51521_at. HG2562-HT2658_s_at. M26708_s_at. M31211_s_at. M11722_at. U05259_rna1_at. U59632_s_at. M28170_at. 33. 33. 33. 33. 30. 30. 29. 29. M20902_at. D42043_at. M65214_s_at. X63097_at. S50223_at. X16546_at. 28. 27. 27. 26. 25. 25. 學. Nat. sit. y. ‧. ‧ 國. 立. 政 35治大35. io. er. 利用基因選擇二 (表 5-8) 的 46 個基因所繪製的雙標圖如圖 5-2 所示，其中. al. PC1 能解釋 39.7%的變異，PC2 能解釋 8.8%的變異，共 48.5% (表 5-9)，解釋力. n. v i n Ch U 5-7 只擷取前 16 個)，因稍差。由於基因個數過多，主成分個數也多 e n g c h(46i 個，表此造成 2 維的雙標圖解釋力較低。然而，從圖中可以發現 ALL 和 AML 的群聚情形都相當明顯，且同樣只需透過 PC1 便能區別兩種類型的白血病。. 34.

(41) 政治大. 立. ‧. ‧ 國. 白血病資料：雙標圖 (II). 學. 圖 5-2. n. PC2. Ch. PC3. PC5. v. PC6. PC7. PC8. 1.324. 1.256. 1.248. 1.159. i n U. PC4. i e n1.767 g c h1.527. 2.012. er. io. al. PC1. 白血病資料：各主成分資訊 (II). sit. y. Nat. 表 5-9. Standard deviation. 4.271. Proportion of Variance. 0.397. 0.088. 0.068. 0.051. 0.038. 0.034. 0.034. 0.029. Cumulative Proportion. 0.397. 0.485. 0.552. 0.603. 0.641. 0.676. 0.709. 0.739. PC9. PC10. PC11. PC12. PC13. PC14. PC15. PC16. Standard deviation. 1.074. 1.004. 0.944. 0.911. 0.896. 0.883. 0.785. 0.747. Proportion of Variance. 0.025. 0.022. 0.019. 0.018. 0.017. 0.017. 0.013. 0.012. Cumulative Proportion. 0.764. 0.786. 0.805. 0.823. 0.841. 0.857. 0.871. 0.883. 35.

(42) 表 5-10. 基因選擇三：Elastic Net (篩選次數 50 次). M19507_at. M23197_at. M27891_at. M84526_at. X95735_at. Y07604_at. 50. 50. 50. 50. 50. 50. 由於 Elastic Net 篩選次數至少 25 次的基因個數過多，這次只選取其篩選次數達 50 次的基因，作為基因選擇三 (表 5-10)。利用這 6 個基因所繪製的雙標圖如圖 5-3 所示，其中 PC1 能解釋 60.2%的變異，PC2 能解釋 18.4%的變異，共 78.6% (表 5-11)。雖然 AML 的分佈較為分散，但同樣只需透過 PC1 便能區別兩種類型的白血病。. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. 36. i n U. v.

(43) 政治大. 立. ‧. ‧ 國. 白血病資料：雙標圖 (III). 學. 圖 5-3. n Standard deviation. 白血病資料：各主成分資訊 (III). PC1. Ch. PC2. i n U. PC3. i e n1.052 g c h0.698. 1.900. er. io. al. sit. y. Nat. 表 5-11. v. PC4. PC5. PC6. 0.641. 0.494. 0.375. Proportion of Variance. 0.602. 0.184. 0.081. 0.069. 0.041. 0.023. Cumulative Proportion. 0.602. 0.786. 0.867. 0.936. 0.977. 1.000. 37.

(44) 接下來，為了再次驗證這三種基因選擇是否為好的變數選擇，將它們直接當作解釋變數建立模型 (不需再篩選變數)。仿照資料背景中所述的操作流程，重複進行模型建立與類型預測，並計算平均預測正確率 (表 5-12)。. 表 5-12. 白血病資料：各基因選擇個數與模型平均預測正確率. 基因選擇. 個數. Avg. CR. 基因選擇一：LASSO (至少 25 次). 8個. 96.9%. 基因選擇二：Elastic Net (至少 25 次). 46 個. 79.0%. 個政治 6大. 基因選擇三：Elastic Net (50 次). 立. 93.8%. ‧ 國. 學. 雖然 Elastic Net 建立的模型有 96.8%的平均預測正確率 (表 5-5)，但將其篩選次數至少一半的 46 個基因 (基因選擇二) 作為解釋變數重新建立模型，卻只. ‧. 有 79.0%的正確率 (表 5-12)，是所有基因選擇中最低的。已知當資料中存在幾個. sit. y. Nat. 高度相關且顯著的解釋變數時，Elastic Net 傾向於篩選一整群。這 46 個基因當. io. er. 中，可能有部分存在高度相關，若將它們同時作為解釋變數，模型會產生共線性問題，導致預測值變異大，因此造成估計失準。顯然地，將這 46 個基因同時作. al. n. v i n C hLASSO 篩選次數至少一半的為解釋變數並非好的做法。使用 8 個基因 (基因選 engchi U 擇一) 作為解釋變數重新建模，有 96.9%的平均預測正確率，是所有基因選擇中. 最高的。縱觀來看，基因選擇一的 8 個基因是最好的變數選擇，亦即全部的 7129 個基因裡，只要選取這 8 個，便能夠完善解釋白血病類型。. 38.

(45) 第六章結論本論文比較了 LASSO 與幾種其衍生方法的特性，包含 LASSO、Elastic Net、 LAD-LASSO、EBLASSO 和 EBENet，這些方法被廣泛應用在線性模型係數的估計，特別是在擁有許多解釋變數且需要篩選部分出來的資料。本文中各方法的調諧參數/超參數皆是透過交叉驗證來尋找最佳值，不過除了運算求得，也可以主觀調整其值。若某種方法估計模型係數後，得到過多的非 0 係數，也就是選取過多的變數，沒有達到理想的篩選，可以逐次將λ調大，重新估計模型係數，直到篩選個數滿意為止。如同研究方法中提到的，λ愈大則估. 政治大模擬結果顯示當反應變數存在明顯的厚尾偏差時，LAD-LASSO 的估計表現立. 計係數縮減程度愈大，變數也就選得愈少。. ‧ 國. 學. 遠優於其它方法，充分展現了其穩健性。另外，當資料中存在幾個高度相關且顯著的解釋變數時，Elastic Net 和 EBENet 有篩選一整群的趨勢。. ‧. 每種估計模型係數方法都沒有絕對的好壞，一筆數據欲找出最適當的方法和. sit. y. Nat. 最理想的結果，可以重複操作資料。攝護腺資料針對同一筆數據重複分配訓練資. al. er. io. 料/測試資料、估計模型係數和預測，得到 EBENet 為最適合的方法。白血病資. v. n. 料也是針對同一筆資料重複運作，得到 Elastic Net 為最適合的方法。不過以基因. Ch. engchi. i n U. 選擇的觀點來看，Elastic Net 篩選的個數過多且可能存在共線性，並非最好的選擇，因此再搭配 PCA、雙標圖等，找出最好的基因選擇。簡言之，沒有一種估計方法使用在所有資料都是最好的，這當中也關係到分析者想要得到的結果。. 39.

(46) 參考文獻黃書彬，攝護腺特異抗原(PSA)過高的意義？？，上網日期 106 年 5 月 17 日，檢自 http://www.kmuh.org.tw/www/kmcj/data/10306/11.htm 蔡政安，2009。《微陣列資料分析(Microarray Data Analysis)》。中國醫藥大學生物統計中心。 Cai, X., Huang, A. and Xu, S. (2011). Fast empirical Bayesian LASSO for multiple quantitative trait locus mapping. BMC Bioinformatics, 12, 211. Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regression.. 政治大 Gao, X.L. and Huang, J. (2010). 立 Asymptotic analysis of high-dimensional LAD Ann. Statist., 32, 407-499.. ‧ 國. 學. regression with Lasso. Statistica Sinica, 20, 1485-1506.. Gill, P., Murray, W. and Wright, M., (1981). Practical optimization. New York:. ‧. Academic Press.. sit. y. Nat. Hoerl, A. and Kennard, R. (1988). Ridge regression. Encyclopedia of Statistical. n. al. er. io. Sciences, 8, 129-136.. i n U. v. Huang, A., Xu, S. and Cai, X. (2015). Empirical Bayesian elastic net for multiple. Ch. engchi. quantitative trait locus mapping. Heredity, 114, 107-115. Tibshirani, R. (1996). Regression Shrinkage and Selection via the Lasso. J. R. Statist. Soc. B, 58, 267-288. Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net. J. R. Statist. Soc. B, 67, 301-320.. 40.

(47)