懲罰廣義線性模型在遺傳關(guān)聯(lián)研究中的應(yīng)用及R軟件實現(xiàn)*

廣東藥學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系(510310)

張俊國　劉　麗　李麗霞　張　敏　郜艷暉△

?

懲罰廣義線性模型在遺傳關(guān)聯(lián)研究中的應(yīng)用及R軟件實現(xiàn)*

廣東藥學(xué)院公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系(510310)

張俊國劉麗李麗霞張敏郜艷暉△

【提要】目的遺傳關(guān)聯(lián)研究中高維數(shù)據(jù)與日俱增。本文探討基于嶺估計、LASSO和彈性網(wǎng)的廣義線性模型在遺傳關(guān)聯(lián)研究的應(yīng)用及軟件實現(xiàn)，為高維關(guān)聯(lián)分析提供方法學(xué)參考。方法介紹懲罰廣義線性模型原理及軟件實現(xiàn)方法，并采用模擬的連鎖平衡和連鎖不平衡的SNPs關(guān)聯(lián)研究數(shù)據(jù)，以懲罰logistic模型例證R軟件glmnet包對廣義線性模型的擬合。結(jié)果對連鎖平衡和連鎖不平衡SNPs模擬數(shù)據(jù)，LASSO與彈性網(wǎng)均給出稀疏解，較好地選擇有關(guān)聯(lián)SNPs而剔除無關(guān)聯(lián)變量；而嶺估計把所有變量都保留在模型中，模型復(fù)雜度高但相應(yīng)的解釋度未增加。結(jié)論LASSO和彈性網(wǎng)可對高維遺傳關(guān)聯(lián)數(shù)據(jù)進(jìn)行有效降維，篩選變量的同時提供參數(shù)估計，從而降低模型的復(fù)雜度。R軟件的glmnet包靈活擬合各類懲罰廣義線性模型，可在高通量遺傳關(guān)聯(lián)分析中推廣應(yīng)用。

懲罰廣義線性模型遺傳關(guān)聯(lián)研究LASSO彈性網(wǎng)glmnet包

近年來隨著高通量測序技術(shù)在遺傳流行病學(xué)研究中的應(yīng)用，大量的高維遺傳關(guān)聯(lián)數(shù)據(jù)不斷涌現(xiàn)，如動輒成百上千甚至數(shù)以萬計的SNPs(single nucleotide polymorphisms)數(shù)據(jù)呈現(xiàn)高維度的特征，給遺傳關(guān)聯(lián)統(tǒng)計分析帶來極大挑戰(zhàn)[1]。在許多情況下，研究中響應(yīng)變量的分布非正態(tài)，且與高維自變量的關(guān)系也非線性，因此統(tǒng)計分析可在廣義線性模型(generalized linear model，GLM)框架下實現(xiàn)[2]。廣義線性模型是經(jīng)典線性模型的自然推廣，包含有一些極具實用價值的模型，如被學(xué)者廣為熟知的logistic回歸、計數(shù)資料中發(fā)揮重要作用的Poisson回歸等。在高維數(shù)據(jù)中運用廣義線性模型，其首要任務(wù)是降維，去除數(shù)據(jù)中的冗余信息或綜合數(shù)據(jù)中相關(guān)信息以減少自變量的個數(shù)[3]；而如何從這些高維數(shù)據(jù)中攫取有價值的信息，即變量選擇顯得尤為重要。

