r語言中對LASSO回歸，Ridge嶺回歸和彈性網絡Elastic Net模型實現
原文鏈接：http://tecdat.cn/?p=3795
Glmnet是一個通過懲罰最大似然關系擬合廣義線性模型的軟件包。正則化路徑是針對正則化參數λ的值網格處的lasso或Elastic Net（彈性網絡）懲罰值計算的。該算法非?？?，并且可以利用輸入矩陣中的稀疏性?x。它適合線性，邏輯和多項式，泊松和Cox回歸模型?？梢詮臄M合模型中做出各種預測。它也可以擬合多元線性回歸。
glmnet?解決以下問題
在覆蓋整個范圍的λ值網格上。這里l（y，η）是觀察i的負對數似然貢獻；例如對于高斯分布是
。?彈性網絡懲罰由α控制，LASSO（α= 1，默認），Ridge（α= 0）。調整參數λ控制懲罰的總強度。
眾所周知，嶺懲罰使相關預測因子的系數彼此縮小，而套索傾向于選擇其中一個而丟棄其他預測因子。彈性網絡則將這兩者混合在一起。
?glmnet?算法使用循環(huán)坐標下降法，該方法在每個參數固定不變的情況下連續(xù)優(yōu)化目標函數，并反復循環(huán)直到收斂，我們的算法可以非常快速地計算求解路徑。
代碼可以處理稀疏的輸入矩陣格式，以及系數的范圍約束，還包括用于預測和繪圖的方法，以及執(zhí)行K折交叉驗證的功能。
?
快速開始
?
首先，我們加載?glmnet?包：
library(glmnet)
包中使用的默認模型是高斯線性模型或“最小二乘”。我們加載一組預先創(chuàng)建的數據以進行說明。用戶可以加載自己的數據，也可以使用工作空間中保存的數據。
該命令?從此保存的R數據中加載輸入矩陣?x?和因向量?y。
我們擬合模型?glmnet。
fit = glmnet(x, y)
可以通過執(zhí)行plot?函數來可視化系數?：
plot(fit)
每條曲線對應一個變量。它顯示了當λ變化時，其系數相對于整個系數向量的?1范數的路徑。上方的軸表示當前λ處非零系數的數量，這是套索的有效自由度（df）。用戶可能還希望對曲線進行注釋。這可以通過label = TRUE?在plot命令中進行設置來完成?。
glmnet?如果我們只是輸入對象名稱或使用print?函數，則會顯示每個步驟的路徑?摘要?：
print(fit)
#### Call: ?glmnet(x = x, y = y)#### ? ? ? Df ? %Dev ?Lambda## ?[1,] ?0 0.0000 1.63000## ?[2,] ?2 0.0553 1.49000## ?[3,] ?2 0.1460 1.35000## ?[4,] ?2 0.2210 1.23000## ?[5,] ?2 0.2840 1.12000## ?[6,] ?2 0.3350 1.02000## ?[7,] ?4 0.3900 0.93300## ?[8,] ?5 0.4560 0.85000## ?[9,] ?5 0.5150 0.77500## [10,] ?6 0.5740 0.70600## [11,] ?6 0.6260 0.64300## [12,] ?6 0.6690 0.58600## [13,] ?6 0.7050 0.53400## [14,] ?6 0.7340 0.48700## [15,] ?7 0.7620 0.44300## [16,] ?7 0.7860 0.40400## [17,] ?7 0.8050 0.36800## [18,] ?7 0.8220 0.33500## [19,] ?7 0.8350 0.30600## [20,] ?7 0.8460 0.27800
它從左到右顯示了非零系數的數量（Df），解釋的（零）偏差百分比（%dev）和λ（Lambda）的值。
我們可以在序列范圍內獲得一個或多個λ處的實際系數：
coef(fit,s=0.1)
## 21 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ? 1## (Intercept) ?0.150928## V1 ? ? ? ? ? 1.320597## V2 ? ? ? ? ? .## V3 ? ? ? ? ? 0.675110## V4 ? ? ? ? ? .## V5 ? ? ? ? ?-0.817412## V6 ? ? ? ? ? 0.521437## V7 ? ? ? ? ? 0.004829## V8 ? ? ? ? ? 0.319416## V9 ? ? ? ? ? .## V10 ? ? ? ? ?.## V11 ? ? ? ? ?0.142499## V12 ? ? ? ? ?.## V13 ? ? ? ? ?.## V14 ? ? ? ? -1.059979## V15 ? ? ? ? ?.## V16 ? ? ? ? ?.## V17 ? ? ? ? ?.## V18 ? ? ? ? ?.## V19 ? ? ? ? ?.## V20 ? ? ? ? -1.021874
還可以使用新的輸入數據在特定的λ處進行預測：
predict(fit,newx=nx,s=c(0.1,0.05))
## ? ? ? ? ? ? 1 ? ? ? 2## ?[1,] ?4.4641 ?4.7001## ?[2,] ?1.7509 ?1.8513## ?[3,] ?4.5207 ?4.6512## ?[4,] -0.6184 -0.6764## ?[5,] ?1.7302 ?1.8451## ?[6,] ?0.3565 ?0.3512## ?[7,] ?0.2881 ?0.2662## ?[8,] ?2.7776 ?2.8209## ?[9,] -3.7016 -3.7773## [10,] ?1.1546 ?1.1067
該函數?glmnet?返回一系列模型供用戶選擇。交叉驗證可能是該任務最簡單，使用最廣泛的方法。
cv.glmnet?是交叉驗證的主要函數。
cv.glmnet?返回一個?cv.glmnet?對象，此處為“ cvfit”，其中包含交叉驗證擬合的所有成分的列表。
我們可以繪制對象。
它包括交叉驗證曲線（紅色虛線）和沿λ序列的上下標準偏差曲線（誤差線）。垂直虛線表示兩個選定的λ。
我們可以查看所選的λ和相應的系數。例如，
cvfit$lambda.min## [1] 0.08307
lambda.min?是給出最小平均交叉驗證誤差的λ值。保存的另一個λ是?lambda.1se，它給出了的模型，使得誤差在最小值的一個標準誤差以內。我們只需要更換?lambda.min?到lambda.1se?以上。
coef(cvfit, s = "lambda.min")
## 21 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ?1## (Intercept) ?0.14936## V1 ? ? ? ? ? 1.32975## V2 ? ? ? ? ? .## V3 ? ? ? ? ? 0.69096## V4 ? ? ? ? ? .## V5 ? ? ? ? ?-0.83123## V6 ? ? ? ? ? 0.53670## V7 ? ? ? ? ? 0.02005## V8 ? ? ? ? ? 0.33194## V9 ? ? ? ? ? .## V10 ? ? ? ? ?.## V11 ? ? ? ? ?0.16239## V12 ? ? ? ? ?.## V13 ? ? ? ? ?.## V14 ? ? ? ? -1.07081## V15 ? ? ? ? ?.## V16 ? ? ? ? ?.## V17 ? ? ? ? ?.## V18 ? ? ? ? ?.## V19 ? ? ? ? ?.## V20 ? ? ? ? -1.04341
注意，系數以稀疏矩陣格式表示。原因是沿著正則化路徑的解通常是稀疏的，因此使用稀疏格式在時間和空間上更為有效。
可以根據擬合的cv.glmnet?對象進行預測?。讓我們看一個示例。
## ? ? ? ? ? ?1## [1,] -1.3647## [2,] ?2.5686## [3,] ?0.5706## [4,] ?1.9682## [5,] ?1.4964
newx?與新的輸入矩陣?s相同，如前所述，是預測的λ值。
?
線性回歸
這里的線性回歸是指兩個模型系列。一個是?gaussian正態(tài)分布，另一個是?mgaussian多元正態(tài)分布。
正態(tài)分布
假設我們有觀測值xi∈Rp并且yi∈R，i = 1，...，N。目標函數是
其中λ≥0是復雜度參數，0≤α≤1在嶺回歸（α=0）和套索LASSO（α=1）之間。
應用坐標下降法解決該問題。具體地說，通過計算βj=β?j處的梯度和簡單的演算，更新為
其中
。
當x?變量標準化為具有單位方差（默認值）時，以上公式適用?。
glmnet?提供各種選項供用戶自定義。我們在這里介紹一些常用的選項，它們可以在glmnet?函數中指定?。
alpha?表示彈性網混合參數α，范圍α∈[0,1]。α=1是套索（默認），α=0是Ridge。weights?用于觀察權重。每個觀察值的默認值為1。nlambda?是序列中λ值的數量。默認值為100。lambda?可以提供，但通常不提供，程序會構建一個序列。自動生成時，λ序列由lambda.max?和?確定?lambda.min.ratio。standardize?是x?在擬合模型序列之前進行變量標準化的邏輯標志?。
例如，我們設置α=0.2，并對后半部分的觀測值賦予兩倍的權重。為了避免在此處顯示太長時間，我們將其設置?nlambda?為20。但是，實際上，建議將λ的數量設置為100（默認值）或更多。
然后我們可以輸出glmnet?對象。
print(fit)
#### Call: ?glmnet(x = x, y = y, weights = c(rep(1, 50), rep(2, 50)), alpha = 0.2, ? ? ?nlambda = 20)#### ? ? ? Df ?%Dev ?Lambda## ?[1,] ?0 0.000 7.94000## ?[2,] ?4 0.179 4.89000## ?[3,] ?7 0.444 3.01000## ?[4,] ?7 0.657 1.85000## ?[5,] ?8 0.785 1.14000## ?[6,] ?9 0.854 0.70300## ?[7,] 10 0.887 0.43300## ?[8,] 11 0.902 0.26700## ?[9,] 14 0.910 0.16400## [10,] 17 0.914 0.10100## [11,] 17 0.915 0.06230## [12,] 17 0.916 0.03840## [13,] 19 0.916 0.02360## [14,] 20 0.916 0.01460## [15,] 20 0.916 0.00896## [16,] 20 0.916 0.00552## [17,] 20 0.916 0.00340
這將顯示生成對象的調用?fit?以及帶有列Df?（非零系數的數量），??%dev?（解釋的偏差百分比）和Lambda?（對應的λ值）?的三列矩陣?。
我們可以繪制擬合的對象。
讓我們針對log-lambda值標記每個曲線來繪制“擬合”。
這是訓練數據中的偏差百分比。我們在這里看到的是，在路徑末端時，該值變化不大，但是系數有點“膨脹”。這使我們可以將注意力集中在重要的擬合部分上。
我們可以提取系數并在某些特定值的情況下進行預測。兩種常用的選項是：
s?指定進行提取的λ值。exact?指示是否需要系數的精確值。
一個簡單的例子是：
## 21 x 2 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ?1 ? ? ? ? 1## (Intercept) ?0.19657 ?0.199099## V1 ? ? ? ? ? 1.17496 ?1.174650## V2 ? ? ? ? ? . ? ? ? ?.## V3 ? ? ? ? ? 0.52934 ?0.531935## V4 ? ? ? ? ? . ? ? ? ?.## V5 ? ? ? ? ?-0.76126 -0.760959## V6 ? ? ? ? ? 0.46627 ?0.468209## V7 ? ? ? ? ? 0.06148 ?0.061927## V8 ? ? ? ? ? 0.38049 ?0.380301## V9 ? ? ? ? ? . ? ? ? ?.## V10 ? ? ? ? ?. ? ? ? ?.## V11 ? ? ? ? ?0.14214 ?0.143261## V12 ? ? ? ? ?. ? ? ? ?.## V13 ? ? ? ? ?. ? ? ? ?.## V14 ? ? ? ? -0.91090 -0.911207## V15 ? ? ? ? ?. ? ? ? ?.## V16 ? ? ? ? ?. ? ? ? ?.## V17 ? ? ? ? ?. ? ? ? ?.## V18 ? ? ? ? ?. ? ? ? ?0.009197## V19 ? ? ? ? ?. ? ? ? ?.## V20 ? ? ? ? -0.86099 -0.863117
左列是，exact = TRUE?右列是?FALSE。從上面我們可以看到，0.01不在序列中，因此盡管沒有太大差異，但還是有一些差異。如果沒有特殊要求，則線性插補就足夠了。
用戶可以根據擬合的對象進行預測。除中的選項外?coef，主要參數是?newx的新值矩陣?x。type?選項允許用戶選擇預測類型：*“鏈接”給出擬合值
因變量與正態(tài)分布的“鏈接”相同。
“系數”計算值為的系數?
s
例如，
## ? ? ? ? ? ?1## [1,] -0.9803## [2,] ?2.2992## [3,] ?0.6011## [4,] ?2.3573## [5,] ?1.7520
給出在λ=0.05時前5個觀測值的擬合值。如果提供的多個值，?s?則會生成預測矩陣。
用戶可以自定義K折交叉驗證。除所有?glmnet?參數外，?cv.glmnet?還有特殊的參數，包括?nfolds?（次數），??foldid?（用戶提供的次數），??type.measure（用于交叉驗證的損失）：*“ deviance”或“ mse”
“ mae”使用平均絕對誤差
舉個例子，
cvfit = cv.glmnet(x, y, type.measure = "mse", nfolds = 20)
根據均方誤差標準進行20折交叉驗證。
并行計算也受?cv.glmnet。為我們在這里給出一個簡單的比較示例。
system.time(cv.glmnet(X, Y))
## ? ?user ?system elapsed## ? 3.591 ? 0.103 ? 3.724
system.time(cv.glmnet(X, Y, parallel = TRUE))## ? ?user ?system elapsed
## ? 4.318 ? 0.391 ? 2.700
從上面的建議可以看出，并行計算可以大大加快計算過程。
“ lambda.min”：達到最小MSE的λ。
cvfit$lambda.min## [1] 0.08307
## 21 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ?1## (Intercept) ?0.14936## V1 ? ? ? ? ? 1.32975## V2 ? ? ? ? ? .## V3 ? ? ? ? ? 0.69096## V4 ? ? ? ? ? .## V5 ? ? ? ? ?-0.83123## V6 ? ? ? ? ? 0.53670## V7 ? ? ? ? ? 0.02005## V8 ? ? ? ? ? 0.33194## V9 ? ? ? ? ? .## V10 ? ? ? ? ?.## V11 ? ? ? ? ?0.16239## V12 ? ? ? ? ?.## V13 ? ? ? ? ?.## V14 ? ? ? ? -1.07081## V15 ? ? ? ? ?.## V16 ? ? ? ? ?.## V17 ? ? ? ? ?.## V18 ? ? ? ? ?.## V19 ? ? ? ? ?.## V20 ? ? ? ? -1.04341
在這里，我們使用相同的k折，為α選擇一個值。
將它們全部放置在同一繪圖上：
我們看到lasso（alpha=1）在這里表現最好。
系數上下限
假設我們要擬合我們的模型，但將系數限制為大于-0.7且小于0.5。這可以通過upper.limits?和?lower.limits?參數實現?：
通常，我們希望系數為正，因此我們只能lower.limit?將其設置?為0。
懲罰因素
此參數允許用戶將單獨的懲罰因子應用于每個系數。每個參數的默認值為1，但可以指定其他值。特別是，任何penalty.factor?等于零的變量?都不會受到懲罰
?
在許多情況下，某些變量可能是重要，我們希望一直保留它們，這可以通過將相應的懲罰因子設置為0來實現：
我們從標簽中看到懲罰因子為0的三個變量始終保留在模型中，而其他變量遵循典型的正則化路徑并最終縮小為0。
?
自定義圖
有時，尤其是在變量數量很少的情況下，我們想在圖上添加變量標簽。
我們首先生成帶有10個變量的一些數據，然后，我們擬合glmnet模型，并繪制標準圖。
我們希望用變量名標記曲線。在路徑的末尾放置系數的位置。
?
多元正態(tài)
使用family = "mgaussian"?option?獲得多元正態(tài)分布glmnet。
顯然，顧名思義，y不是向量，而是矩陣。結果，每個λ值的系數也是一個矩陣。
在這里，我們解決以下問題：
這里，βj是p×K系數矩陣β的第j行，對于單個預測變量xj，我們用每個系數K向量βj的組套索罰分代替每個單一系數的絕對罰分。
我們使用預先生成的一組數據進行說明。
我們擬合數據，并返回對象“ mfit”。
mfit = glmnet(x, y, family = "mgaussian")
如果為?standardize.response = TRUE，則將因變量標準化。
為了可視化系數，我們使用?plot?函數。
注意我們設置了?type.coef = "2norm"。在此設置下，每個變量繪制一條曲線，其值等于?2范數。默認設置為?type.coef = "coef"，其中為每個因變量創(chuàng)建一個系數圖。
通過使用該函數coef?，我們可以提取要求的λ值的系數，?并通過進行預測?。
## , , 1#### ? ? ? ? ? y1 ? ? ?y2 ? ? ?y3 ? ?y4## [1,] -4.7106 -1.1635 ?0.6028 3.741## [2,] ?4.1302 -3.0508 -1.2123 4.970## [3,] ?3.1595 -0.5760 ?0.2608 2.054## [4,] ?0.6459 ?2.1206 -0.2252 3.146## [5,] -1.1792 ?0.1056 -7.3353 3.248#### , , 2#### ? ? ? ? ? y1 ? ? ?y2 ? ? ?y3 ? ?y4## [1,] -4.6415 -1.2290 ?0.6118 3.780## [2,] ?4.4713 -3.2530 -1.2573 5.266## [3,] ?3.4735 -0.6929 ?0.4684 2.056## [4,] ?0.7353 ?2.2965 -0.2190 2.989## [5,] -1.2760 ?0.2893 -7.8259 3.205
預測結果保存在三維數組中，其中前兩個維是每個因變量的預測矩陣，第三個維表示因變量。
我們還可以進行k折交叉驗證。
我們繪制結果?cv.glmnet?對象“ cvmfit”。
顯示選定的λ最佳值
cvmfit$lambda.min## [1] 0.04732cvmfit$lambda.1se## [1] 0.1317
邏輯回歸
當因變量是分類的時，邏輯回歸是另一個廣泛使用的模型。如果有兩個可能的結果，則使用二項式分布，否則使用多項式。
二項式模型
對于二項式模型，假設因變量的取值為G = {1,2} 。表示yi = I（gi = 1）。我們建模
可以用以下形式寫
?
懲罰邏輯回歸的目標函數使用負二項式對數似然
我們的算法使用對數似然的二次逼近，然后對所得的懲罰加權最小二乘問題進行下降。這些構成了內部和外部循環(huán)。
出于說明目的，我們?從數據文件加載預生成的輸入矩陣?x?和因變量?y。
對于二項式邏輯回歸，因變量y可以是兩個級別的因子，也可以是計數或比例的兩列矩陣。
glmnet?二項式回歸的其他可選參數與正態(tài)分布的參數?幾乎相同。不要忘記將family?選項設置?為“ binomial”。
fit = glmnet(x, y, family = "binomial")
像以前一樣，我們可以輸出和繪制擬合的對象，提取特定λ處的系數，并進行預測。
邏輯回歸略有不同，主要體現在選擇上?type?！版溄印焙汀耙蜃兞俊辈坏葍r，“類”僅可用于邏輯回歸?？傊?，*“鏈接”給出了線性預測變量
“因變量”給出合適的概率
“類別”產生對應于最大概率的類別標簽。
“系數”計算值為的系數?
s
在下面的示例中，我們在λ=0.05,0.01的情況下對類別標簽進行了預測。
## ? ? ?1 ? 2## [1,] "0" "0"## [2,] "1" "1"## [3,] "1" "1"## [4,] "0" "0"## [5,] "1" "1"
對于邏輯回歸，type.measure：
“偏差”使用實際偏差。
“ mae”使用平均絕對誤差。
“class”給出錯誤分類錯誤。
“ auc”（僅適用于兩類邏輯回歸）給出了ROC曲線下的面積。
例如，
它使用分類誤差作為10倍交叉驗證的標準。
我們繪制對象并顯示λ的最佳值。
cvfit$lambda.min## [1] 0.01476cvfit$lambda.1se## [1] 0.02579
coef?并且?predict?類似于正態(tài)分布案例，因此我們省略了細節(jié)。我們通過一些例子進行回顧。
## 31 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ?1## (Intercept) ?0.24371## V1 ? ? ? ? ? 0.06897## V2 ? ? ? ? ? 0.66252## V3 ? ? ? ? ?-0.54275## V4 ? ? ? ? ?-1.13693## V5 ? ? ? ? ?-0.19143## V6 ? ? ? ? ?-0.95852## V7 ? ? ? ? ? .## V8 ? ? ? ? ?-0.56529## V9 ? ? ? ? ? 0.77454## V10 ? ? ? ? -1.45079## V11 ? ? ? ? -0.04363## V12 ? ? ? ? -0.06894## V13 ? ? ? ? ?.## V14 ? ? ? ? ?.## V15 ? ? ? ? ?.## V16 ? ? ? ? ?0.36685## V17 ? ? ? ? ?.## V18 ? ? ? ? -0.04014## V19 ? ? ? ? ?.## V20 ? ? ? ? ?.## V21 ? ? ? ? ?.## V22 ? ? ? ? ?0.20882## V23 ? ? ? ? ?0.34014## V24 ? ? ? ? ?.## V25 ? ? ? ? ?0.66310## V26 ? ? ? ? -0.33696## V27 ? ? ? ? -0.10570## V28 ? ? ? ? ?0.24318## V29 ? ? ? ? -0.22445## V30 ? ? ? ? ?0.11091
如前所述，此處返回的結果僅針對因子因變量的第二類。
## ? ? ? 1## ?[1,] "0"## ?[2,] "1"## ?[3,] "1"## ?[4,] "0"## ?[5,] "1"## ?[6,] "0"## ?[7,] "0"## ?[8,] "0"## ?[9,] "1"## [10,] "1"
多項式模型
對于多項式模型，假設因變量變量的K級別為G = {1,2，…，K}。在這里我們建模
設Y為N×K指標因變量矩陣，元素yi?= I（gi =?）。然后彈性網懲罰的負對數似然函數變?yōu)?/p>
β是系數的p×K矩陣。βk指第k列（對于結果類別k），βj指第j行（變量j的K個系數的向量）。最后一個懲罰項是||βj|| q ，我們對q有兩個選擇：q∈{1,2}。當q = 1時，這是每個參數的套索懲罰。當q = 2時，這是對特定變量的所有K個系數的分組套索懲罰，這使它們在一起全為零或非零。
對于多項式情況，用法類似于邏輯回歸，我們加載一組生成的數據。
glmnet?除少數情況外，多項式邏輯回歸中的可選參數?與二項式回歸基本相似。
多項式回歸的一個特殊選項是?type.multinomial，如果允許，則允許使用分組的套索罰分?type.multinomial = "grouped"。這將確保變量的多項式系數全部一起輸入或輸出，就像多元因變量一樣。
我們繪制結果。
我們還可以進行交叉驗證并繪制返回的對象。
?
預測最佳選擇的λ：
## ? ? ? 1## ?[1,] "3"## ?[2,] "2"## ?[3,] "2"## ?[4,] "1"## ?[5,] "1"## ?[6,] "3"## ?[7,] "3"## ?[8,] "1"## ?[9,] "1"## [10,] "2"
?
泊松模型
Poisson回歸用于在假設Poisson誤差的情況下對計數數據進行建模，或者在均值和方差成比例的情況下使用非負數據進行建模。泊松也是指數分布族的成員。我們通常以對數建模：
。
給定觀測值
的對數似然
?
和以前一樣，我們優(yōu)化了懲罰對數：
Glmnet使用外部牛頓循環(huán)和內部加權最小二乘循環(huán)（如邏輯回歸）來優(yōu)化此標準。
首先，我們加載一組泊松數據。
?
再次，繪制系數。
像以前一樣，我們可以?分別使用coef?和?提取系數并在特定的λ處進行預測?predict。
例如，我們可以
## 21 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ?1## (Intercept) ?0.61123## V1 ? ? ? ? ? 0.45820## V2 ? ? ? ? ?-0.77061## V3 ? ? ? ? ? 1.34015## V4 ? ? ? ? ? 0.04350## V5 ? ? ? ? ?-0.20326## V6 ? ? ? ? ? .## V7 ? ? ? ? ? .## V8 ? ? ? ? ? .## V9 ? ? ? ? ? .## V10 ? ? ? ? ?.## V11 ? ? ? ? ?.## V12 ? ? ? ? ?0.01816## V13 ? ? ? ? ?.## V14 ? ? ? ? ?.## V15 ? ? ? ? ?.## V16 ? ? ? ? ?.## V17 ? ? ? ? ?.## V18 ? ? ? ? ?.## V19 ? ? ? ? ?.## V20 ? ? ? ? ?.
## ? ? ? ? ? ?1 ? ? ? 2
## [1,] ?2.4944 ?4.4263
## [2,] 10.3513 11.0586
## [3,] ?0.1180 ?0.1782
## [4,] ?0.9713 ?1.6829
## [5,] ?1.1133 ?1.9935
我們還可以使用交叉驗證來找到最佳的λ，從而進行推斷。
選項幾乎與正態(tài)族相同，不同之處在于?type.measure*，“ mse”代表均方誤差*，“ mae”代表均值絕對誤差。
我們可以繪制?cv.glmnet?對象。
我們還可以顯示最佳的λ和相應的系數。
## 21 x 2 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ? ? ? ? ? 1 ? ? ? ?2## (Intercept) ?0.031263 ?0.18570## V1 ? ? ? ? ? 0.619053 ?0.57537## V2 ? ? ? ? ?-0.984550 -0.93212## V3 ? ? ? ? ? 1.525234 ?1.47057## V4 ? ? ? ? ? 0.231591 ?0.19692## V5 ? ? ? ? ?-0.336659 -0.30469## V6 ? ? ? ? ? 0.001026 ?.## V7 ? ? ? ? ?-0.012830 ?.## V8 ? ? ? ? ? . ? ? ? ? .## V9 ? ? ? ? ? . ? ? ? ? .## V10 ? ? ? ? ?0.015983 ?.## V11 ? ? ? ? ?. ? ? ? ? .## V12 ? ? ? ? ?0.030867 ?0.02585## V13 ? ? ? ? -0.027971 ?.## V14 ? ? ? ? ?0.032750 ?.## V15 ? ? ? ? -0.005933 ?.## V16 ? ? ? ? ?0.017506 ?.## V17 ? ? ? ? ?. ? ? ? ? .## V18 ? ? ? ? ?0.004026 ?.## V19 ? ? ? ? -0.033579 ?.## V20 ? ? ? ? ?0.012049 ?0.00993
?
Cox模型
Cox比例風險模型通常用于研究預測變量與生存時間之間的關系。
Cox比例風險回歸模型，它不是直接考察
?與X的關系，而是用
?作為因變量，模型的基本形式為：
式中，
?為自變量的偏回歸系數，它是須從樣本數據作出估計的參數；
?是當X向量為0時，
?的基準危險率，它是有待于從樣本數據作出估計的量。簡稱為Cox回歸模型。
由于Cox回歸模型對
?未作任何假定，因此Cox回歸模型在處理問題時具有較大的靈活性；另一方面，在許多情況下，我們只需估計出參數
?(如因素分析等)，即使在
?未知的情況下，仍可估計出參數
?。這就是說，Cox回歸模型由于含有
?，因此它不是完全的參數模型，但仍可根據公式(1)作出參數
?的估計，故Cox回歸模型屬于半參數模型。
公式可以轉化為：
?
我們使用一組預先生成的樣本數據。用戶可以加載自己的數據并遵循類似的過程。在這種情況下，x必須是協變量值的n×p矩陣-每行對應一個患者，每列對應一個協變量。y是一個n×2矩陣。
## ? ? ? ? time status## [1,] 1.76878 ? ? ?1## [2,] 0.54528 ? ? ?1## [3,] 0.04486 ? ? ?0## [4,] 0.85032 ? ? ?0## [5,] 0.61488 ? ? ?1
Surv?包中的?函數?survival?可以創(chuàng)建這樣的矩陣。
我們計算默認設置下的求解路徑。
繪制系數。
提取特定值λ處的系數。
## 30 x 1 sparse Matrix of class "dgCMatrix"## ? ? ? ? ? ?1## V1 ? 0.37694## V2 ?-0.09548## V3 ?-0.13596## V4 ? 0.09814## V5 ?-0.11438## V6 ?-0.38899## V7 ? 0.24291## V8 ? 0.03648## V9 ? 0.34740## V10 ?0.03865## V11 ?.## V12 ?.## V13 ?.## V14 ?.## V15 ?.## V16 ?.## V17 ?.## V18 ?.## V19 ?.## V20 ?.## V21 ?.## V22 ?.## V23 ?.## V24 ?.## V25 ?.## V26 ?.## V27 ?.## V28 ?.## V29 ?.## V30 ?.
函數?cv.glmnet?可用于計算Cox模型的k折交叉驗證。
擬合后，我們可以查看最佳λ值和交叉驗證的誤差圖，幫助評估我們的模型。
如前所述，圖中的左垂直線向我們顯示了CV誤差曲線達到最小值的位置。右邊的垂直線向我們展示了正則化的模型，其CV誤差在最小值的1個標準偏差之內。我們還提取了最優(yōu)λ。
cvfit$lambda.min## [1] 0.01594cvfit$lambda.1se## [1] 0.04869
我們可以檢查模型中的協變量并查看其系數。
index.min
## ?[1] ?0.491297 -0.174601 -0.218649 ?0.175112 -0.186673 -0.490250 ?0.335197## ?[8] ?0.091587 ?0.450169 ?0.115922 ?0.017595 -0.018365 -0.002806 -0.001423## [15] -0.023429 ?0.001688 -0.008236
coef.min## 30 x 1 sparse Matrix of class "dgCMatrix"
## ? ? ? ? ? ? 1
## V1 ? 0.491297
## V2 ?-0.174601
## V3 ?-0.218649
## V4 ? 0.175112
## V5 ?-0.186673
## V6 ?-0.490250
## V7 ? 0.335197
## V8 ? 0.091587
## V9 ? 0.450169
## V10 ?0.115922
## V11 ?.
## V12 ?.
## V13 ?0.017595
## V14 ?.
## V15 ?.
## V16 ?.
## V17 -0.018365
## V18 ?.
## V19 ?.
## V20 ?.
## V21 -0.002806
## V22 -0.001423
## V23 ?.
## V24 ?.
## V25 -0.023429
## V26 ?.
## V27 ?0.001688
## V28 ?.
## V29 ?.
## V30 -0.008236
?
稀疏矩陣
我們的程序包支持稀疏的輸入矩陣，該矩陣可以高效地存儲和操作大型矩陣，但只有少數幾個非零條目。
我們加載一組預先創(chuàng)建的樣本數據。
加載100 * 20的稀疏矩陣和?y因向量。
## [1] "dgCMatrix"## attr(,"package")## [1] "Matrix"
我們可以像以前一樣擬合模型。
fit = glmnet(x, y)
進行交叉驗證并繪制結果對象。
預測新輸入矩陣?。例如，
## ? ? ? ? ? ?1## [1,] ?0.3826## [2,] -0.2172## [3,] -1.6622## [4,] -0.4175## [5,] -1.3941
?
參考文獻
Jerome Friedman, Trevor Hastie and Rob Tibshirani. (2008).
Regularization Paths for Generalized Linear Models via Coordinate Descent
參考文獻
1.matlab偏最小二乘回歸(PLSR)和主成分回歸(PCR)
2.R語言高維數據的主成分pca、 t-SNE算法降維與可視化分析
3.主成分分析(PCA)基本原理及分析實例
4.基于R語言實現LASSO回歸分析
5.使用LASSO回歸預測股票收益數據分析
6.r語言中對lasso回歸，ridge嶺回歸和elastic-net模型
7.r語言中的偏最小二乘回歸pls-da數據分析
8.r語言中的偏最小二乘pls回歸算法
9.R語言線性判別分析（LDA），二次判別分析（QDA）和正則判別分析（RDA）
