原文鏈接：http://tecdat.cn/?p=3795

Glmnet是一個通過懲罰最大似然關(guān)系擬合廣義線性模型的軟件包。正則化路徑是針對正則化參數(shù)λ的值網(wǎng)格處的lasso或Elastic Net（彈性網(wǎng)絡(luò)）懲罰值計算的?（?點擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)********?）。

該算法非常快，并且可以利用輸入矩陣中的稀疏性?x。它適合線性，邏輯和多項式，泊松和Cox回歸模型。可以從擬合模型中做出各種預(yù)測。

它也可以擬合多元線性回歸。

glmnet?解決以下問題

在覆蓋整個范圍的λ值網(wǎng)格上。這里l（y，η）是觀察i的負(fù)對數(shù)似然貢獻；例如對于高斯分布是

。?_彈性網(wǎng)絡(luò)_懲罰由α控制，LASSO（α= 1，默認(rèn)），Ridge（α= 0）。調(diào)整參數(shù)λ控制懲罰的總強度。

眾所周知，嶺懲罰使相關(guān)預(yù)測因子的系數(shù)彼此縮小，而套索傾向于選擇其中一個而丟棄其他預(yù)測因子。_彈性網(wǎng)絡(luò)_則將這兩者混合在一起。

glmnet?算法使用循環(huán)坐標(biāo)下降法，該方法在每個參數(shù)固定不變的情況下連續(xù)優(yōu)化目標(biāo)函數(shù)，并反復(fù)循環(huán)直到收斂，我們的算法可以非常快速地計算求解路徑。

代碼可以處理稀疏的輸入矩陣格式，以及系數(shù)的范圍約束，還包括用于預(yù)測和繪圖的方法，以及執(zhí)行K折交叉驗證的功能。

快速開始

首先，我們加載?glmnet?包：

library(glmnet)

包中使用的默認(rèn)模型是高斯線性模型或“最小二乘”。我們加載一組預(yù)先創(chuàng)建的數(shù)據(jù)以進行說明。用戶可以加載自己的數(shù)據(jù)，也可以使用工作空間中保存的數(shù)據(jù)。

該命令?從此保存的R數(shù)據(jù)中加載輸入矩陣?x?和因向量?y。

我們擬合模型?glmnet。

fit?=?glmnet(x,?y)

可以通過執(zhí)行plot?函數(shù)來可視化系數(shù)?：

plot(fit)

每條曲線對應(yīng)一個變量。它顯示了當(dāng)λ變化時，其系數(shù)相對于整個系數(shù)向量的?1范數(shù)的路徑。上方的軸表示當(dāng)前λ處非零系數(shù)的數(shù)量，這是套索的有效自由度（df）。用戶可能還希望對曲線進行注釋。這可以通過label = TRUE?在plot命令中進行設(shè)置來完成?。

點擊標(biāo)題查閱往期內(nèi)容

R語言自適應(yīng)LASSO 多項式回歸、二元邏輯回歸和嶺回歸應(yīng)用分析

左右滑動查看更多

01

02

03

04

glmnet?如果我們只是輸入對象名稱或使用print?函數(shù)，則會顯示每個步驟的路徑?摘要?：

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

它從左到右顯示了非零系數(shù)的數(shù)量（Df），解釋的（零）偏差百分比（%dev）和λ（Lambda）的值。

我們可以在序列范圍內(nèi)獲得一個或多個λ處的實際系數(shù)：

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

還可以使用新的輸入數(shù)據(jù)在特定的λ處進行預(yù)測：

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

該函數(shù)?glmnet?返回一系列模型供用戶選擇。交叉驗證可能是該任務(wù)最簡單，使用最廣泛的方法。

cv.glmnet?是交叉驗證的主要函數(shù)。

cv.glmnet?返回一個?cv.glmnet?對象，此處為“ cvfit”，其中包含交叉驗證擬合的所有成分的列表。

我們可以繪制對象。

它包括交叉驗證曲線（紅色虛線）和沿λ序列的上下標(biāo)準(zhǔn)偏差曲線（誤差線）。垂直虛線表示兩個選定的λ。