傳統(tǒng)廣義線性模型通過對數(shù)似然函數(shù)進(jìn)行參數(shù)估計，系數(shù)估計值一般非零，自變量多時模型復(fù)雜度大大提高，而存在多重共線性時參數(shù)估計更易于失效。從模型簡潔度和解釋力的角度看，研究者希望從大量的預(yù)測變量中挑選盡量小的子集，用盡可能簡單的模型獲得最好的解釋效果[4]。針對傳統(tǒng)參數(shù)估計的缺陷，1970年Hoerl和Kennard提出嶺估計，在似然函數(shù)基礎(chǔ)上對系數(shù)平方進(jìn)行懲罰(稱L2懲罰)，可在一定程度上解決共線性問題來提高模型的精確度[5]；但是它將系數(shù)壓縮趨于0卻不為0，因此不能給出簡單且解釋性強(qiáng)的模型。1996年Tibshirani提出稱之為LASSO(least absolute shrinkage and selection operator)的變量選擇方法，用模型系數(shù)的絕對值函數(shù)作為懲罰來壓縮模型(稱L1懲罰)，使絕對值較小的系數(shù)自動壓縮為0，同時實現(xiàn)變量選擇和參數(shù)估計[6]。然而在變量數(shù)p遠(yuǎn)遠(yuǎn)大于觀測值n的情形下，LASSO最多只能選擇n個變量，導(dǎo)致模型過于稀疏。因此Zou等學(xué)者在2005年提出將嶺估計和LASSO凸結(jié)合的懲罰方法，稱為“彈性網(wǎng)(elastic net)”，可有效解決該問題[7]。目前懲罰廣義線性模型可在R軟件中實現(xiàn)，運用Glmnet包可有效、靈活地擬合各類懲罰模型，為高維數(shù)據(jù)提供了有力的分析工具。本文采用模擬的SNPs關(guān)聯(lián)研究數(shù)據(jù)，以懲罰logistic模型例證glmnet包對廣義線性模型的擬合。

模型原理和方法

假定觀測值xi∈Rp，因變量yi∈R,i=1,…,N，并假定因變量取值為1和0，構(gòu)建logistic模型為：

(1)

式(1)通常采用最大似然估計。而懲罰logistic回歸是在對數(shù)似然函數(shù)的基礎(chǔ)上，加入懲罰項[8]，通過最小化懲罰目標(biāo)函數(shù)得到參數(shù)估計值，即：

(2)

研究表明，如果自變量間存在相關(guān)，嶺估計傾向于同時壓縮相關(guān)自變量的系數(shù)，而LASSO傾向于僅選擇相關(guān)自變量中的一個，彈性網(wǎng)則是嶺估計和LASSO的有效結(jié)合。設(shè)置α=0.5將傾向于將相關(guān)變量同時納入或同時排除，當(dāng)α=1-ε(ε為雖小卻大于0的數(shù)值)時，會讓彈性網(wǎng)模型執(zhí)行起來更傾向于LASSO，可以排除許多存在高度相關(guān)的變量。

懲罰logistic回歸的算法多采用坐標(biāo)下降法(coordinate descent)，這是目前關(guān)于LASSO的最快計算方法[9]。算法原理為對每一個參數(shù)在保持其他參數(shù)固定的情況下進(jìn)行優(yōu)化循環(huán)，直到系數(shù)穩(wěn)定為止。計算在λ的格點值上進(jìn)行，而λ的確定可以通過k折交叉驗證(k-fold cross-validation)得出。

R軟件glmnet包的基本語法

R軟件中可使用glmnet包實現(xiàn)懲罰廣義線性模型[8]。glmnet能夠處理各種數(shù)據(jù)類型，擬合線性、多項式、logistic和Poisson回歸等廣義線性模型。它的正則化路徑是在λ的格點值上計算LASSO或者彈性網(wǎng)的懲罰項。主要函數(shù)為：

glmnet(x,y,family=c(“gaussian”,“binomial”,“Poisson”,“multinomial”,“mgaussian”),alpha=1,…)

上述語句可擬合廣義線性模型，得到正則化的路徑；其中，x為自變量矩陣，y為因變量，family中指定因變量的分布類型，默認(rèn)各類分布的典型鏈接函數(shù)，如上述選擇項分別對應(yīng)正態(tài)分布(恒等鏈接)，二項分布(logit鏈接)，Poisson分布(log鏈接)，多項分布(累積logit鏈接)和多元正態(tài)分布(恒等鏈接)；alpha是彈性網(wǎng)的混合參數(shù)，取值范圍為[0,1]。當(dāng)設(shè)置alpha=1時擬合LASSO(缺省值)，alpha=0時擬合嶺估計；在(0,1)區(qū)間內(nèi)擬合彈性網(wǎng)。實際應(yīng)用中，alpha取值的選擇也可借助cv.glmnet函數(shù)，通過設(shè)定不同的alpha值進(jìn)行交叉驗證：

cv.glmnet(x,y,nfolds,alpha,…)

上述函數(shù)可完成模型的K折交叉驗證，除選擇alpha值外，還主要用于確定模型的調(diào)整參數(shù)λ。其中x和y意義同前，nfolds指定交叉驗證的K值，表示將所有觀測隨機(jī)分為K份，循環(huán)采用K-1份作為訓(xùn)練樣本，剩余1份作為驗證樣本，默認(rèn)值為10，最小可設(shè)定為3。

模型擬合后，可用coef和predict函數(shù)顯示懲罰模型的參數(shù)估計值及據(jù)此模型進(jìn)行預(yù)測：

coef(object,s=c(“l(fā)ambda.1se”,“l(fā)ambda.min”)or null,…)

predict(object,newx,s=c(“l(fā)ambda.1se”,“l(fā)ambda.min”)or null,…)

上述函數(shù)中，如對參數(shù)s進(jìn)行設(shè)定，object定義為cv.glmnet交叉驗證選擇的模型，否則為glmnet估計的所有模型；s中l(wèi)ambda.min為交叉驗證結(jié)果中平均誤差最小的模型對應(yīng)的值，lambda.1se則為最小平均誤差1倍標(biāo)準(zhǔn)誤內(nèi)對應(yīng)λ中的最大值；一般而言，定義Lambda.min輸出最優(yōu)模型(平均誤差最小)，但可能存在模型復(fù)雜及過擬合的問題；定義Lambda.1se可輸出更簡潔的模型，但較最優(yōu)模型存在相對誤差。此外，predict函數(shù)中的newx為參與預(yù)測的自變量矩陣。

實例分析

為說明懲罰廣義線性模型及R軟件glmnet包在遺傳關(guān)聯(lián)研究中的應(yīng)用，本文采用病例對照設(shè)計的遺傳關(guān)聯(lián)研究模擬數(shù)據(jù)，探討不同的遺傳情景下各類懲罰模型的參數(shù)估計效果。數(shù)據(jù)模擬過程中，設(shè)置了連鎖不平衡(linkage disequilibrium，LD)參數(shù)，以反映群體中不同基因座上等位基因間的關(guān)聯(lián)。在模擬數(shù)據(jù)中LD定義為高維SNPs間的相關(guān)系數(shù)，分兩種情況：(1)SNPs變異間獨立(LD=0)；(2)SNPs間高度相關(guān)(LD=0.9)。病例組和對照組樣本量各為500，為便于說明glmnet包的輸出結(jié)果，模擬數(shù)據(jù)集中只包含20個SNPs，并隨機(jī)抽取其中5個(X1，X4，X5，X12，X15)設(shè)置為與因變量有關(guān)聯(lián)的SNPs，剩余15個設(shè)置為與因變量無關(guān)聯(lián)的噪聲變異；與疾病有關(guān)聯(lián)的變異效應(yīng)值ORs均設(shè)為2。模擬數(shù)據(jù)在R軟件中完成，具體模擬方法參見文獻(xiàn)[10]。數(shù)據(jù)集形式為1000個觀測行，22個變量列，如表1所示。

表1　數(shù)據(jù)集形式

實例分析時，先在R軟件中載入glmnet包，讀取R軟件工作目錄的數(shù)據(jù)，建立自變量和因變量矩陣；設(shè)定隨機(jī)數(shù)種子，可以保證結(jié)果的唯一性。例如對連鎖平衡數(shù)據(jù)集SAMPLE_LD0.csv擬合懲罰logistic回歸模型，具體程序為：

library(glmnet)#載入glmnet包

linkage<-read.csv(“SAMPLE_LD0.csv”)#讀入數(shù)據(jù)

linkage.x<-as.matrix(linkage[,3:22])#生成自變量矩陣

linkage.y<-as.matrix(linkage[,2])#生成因變量矩陣

set.seed(2015)#設(shè)定隨機(jī)數(shù)種子

link_glm<-glmnet(linkage.x,linkage.y,family=“binomial”)#擬合LASSO模型

plot(link_glm)#繪制系數(shù)正則化的路徑圖