我們可以查看所選的λ和相應(yīng)的系數(shù)。例如，

cvfit$lambda.min##?[1]?0.08307

lambda.min?是給出最小平均交叉驗證誤差的λ值。保存的另一個λ是?lambda.1se，它給出了的模型，使得誤差在最小值的一個標(biāo)準(zhǔn)誤差以內(nèi)。我們只需要更換?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

注意，系數(shù)以稀疏矩陣格式表示。原因是沿著正則化路徑的解通常是稀疏的，因此使用稀疏格式在時間和空間上更為有效。

可以根據(jù)擬合的cv.glmnet?對象進行預(yù)測?。讓我們看一個示例。

##????????????1##?[1,]?-1.3647##?[2,]??2.5686##?[3,]??0.5706##?[4,]??1.9682##?[5,]??1.4964

newx?與新的輸入矩陣?s相同，如前所述，是預(yù)測的λ值。

線性回歸

這里的線性回歸是指兩個模型系列。一個是?gaussian正態(tài)_分布_，另一個是?mgaussian多元正態(tài)_分布_。

正態(tài)_分布_

假設(shè)我們有觀測值xi∈Rp并且yi∈R，i = 1，...，N。目標(biāo)函數(shù)是

其中λ≥0是復(fù)雜度參數(shù)，0≤α≤1在嶺回歸（α=0）和套索LASSO（α=1）之間。

應(yīng)用坐標(biāo)下降法解決該問題。具體地說，通過計算βj=β?j處的梯度和簡單的演算，更新為

其中

。

當(dāng)x?變量標(biāo)準(zhǔn)化為具有單位方差（默認(rèn)值）時，以上公式適用?。

glmnet?提供各種選項供用戶自定義。我們在這里介紹一些常用的選項，它們可以在glmnet?函數(shù)中指定?。

• alpha?表示彈性網(wǎng)混合參數(shù)α，范圍α∈[0,1]。α=1是套索（默認(rèn)），α=0是Ridge。

• weights?用于觀察權(quán)重。每個觀察值的默認(rèn)值為1。

• nlambda?是序列中λ值的數(shù)量。默認(rèn)值為100。

• lambda?可以提供，但通常不提供，程序會構(gòu)建一個序列。自動生成時，λ序列由lambda.max?和?確定?lambda.min.ratio。

• standardize?是x?在擬合模型序列之前進行變量標(biāo)準(zhǔn)化的邏輯標(biāo)志?。

例如，我們設(shè)置α=0.2，并對后半部分的觀測值賦予兩倍的權(quán)重。為了避免在此處顯示太長時間，我們將其設(shè)置?nlambda?為20。但是，實際上，建議將λ的數(shù)量設(shè)置為100（默認(rèn)值）或更多。

然后我們可以輸出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

這將顯示生成對象的調(diào)用?fit?以及帶有列Df?（非零系數(shù)的數(shù)量），??%dev?（解釋的偏差百分比）和Lambda?（對應(yīng)的λ值）?的三列矩陣?。

我們可以繪制擬合的對象。

讓我們針對log-lambda值標(biāo)記每個曲線來繪制“擬合”。

這是訓(xùn)練數(shù)據(jù)中的偏差百分比。我們在這里看到的是，在路徑末端時，該值變化不大，但是系數(shù)有點“膨脹”。這使我們可以將注意力集中在重要的擬合部分上。

我們可以提取系數(shù)并在某些特定值的情況下進行預(yù)測。兩種常用的選項是：

• s?指定進行提取的λ值。

• exact?指示是否需要系數(shù)的精確值。

一個簡單的例子是：

##?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不在序列中，因此盡管沒有太大差異，但還是有一些差異。如果沒有特殊要求，則線性插補就足夠了。

用戶可以根據(jù)擬合的對象進行預(yù)測。除中的選項外?coef，主要參數(shù)是?newx的新值矩陣?x。type?選項允許用戶選擇預(yù)測類型：*“鏈接”給出擬合值

• 因變量與正態(tài)分布的“鏈接”相同。

• “系數(shù)”計算值為的系數(shù)?s

例如，