由于因變量為二分類，故在glmnet函數(shù)的family選項中選擇binomial，alpha缺省時為1，擬合基于LASSO的懲罰logistic回歸；如定義為0或(0，1)之間的數(shù)值，則擬合基于嶺估計和彈性網(wǎng)的懲罰logistic回歸。Link_glm中包含不同取值時的所有模型。圖1為alpha=1，0和0.5時分別擬合LASSO，嶺估計和彈性網(wǎng)系數(shù)懲罰的收縮情況；其中橫坐標(biāo)為L1正則化，即系數(shù)的絕對值之和，又稱曼哈頓距離；縱坐標(biāo)為系數(shù)，上方數(shù)字是模型中保留的變量個數(shù)，圖中可看到系數(shù)在正則化路徑上的選擇。此例中LASSO和彈性網(wǎng)壓縮過程相似，可將部分系數(shù)連續(xù)壓縮為0，從而實現(xiàn)變量選擇的目的，而嶺估計壓縮系數(shù)的過程中仍保留所有變量在模型中。

為選擇適宜的λ值，使用cv.glmnet函數(shù)，通過交叉驗證觀察不同λ取值時的模型誤差，具體程序為：

link_cv <-cv.glmnet(linkage.x,linkage.y,nfolds=10,family=“binomial”)#10折交叉驗證

plot(link_cv)#繪制交叉驗證曲線圖

link_cv$lambda.min#最小誤差模型的lambda

link_cv$lambda.1se#1倍SE內(nèi)的最大lambda

圖1　參數(shù)估計變化趨勢圖

圖2為分別擬合LASSO，嶺估計和彈性網(wǎng)時取不同log(λ)時模型誤差的變化趨勢。其中左側(cè)虛線對應(yīng)平均誤差最小時(最優(yōu)模型)的log(λ)(lambda.min)，右側(cè)另一條虛線是其1倍SE時對應(yīng)的log(λ)(lambda.1se)，對應(yīng)更簡潔的模型；圖形上方數(shù)字為模型對應(yīng)的變量個數(shù)。本例中LASSO和彈性網(wǎng)各模型平均誤差接近，而嶺估計較高。如基于LASSO懲罰回歸，由于lambda.min和lambda.1se 對應(yīng)的模型誤差變化不大，更偏好簡潔模型時，選擇lambda.1se對應(yīng)的模型，此時λ值為0.0255，交叉驗證模型的平均誤差為1.3558。

選擇適宜的λ后，可提取最終模型參數(shù)估計的結(jié)果，包括回歸系數(shù)(beta)，自由度(df)，λ和決定系數(shù)R2(dev.ratio)等；其中決定系數(shù)是評判自變量對因變量解釋程度的重要指標(biāo)，可進(jìn)行模型之間的比較，具體程序為：

link_cv.best <-link_cv$glmnet.fit#提取最佳模型

link_cv.coef <-coef(link_cv$glmnet.fit,s=link_cv$lambda.1se)#提取模型系數(shù)

link_cv.coef #輸出模型系數(shù)

R<-link_cv.best$dev.ratio[which(link_cv.best$lambda==link_cv$lambda.1se)]#輸出R2

圖2　交叉驗證λ和模型誤差的關(guān)系

根據(jù)以上的語句，可輸出模型最終參數(shù)估計結(jié)果及R2。本例中，調(diào)整glmnet函數(shù)和cv.glmnet函數(shù)內(nèi)的alpha值依次為1,0.5,0，分別對連鎖平衡數(shù)據(jù)和連鎖不平衡數(shù)據(jù)擬合LASSO、彈性網(wǎng)和嶺估計，各模型參數(shù)估計值及R2值如表2所示。

表2　連鎖平衡與連鎖不平衡數(shù)據(jù)擬合模型的參數(shù)估計值及R2

由表2可得，本例中，對連鎖平衡和連鎖不平衡數(shù)據(jù)，LASSO均給出稀疏解，較好地選擇有關(guān)聯(lián)SNPs而剔除無關(guān)聯(lián)變量；除關(guān)聯(lián)SNPs外，彈性網(wǎng)多保留一個無關(guān)聯(lián)SNP；而嶺估計因無法把參數(shù)壓縮為0，所有變量都保留在模型中，模型的復(fù)雜度高，但解釋度未增加。

討　　論

由于生物學(xué)測序技術(shù)的迅猛發(fā)展，遺傳關(guān)聯(lián)數(shù)據(jù)趨于高維，需要處理的信息量越來越大且存在大量的冗余信息和噪音，選擇顯著變量、擬合較簡單卻解釋程度好的模型是研究者面臨的共同問題，相應(yīng)的統(tǒng)計分析方法也越來越被關(guān)注?；趹土P的廣義線性模型在處理共線性，降維以及變量選擇方面有著獨特的優(yōu)勢，在各類科研領(lǐng)域尤其是遺傳流行病學(xué)研究中的應(yīng)用也隨之逐步興盛起來[3]。本文以遺傳關(guān)聯(lián)分析中常用的病例對照設(shè)計模擬數(shù)據(jù)為例，擬合LASSO，嶺估計和彈性網(wǎng)三種懲罰模型，運用glmnet包進(jìn)行關(guān)聯(lián)性SNPs位點選擇，說明懲罰模型在遺傳關(guān)聯(lián)分析中應(yīng)用的可行性及科學(xué)性。

實例分析可見，嶺估計能夠有效克服自變量間的高度相關(guān)，得到穩(wěn)定的模型，但是它沒有稀疏性，不能有效地篩選自變量；但有研究表明，在嶺估計基礎(chǔ)上引入“Boosting”后可得到與LASSO同樣的估計，實現(xiàn)篩選變量的目的[11]。LASSO即L1懲罰可很好地降低維度，輕微解決共線性問題，但也有研究表明，LASSO有時可能過度壓縮參數(shù)[6]。彈性網(wǎng)對系數(shù)的懲罰介于嶺估計與LASSO之間，在本例結(jié)果中的表現(xiàn)更接近LASSO，且和LASSO一樣能同時進(jìn)行模型選擇和參數(shù)估計。以往研究顯示，彈性網(wǎng)一般不會過度壓縮參數(shù)，能有效處理組群效應(yīng)(group effect)，特別是能直接處理變量數(shù)遠(yuǎn)遠(yuǎn)大于樣本數(shù)的問題，這些是LASSO所不能達(dá)到的[12]。因此在實際分析中，還需根據(jù)具體情況選擇適宜的方法，如自變量弱共線時可考慮LASSO，強(qiáng)共線、高維小樣本數(shù)據(jù)時傾向于彈性網(wǎng)。本研究為便于說明R軟件結(jié)果，只采用了單個模擬數(shù)據(jù)集，無法體現(xiàn)懲罰模型更進(jìn)一步的性質(zhì)，如不同樣本量和連鎖不平衡參數(shù)時參數(shù)估計的準(zhǔn)確度、混雜無關(guān)聯(lián)變異比例對變量選擇的影響等；更深入研究和比較各種懲罰方法，還需設(shè)置更多遺傳情境的模擬研究和實例應(yīng)用。

glmnet包作為R數(shù)據(jù)挖掘中的三駕馬車之一，過程靈活，功能性強(qiáng)，可以擬合各類基于懲罰的廣義線性模型，通過最小化懲罰目標(biāo)函數(shù)得到參數(shù)估計值。但是與一般的廣義線性模型不同，glmnet包無法直接計算Fisher信息陣進(jìn)而估計參數(shù)的方差。有研究者提出可采用bootstrap方法計算估計值的標(biāo)準(zhǔn)誤[6]，或采用Bayesian LASSO獲得有效的β及標(biāo)準(zhǔn)誤[13]。還有學(xué)者提出可在似然比中設(shè)定“三明治”(sandwich formula)用于協(xié)方差估計[14]，這些方法及應(yīng)用程序?qū)⑦M(jìn)一步補充和完善glmnet包的功能。此外，在研究SNPs與疾病預(yù)后常用的生存時間資料時，glmnet包還可擬合基于三種懲罰的Cox模型。除此之外，R中還有許多其他擬合LASSO的相關(guān)軟件包：如較早時基于最小角回歸擬合傳統(tǒng)線性模型的lars；運用同倫算法擬合廣義線性模型的LASSO 2；還有和glmnet類似的penalized，也可選擇L1，L2或彈性網(wǎng)懲罰擬合各種廣義線性模型；而基于LARS-EN算法的elasticnet可以進(jìn)行彈性網(wǎng)懲罰，能夠估計稀疏主成分。實際應(yīng)用中，研究者可結(jié)合數(shù)據(jù)特征和實際需求加以選擇。