##????????????1##?[1,]?-0.9803##?[2,]??2.2992##?[3,]??0.6011##?[4,]??2.3573##?[5,]??1.7520

給出在λ=0.05時前5個觀測值的擬合值。如果提供的多個值，?s?則會生成預(yù)測矩陣。

用戶可以自定義K折交叉驗證。除所有?glmnet?參數(shù)外，?cv.glmnet?還有特殊的參數(shù)，包括?nfolds?（次數(shù)），??foldid?（用戶提供的次數(shù)），??type.measure（用于交叉驗證的損失）：*“ deviance”或“ mse”

• “ mae”使用平均絕對誤差

舉個例子，

cvfit?=?cv.glmnet(x,?y,?type.measure?=?"mse",?nfolds?=?20)

根據(jù)均方誤差標(biāo)準(zhǔn)進行20折交叉驗證。

并行計算也受?cv.glmnet。為我們在這里給出一個簡單的比較示例。

system.time(cv.glmnet(X,?Y))##????user??system?elapsed?##???3.591???0.103???3.724system.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折，為α選擇一個值。

將它們?nèi)糠胖迷谕焕L圖上：

我們看到lasso（alpha=1）在這里表現(xiàn)最好。

系數(shù)上下限

假設(shè)我們要擬合我們的模型，但將系數(shù)限制為大于-0.7且小于0.5。這可以通過upper.limits?和?lower.limits?參數(shù)實現(xiàn)?：

通常，我們希望系數(shù)為正，因此我們只能lower.limit?將其設(shè)置?為0。

懲罰因素

此參數(shù)允許用戶將單獨的懲罰因子應(yīng)用于每個系數(shù)。每個參數(shù)的默認(rèn)值為1，但可以指定其他值。特別是，任何penalty.factor?等于零的變量?都不會受到懲罰

在許多情況下，某些變量可能是重要，我們希望一直保留它們，這可以通過將相應(yīng)的懲罰因子設(shè)置為0來實現(xiàn)：

我們從標(biāo)簽中看到懲罰因子為0的三個變量始終保留在模型中，而其他變量遵循典型的正則化路徑并最終縮小為0。

自定義圖

有時，尤其是在變量數(shù)量很少的情況下，我們想在圖上添加變量標(biāo)簽。

我們首先生成帶有10個變量的一些數(shù)據(jù)，然后，我們擬合glmnet模型，并繪制標(biāo)準(zhǔn)圖。

我們希望用變量名標(biāo)記曲線。在路徑的末尾放置系數(shù)的位置。

多元正態(tài)

使用family = "mgaussian"?option?獲得多元正態(tài)分布glmnet。

顯然，顧名思義，y不是向量，而是矩陣。結(jié)果，每個λ值的系數(shù)也是一個矩陣。

在這里，我們解決以下問題：

這里，βj是p×K系數(shù)矩陣β的第j行，對于單個預(yù)測變量xj，我們用每個系數(shù)K向量βj的組套索罰分代替每個單一系數(shù)的絕對罰分。

我們使用預(yù)先生成的一組數(shù)據(jù)進行說明。

我們擬合數(shù)據(jù)，并返回對象“ mfit”。

mfit?=?glmnet(x,?y,?family?=?"mgaussian")

如果為?standardize.response = TRUE，則將因變量標(biāo)準(zhǔn)化。

為了可視化系數(shù)，我們使用?plot?函數(shù)。

注意我們設(shè)置了?type.coef = "2norm"。在此設(shè)置下，每個變量繪制一條曲線，其值等于?2范數(shù)。默認(rèn)設(shè)置為?type.coef = "coef"，其中為每個因變量創(chuàng)建一個系數(shù)圖。

通過使用該函數(shù)coef?，我們可以提取要求的λ值的系數(shù)，?并通過進行預(yù)測?。

##?,?,?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

預(yù)測結(jié)果保存在三維數(shù)組中，其中前兩個維是每個因變量的預(yù)測矩陣，第三個維表示因變量。

我們還可以進行k折交叉驗證。

我們繪制結(jié)果?cv.glmnet?對象“ cvmfit”。

顯示選定的λ最佳值