總之，運用R對高維數(shù)據(jù)建模時，懲罰廣義線性模型克服了傳統(tǒng)線性模型存在的預(yù)測精度和模型解釋度較差的缺陷，基于懲罰的變量選擇方法能夠在選擇變量的同時得到參數(shù)估計，加之計算量小，因此在處理高維數(shù)據(jù)時，有著傳統(tǒng)變量選擇方法無可比擬的優(yōu)越性[15]。相信隨著懲罰廣義線性模型及軟件的推廣使用，它將在高通量遺傳關(guān)聯(lián)分析中發(fā)揮極其重要的作用。

[1]Ayers KL,Cordell HJ.Analysis of Genetic Analysis Workshop 18 data with gene-based penalized regression.BioMed Central Ltd,2014,S43.

[2]Zhang Z,Ersoz E,Lai C,et al.Mixed linear model approach adapted for genome-wide association studies.Nature genetics,2010,42(4):355-360.

[3]龔建朝.Lasso及其相關(guān)方法在廣義線性模型模型選擇中的應(yīng)用.中南大學(xué),2008.

[4]馬文浩.各種L_q懲罰在變量選擇中的應(yīng)用及其比較.山東大學(xué),2012.

[5]Hoerl AE,Kennard RW.Ridge Regression:Biased Estimation for Nonorthogonal Problems.Technometrics,1970(12):55-67.

[6]Tibshirani R.Regression shrinkage and selection via the lasso.Journal of the Royal Statistical Society.Series B(Methodological),1996:267-288.

[7]Zou H,Hastie T.Regularization and variable selection via the elastic net.J R Statist,2005(67):301-320.

[8]Friedman J,Hastie T,Simon N,et al.glmnet:Lasso and elastic-net regularized generalized linear models.CRAN，2009.

[9]Friedman J,Hastie T,Tibshirani R.Regularization paths for generalized linear models via coordinate descent.Journal of Statistical Software,2009,33(1):1-22.

[10]梁融.稀有變異關(guān)聯(lián)性分析中折疊與非折疊法的模擬比較研究.廣東藥學(xué)院,2014.

[11]Tutz G,Binder H.Boosting ridge regression.Computational Statistics & Data Analysis,2007,51(12):6044-6059.

[12]Bhlmann P,Sara VDG.Statistics for high-dimensional data:methods,theory and applications.Springer Science & Business Media,2011.

[13]Casella TPG.The Bayesian Lasso.Journal of the American Statistical Association,2008,103(482):681-686.

[14]Fan J,Li R.Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties.Journal of the American Statistical Association,2001,96(456):1348-1360.

[15]張秀秀,王慧,田雙雙,等.高維數(shù)據(jù)回歸分析中基于LASSO的自變量選擇.中國衛(wèi)生統(tǒng)計,2013,30(6):922-926.

(責(zé)任編輯：鄧妍)

Application of Penalized Generalized Linear Model in Genetic Association Study and its Software Implementation in R

Zhang Junguo,Liu Li,Li Lixia,et al

(Department of Epidemiology and Biostatistics,School of Public Health,Guangdong Pharmaceutical University(510310),Guangzhou)

ObjectiveWith the development of high-dimensional data in genetic association studies,this study was aimed to depict the application of penalized generalized linear model based on ridge,LASSO and elastic net in genetic association study and its software implementation.MethodsWe introduced the penalized generalized linear model in detail,and performed the analyses in the simulated SNPs data in linkage equilibrium and linkage disequilibrium,respectively.At last,we used the ‘glmnet’ package in R software to fit the penalized generalized linear model with the example of logistic model.ResultsFor both of the SNPs data in linkage equilibrium and linkage disequilibrium,LASSO and elastic net could worked out sparse solution and effectively detected out the associated SNPs,as well as fairly excluded the non-associated variables.By contrast,the ridge model included all variables,and thus promoted the complexity of model,but meanwhile its explanation was not added.ConclusionFor the high-dimensional association data,LASSO and Elastic net can achieve effective dimensionality reduction,variables selection,and parameters estimation simultaneously,and thus reduce the complexity of the model.R glmnet package can flexibly fit different types of generalized linear model.Therefore,it can be wildly used in high-throughput genetic association studies.

Penalized generalized linear model;Genetic association studies;LASSO;Elastic net;Glmnet package

國家自然科學(xué)基金(81302493)；廣東省科技廳社會發(fā)展基金(2014A020212307)；廣東省自然科學(xué)基金(S2013040013590)

郜艷暉，E-mail:gao_yanhui@163.com