cvmfit$lambda.min##?[1]?0.04732cvmfit$lambda.1se##?[1]?0.1317

邏輯回歸

當(dāng)因變量是分類的時，邏輯回歸是另一個廣泛使用的模型。如果有兩個可能的結(jié)果，則使用二項式分布，否則使用多項式。

二項式模型

對于二項式模型，假設(shè)因變量的取值為G = {1,2} 。表示yi = I（gi = 1）。我們建模

可以用以下形式寫

懲罰邏輯回歸的目標(biāo)函數(shù)使用負(fù)二項式對數(shù)似然

我們的算法使用對數(shù)似然的二次逼近，然后對所得的懲罰加權(quán)最小二乘問題進行下降。這些構(gòu)成了內(nèi)部和外部循環(huán)。

出于說明目的，我們?從數(shù)據(jù)文件加載預(yù)生成的輸入矩陣?x?和因變量?y。

對于二項式邏輯回歸，因變量y可以是兩個級別的因子，也可以是計數(shù)或比例的兩列矩陣。

glmnet?二項式回歸的其他可選參數(shù)與正態(tài)分布的參數(shù)?幾乎相同。不要忘記將family?選項設(shè)置?為“ binomial”。

fit?=?glmnet(x,?y,?family?=?"binomial")

像以前一樣，我們可以輸出和繪制擬合的對象，提取特定λ處的系數(shù)，并進行預(yù)測。

邏輯回歸略有不同，主要體現(xiàn)在選擇上?type?！版溄印焙汀耙蜃兞俊辈坏葍r，“類”僅可用于邏輯回歸。總之，*“鏈接”給出了線性預(yù)測變量

• “因變量”給出合適的概率

• “類別”產(chǎn)生對應(yīng)于最大概率的類別標(biāo)簽。

• “系數(shù)”計算值為的系數(shù)?s

在下面的示例中，我們在λ=0.05,0.01的情況下對類別標(biāo)簽進行了預(yù)測。

##??????1???2??##?[1,]?"0"?"0"##?[2,]?"1"?"1"##?[3,]?"1"?"1"##?[4,]?"0"?"0"##?[5,]?"1"?"1"

對于邏輯回歸，type.measure：

• “偏差”使用實際偏差。

• “ mae”使用平均絕對誤差。

• “class”給出錯誤分類錯誤。

• “ auc”（僅適用于兩類邏輯回歸）給出了ROC曲線下的面積。

例如，

它使用分類誤差作為10倍交叉驗證的標(biāo)準(zhǔn)。

我們繪制對象并顯示λ的最佳值。

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

如前所述，此處返回的結(jié)果僅針對因子因變量的第二類。

##???????1??##??[1,]?"0"##??[2,]?"1"##??[3,]?"1"##??[4,]?"0"##??[5,]?"1"##??[6,]?"0"##??[7,]?"0"##??[8,]?"0"##??[9,]?"1"##?[10,]?"1"

多項式模型

對于多項式模型，假設(shè)因變量變量的K級別為G = {1,2，…，K}。在這里我們建模

設(shè)Y為N×K指標(biāo)因變量矩陣，元素yi?= I（gi =?）。然后彈性網(wǎng)懲罰的負(fù)對數(shù)似然函數(shù)變?yōu)?/p>

β是系數(shù)的p×K矩陣。βk指第k列（對于結(jié)果類別k），βj指第j行（變量j的K個系數(shù)的向量）。最后一個懲罰項是||βj|| q ，我們對q有兩個選擇：q∈{1,2}。當(dāng)q = 1時，這是每個參數(shù)的套索懲罰。當(dāng)q = 2時，這是對特定變量的所有K個系數(shù)的分組套索懲罰，這使它們在一起全為零或非零。

對于多項式情況，用法類似于邏輯回歸，我們加載一組生成的數(shù)據(jù)。

glmnet?除少數(shù)情況外，多項式邏輯回歸中的可選參數(shù)?與二項式回歸基本相似。

多項式回歸的一個特殊選項是?type.multinomial，如果允許，則允許使用分組的套索罰分?type.multinomial = "grouped"。這將確保變量的多項式系數(shù)全部一起輸入或輸出，就像多元因變量一樣。

我們繪制結(jié)果。

我們還可以進行交叉驗證并繪制返回的對象。

預(yù)測最佳選擇的λ：

##???????1??##??[1,]?"3"##??[2,]?"2"##??[3,]?"2"##??[4,]?"1"##??[5,]?"1"##??[6,]?"3"##??[7,]?"3"##??[8,]?"1"##??[9,]?"1"##?[10,]?"2"

泊松模型

Poisson回歸用于在假設(shè)Poisson誤差的情況下對計數(shù)數(shù)據(jù)進行建模，或者在均值和方差成比例的情況下使用非負(fù)數(shù)據(jù)進行建模。泊松也是指數(shù)分布族的成員。我們通常以對數(shù)建模：。
給定觀測值的對數(shù)似然

和以前一樣，我們優(yōu)化了懲罰對數(shù)：

Glmnet使用外部牛頓循環(huán)和內(nèi)部加權(quán)最小二乘循環(huán)（如邏輯回歸）來優(yōu)化此標(biāo)準(zhǔn)。

首先，我們加載一組泊松數(shù)據(jù)。

再次，繪制系數(shù)。

像以前一樣，我們可以?分別使用coef?和?提取系數(shù)并在特定的λ處進行預(yù)測?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?對象。

我們還可以顯示最佳的λ和相應(yīng)的系數(shù)。

##?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比例風(fēng)險模型通常用于研究預(yù)測變量與生存時間之間的關(guān)系。

Cox比例風(fēng)險回歸模型，它不是直接考察?與X的關(guān)系，而是用?作為因變量，模型的基本形式為：

式中，?為自變量的偏回歸系數(shù)，它是須從樣本數(shù)據(jù)作出估計的參數(shù)；?是當(dāng)X向量為0時，?的基準(zhǔn)危險率，它是有待于從樣本數(shù)據(jù)作出估計的量。簡稱為Cox回歸模型。

由于Cox回歸模型對?未作任何假定，因此Cox回歸模型在處理問題時具有較大的靈活性；另一方面，在許多情況下，我們只需估計出參數(shù)?(如因素分析等)，即使在?未知的情況下，仍可估計出參數(shù)?。這就是說，Cox回歸模型由于含有?，因此它不是完全的參數(shù)模型，但仍可根據(jù)公式(1)作出參數(shù)?的估計，故Cox回歸模型屬于半?yún)?shù)模型。

公式可以轉(zhuǎn)化為：

我們使用一組預(yù)先生成的樣本數(shù)據(jù)。用戶可以加載自己的數(shù)據(jù)并遵循類似的過程。在這種情況下，x必須是協(xié)變量值的n×p矩陣-每行對應(yīng)一個患者，每列對應(yīng)一個協(xié)變量。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?包中的?函數(shù)?survival?可以創(chuàng)建這樣的矩陣。

我們計算默認(rèn)設(shè)置下的求解路徑。

繪制系數(shù)。

提取特定值λ處的系數(shù)。

##?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??.

函數(shù)?cv.glmnet?可用于計算Cox模型的k折交叉驗證。

擬合后，我們可以查看最佳λ值和交叉驗證的誤差圖，幫助評估我們的模型。

如前所述，圖中的左垂直線向我們顯示了CV誤差曲線達到最小值的位置。右邊的垂直線向我們展示了正則化的模型，其CV誤差在最小值的1個標(biāo)準(zhǔn)偏差之內(nèi)。我們還提取了最優(yōu)λ。

cvfit$lambda.min##?[1]?0.01594cvfit$lambda.1se##?[1]?0.04869

我們可以檢查模型中的協(xié)變量并查看其系數(shù)。

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.008236coef.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

稀疏矩陣

我們的程序包支持稀疏的輸入矩陣，該矩陣可以高效地存儲和操作大型矩陣，但只有少數(shù)幾個非零條目。

我們加載一組預(yù)先創(chuàng)建的樣本數(shù)據(jù)。

加載100 * 20的稀疏矩陣和?y因向量。

##?[1]?"dgCMatrix"##?attr(,"package")##?[1]?"Matrix"

我們可以像以前一樣擬合模型。

fit?=?glmnet(x,?y)

進行交叉驗證并繪制結(jié)果對象。

預(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

欲獲取全文文件，請點擊左下角“閱讀原文”。

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《r語言中對LASSO回歸，Ridge嶺回歸和彈性網(wǎng)絡(luò)Elastic Net模型實現(xiàn)》。

點擊標(biāo)題查閱往期內(nèi)容

基于R語言實現(xiàn)LASSO回歸分析
R語言Lasso回歸模型變量選擇和糖尿病發(fā)展預(yù)測模型
【視頻】Lasso回歸、嶺回歸正則化回歸數(shù)學(xué)原理及R軟件實例
群組變量選擇、組懲罰group lasso套索模型預(yù)測新生兒出生體重風(fēng)險因素數(shù)據(jù)和交叉驗證、可視化
【視頻】Lasso回歸、嶺回歸等正則化回歸數(shù)學(xué)原理及R語言實例R語言Lasso回歸模型變量選擇和糖尿病發(fā)展預(yù)測模型
用LASSO，adaptive LASSO預(yù)測通貨膨脹時間序列MATLAB用Lasso回歸擬合高維數(shù)據(jù)和交叉驗證
群組變量選擇、組懲罰group lasso套索模型預(yù)測新生兒出生體重風(fēng)險因素數(shù)據(jù)和交叉驗證、可視化
高維數(shù)據(jù)懲罰回歸方法：主成分回歸PCR、嶺回歸、lasso、彈性網(wǎng)絡(luò)elastic net分析基因數(shù)據(jù)
Python高維變量選擇:SCAD平滑剪切絕對偏差懲罰、Lasso懲罰函數(shù)比較
R使用LASSO回歸預(yù)測股票收益
廣義線性模型glm泊松回歸的lasso、彈性網(wǎng)絡(luò)分類預(yù)測學(xué)生考試成績數(shù)據(jù)和交叉驗證
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)
R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
r語言中對LASSO回歸，Ridge嶺回歸和彈性網(wǎng)絡(luò)Elastic Net模型實現(xiàn)
R語言高維數(shù)據(jù)懲罰回歸方法：主成分回歸PCR、嶺回歸、lasso、彈性網(wǎng)絡(luò)elastic net分析基因數(shù)據(jù)（含練習(xí)題）
廣義線性模型glm泊松回歸的lasso、彈性網(wǎng)絡(luò)分類預(yù)測學(xué)生考試成績數(shù)據(jù)和交叉驗證
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)
R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
r語言中對LASSO回歸，Ridge嶺回歸和彈性網(wǎng)絡(luò)Elastic Net模型實現(xiàn)
R語言高維數(shù)據(jù)懲罰回歸方法：主成分回歸PCR、嶺回歸、lasso、彈性網(wǎng)絡(luò)elastic net分析基因數(shù)據(jù)（含練習(xí)題）
Python中LARS和Lasso回歸之最小角算法Lars分析波士頓住房數(shù)據(jù)實例
R語言自適應(yīng)LASSO 多項式回歸、二元邏輯回歸和嶺回歸應(yīng)用分析
R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例
Python中的Lasso回歸之最小角算法LARS
r語言中對LASSO回歸，Ridge嶺回歸和彈性網(wǎng)絡(luò)Elastic Net模型實現(xiàn)
R語言實現(xiàn)LASSO回歸——自己編寫LASSO回歸算法
r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現(xiàn)
R使用LASSO回歸預(yù)測股票收益
R語言如何和何時使用glmnet嶺回歸
R語言中的嶺回歸、套索回歸、主成分回歸：線性模型選擇和正則化
Python中的ARIMA模型、SARIMA模型和SARIMAX模型對時間序列預(yù)測
R語言arima，向量自回歸（VAR），周期自回歸(PAR)模型分析溫度時間序列
【視頻】Python和R語言使用指數(shù)加權(quán)平均(EWMA)，ARIMA自回歸移動平均模型預(yù)測時間序列
Python用ARIMA和SARIMA模型預(yù)測銷量時間序列數(shù)據(jù)