哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)
謝宏宇 侯 艷 李 康△
?
基于正則化回歸的組學(xué)數(shù)據(jù)變量篩選方法*
哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)
謝宏宇侯艷李康△
近年來,隨著各種生物檢測(cè)技術(shù)的發(fā)展,醫(yī)學(xué)研究中出現(xiàn)了各種高通量數(shù)據(jù),如基因組、蛋白質(zhì)組和代謝組學(xué)數(shù)據(jù)等,變量選擇是生物標(biāo)志物識(shí)別和建立分類模型的重要環(huán)節(jié),由于高維組學(xué)數(shù)據(jù)中的絕大多數(shù)變量對(duì)分類并不起作用,并且存在多重共線性、模型過擬合等問題,傳統(tǒng)的基于最小二乘方法估計(jì)的線性回歸并不適用于高維數(shù)據(jù)[1]。在高維組學(xué)數(shù)據(jù)特征變量篩選過程中,由于變量數(shù)目很多,子集篩選方法計(jì)算量巨大,并可能由于選擇不同的篩選準(zhǔn)則導(dǎo)致篩選的結(jié)果有很大差異[2];維數(shù)縮減方法雖然能夠避免計(jì)算量大的問題,但是由于模型中的變量不再是原始變量,模型的可解釋性差;而正則化回歸方法由于在解回歸方程的過程中可以同時(shí)實(shí)現(xiàn)參數(shù)估計(jì)和變量篩選,且計(jì)算速度快,對(duì)變量數(shù)目沒有限制,因此受到研究者的關(guān)注[3]。這類方法不僅能夠用于單一組學(xué)數(shù)據(jù)的變量篩選,同時(shí)也能拓展到多組學(xué)數(shù)據(jù)融合的情況,因此在實(shí)際中具有很好的應(yīng)用前景。本文將對(duì)正則化回歸方法及在高維組學(xué)數(shù)據(jù)中的應(yīng)用做一綜述。
正則化是指在原有的損失函數(shù)的基礎(chǔ)上增加懲罰回歸系數(shù)的正則項(xiàng)。記β=(β1,β2,…,βm)為回歸系數(shù),m為總自變量的個(gè)數(shù),則在線性回歸中,通過最小化損失函數(shù)可以得出對(duì)應(yīng)模型的回歸系數(shù)估計(jì)值
(1)
其中‖·‖2是L2范數(shù),表示向量各元素平方和的平方根,該式表示取右端函數(shù)最小值的模型參數(shù),這實(shí)際是傳統(tǒng)的最小二乘估計(jì)。當(dāng)變量個(gè)數(shù)較多時(shí),利用該式估計(jì)得出的回歸模型存在過擬合的風(fēng)險(xiǎn),正則化則可以在保留所有特征變量的情況下,避免過擬合的發(fā)生,其基本原理是通過增加的正則項(xiàng),減少所有特征變量回歸系數(shù)估計(jì)值的數(shù)量級(jí),具體的表現(xiàn)形式如下:
(2)
其中,λP(β)表示正則化項(xiàng),λ為正則化參數(shù),P(β)為回歸系數(shù)的懲罰函數(shù),主要目的是用來平衡模型對(duì)樣本數(shù)據(jù)的擬合程度以及回歸模型的預(yù)測(cè)能力。在正則化項(xiàng)中,如果正則化參數(shù)設(shè)定較大會(huì)使得每個(gè)回歸系數(shù)估計(jì)值偏??;如果回歸系數(shù)估計(jì)值小到一定程度時(shí),相當(dāng)于因變量只等于常數(shù)項(xiàng),類似于擬合了一條水平直線,導(dǎo)致欠擬合,產(chǎn)生過高的偏差。如果模型中涉及到高階項(xiàng),則回歸系數(shù)的估計(jì)值越小,對(duì)應(yīng)的曲線越光滑,從而使函數(shù)得到簡(jiǎn)化,實(shí)際中需要選擇合適的正則化參數(shù)值。目前正則化參數(shù)的選擇可以通過偏差原理、Engl誤差極小原理、Hansen 的L曲線準(zhǔn)則、擬最優(yōu)準(zhǔn)則和交叉驗(yàn)證等方法進(jìn)行確定[4]。回歸模型的系數(shù)估計(jì)值可以通過梯度下降等方法進(jìn)行求解。
1.嶺回歸
嶺回歸(ridge regression)方法由Hoerl 和Kennard提出[5],其基本思想是在傳統(tǒng)最小化殘差平方和基礎(chǔ)上加入回歸系數(shù)的L2范數(shù)懲罰項(xiàng)從而收縮回歸系數(shù)。最小化回歸系數(shù)的L2范數(shù),會(huì)使稀疏矩陣中每個(gè)元素的值都很小,但并不一定為0。回歸系數(shù)估計(jì)值的表達(dá)式如下:
(3)
其中λ是正則化參數(shù)。由于L2范數(shù)可以收縮回歸系數(shù)估計(jì)值,因此能夠在一定程度上避免模型的過擬合。嶺回歸的主要特點(diǎn)是通過L2范數(shù)對(duì)回歸系數(shù)的連續(xù)收縮,能夠使每個(gè)變量的系數(shù)變小,從而通過損失無偏性提高了模型的預(yù)測(cè)能力。主要缺點(diǎn)是,嶺回歸將所有的預(yù)測(cè)變量均保留在模型中,因此在分析高維組學(xué)數(shù)據(jù)時(shí)會(huì)導(dǎo)致模型的可解釋性較差。
2.lasso回歸
Tibshirani于1996年提出了基于線性回歸的最小化的絕對(duì)收縮和選擇算子(least absolute shrinkage and selection operator,lasso)來收縮回歸系數(shù),這種方法在損失函數(shù)中增加了回歸系數(shù)的L1范數(shù)懲罰項(xiàng),表示為‖·‖1,代表向量中各個(gè)元素絕對(duì)值之和。在回歸系數(shù)的絕對(duì)值之和小于一個(gè)常數(shù)的約束條件下,使殘差平方和最小化,能夠使部分回歸系數(shù)等于0,同時(shí)實(shí)現(xiàn)回歸系數(shù)收縮和變量篩選,從而提高了模型的可解釋性[6]。回歸系數(shù)估計(jì)值的表達(dá)式如下:
(4)
隨著正則化參數(shù)λ的增大,lasso方法能夠不斷地縮小回歸系數(shù)的估計(jì)值,使其趨近于0,實(shí)現(xiàn)回歸系數(shù)的稀疏化。在高維組學(xué)數(shù)據(jù)中,最常用于估計(jì)lasso回歸系數(shù)的方法為最小角算法(least angle regression,LARS)[7],這種算法相對(duì)于最小二乘回歸能夠很好地解決lasso回歸的計(jì)算問題。lasso回歸存在一定的局限性,即在自變量個(gè)數(shù)m遠(yuǎn)大于樣本量n時(shí),只能保證lasso回歸中最多選擇n個(gè)變量;同時(shí),如果一組變量高度相關(guān)時(shí),這種算法只傾向于選擇其中之一,而不關(guān)心選擇的究竟是哪個(gè)變量[5]。
3.自適應(yīng)lasso回歸
Zou(2006)發(fā)現(xiàn)lasso回歸中L1范數(shù)懲罰項(xiàng)對(duì)所有回歸系數(shù)懲罰強(qiáng)度相同,從而導(dǎo)致了回歸系數(shù)估計(jì)值不具有漸進(jìn)正態(tài)性[8]。另外,lasso回歸只有在兩種特定條件下的變量篩選才具有相合性[5]。為了解決這個(gè)問題,Zou(2006)提出了一種新的lasso回歸方法,即自適應(yīng)lasso回歸,回歸系數(shù)估計(jì)值的表達(dá)式為
(5)
其中wj代表回歸系數(shù)的權(quán)重。該方法的特點(diǎn)是可以對(duì)不同的系數(shù)設(shè)置不同的懲罰權(quán)重,從而改進(jìn)其估計(jì)值的準(zhǔn)確性。如果依賴于數(shù)據(jù)本身對(duì)權(quán)重做出恰當(dāng)?shù)倪x擇,則自適應(yīng)lasso回歸具有相合性和漸進(jìn)正態(tài)性,并且能夠避免局部最優(yōu)化的問題[9]。Breheny在2013年將自適應(yīng)lasso回歸方法用于微陣列數(shù)據(jù)進(jìn)行變量篩選[10],并采用預(yù)測(cè)誤差均方作為評(píng)價(jià)指標(biāo)與lasso回歸進(jìn)行比較,結(jié)果表明自適應(yīng)lasso回歸方法的預(yù)測(cè)誤差均方小于lasso回歸,并且篩選出的差異變量更少[11]。
4.樸素彈性網(wǎng)和彈性網(wǎng)算法
高維組學(xué)數(shù)據(jù)往往具有高度相關(guān)性和分組特征(例如來自于同一通路的基因),如前所述lasso回歸方法針對(duì)以上兩種情況進(jìn)行變量篩選效果不理想。因此Zou于2003年提出了彈性網(wǎng)算法(elastic net)[12],該方法既能夠同時(shí)實(shí)現(xiàn)變量的自動(dòng)篩選和回歸系數(shù)的連續(xù)收縮,又能夠保證選擇出同一分組內(nèi)與因變量相關(guān)性大的變量。樸素彈性網(wǎng)(naive elastic net)算法是最基礎(chǔ)的彈性網(wǎng)算法,主要是將lasso回歸的懲罰項(xiàng)和嶺懲罰項(xiàng)相結(jié)合,其表達(dá)式為
(6)
其中,λ1和λ2均為非負(fù)的正則化參數(shù)。若記α=λ2/(λ1+λ2),上式等價(jià)于
(7)
并且(1-α)‖β‖1+α‖β‖2≤t
其中,t為一個(gè)常數(shù)界值,(1-α)‖β‖1+α‖β‖2稱為彈性網(wǎng)懲罰,α∈[0,1]。當(dāng)α=0時(shí),該式為lasso回歸懲罰項(xiàng);當(dāng)α=1時(shí),為嶺回歸懲罰項(xiàng)。
樸素彈性網(wǎng)算法的參數(shù)估計(jì)分為兩個(gè)階段:首先固定λ2找到嶺回歸系數(shù),然后通過λ1進(jìn)行系數(shù)壓縮。雖然樸素彈性網(wǎng)算法能夠克服傳統(tǒng)lasso回歸的部分不足,但模擬實(shí)驗(yàn)表明只有當(dāng)它接近嶺回歸或者lasso回歸時(shí),才能獲得較理想的變量篩選結(jié)果。因此,Zou于2005年又提出了對(duì)樸素彈性網(wǎng)系數(shù)進(jìn)行重縮放,這種方法即為目前的彈性網(wǎng)算法。當(dāng)對(duì)預(yù)測(cè)變量進(jìn)行標(biāo)準(zhǔn)化后,彈性網(wǎng)方法的回歸系數(shù)與樸素彈性網(wǎng)的回歸系數(shù)之間具有如下關(guān)系:
(8)
其中1+λ2為收縮因子。
研究結(jié)果表明,彈性網(wǎng)算法與lasso回歸和嶺回歸相比具有較好的篩選變量的性能,L1范數(shù)可以實(shí)現(xiàn)自動(dòng)變量篩選,L2范數(shù)可以實(shí)現(xiàn)連續(xù)收縮,尤其在自變量之間存在較強(qiáng)的相關(guān)性時(shí),彈性網(wǎng)算法能夠明顯的提高預(yù)測(cè)的準(zhǔn)確性[3]。Zou將幾種變量篩選方法應(yīng)用于白血病患者的基因表達(dá)數(shù)據(jù),目的是篩選用于診斷和預(yù)測(cè)白血病分型的基因。結(jié)果表明,彈性網(wǎng)算法構(gòu)建的模型分類效果優(yōu)于支持向量機(jī)和懲罰logistic回歸等方法,并且能對(duì)組內(nèi)基因進(jìn)行篩選。由于彈性網(wǎng)算法估計(jì)出的系數(shù)不具有相合性和漸進(jìn)正態(tài)性[5,8],因此Zou于2009年提出將自適應(yīng)加權(quán)L1懲罰納入到彈性網(wǎng)算法正則項(xiàng)中提高估計(jì)準(zhǔn)確性,即自適應(yīng)彈性網(wǎng)算法(adaptive elastic-net),該方法可以將其視為彈性網(wǎng)算法和自適應(yīng)lasso的結(jié)合,具有相合性和漸進(jìn)正態(tài)性[13]。
正則化參數(shù)λ決定了模型中回歸系數(shù)估計(jì)值的大小和稀疏化的程度。確定正則化參數(shù)的基本方法有交叉驗(yàn)證[14]、貝葉斯信息準(zhǔn)則(BIC)[15]、Cp統(tǒng)計(jì)量[7]和赤池信息量準(zhǔn)則(AIC)[16]。Zou于2007年從變量篩選的角度證明了BIC相對(duì)于其他方法更適用于參數(shù)值的選擇,該方法能夠產(chǎn)生一個(gè)更加稀疏的模型[12]。Chen等認(rèn)為使用BIC準(zhǔn)則在高維數(shù)據(jù)中篩選變量的標(biāo)準(zhǔn)具有一定的任意性,因此提出了擴(kuò)展的BIC方法(EBIC),這種方法既考慮了未知參數(shù)的個(gè)數(shù),也考慮了模型空間的復(fù)雜性,并且能夠更加嚴(yán)格地控制差異變量錯(cuò)誤發(fā)現(xiàn)率[17]。
5.分組lasso回歸
(9)
其中L為變量的組別數(shù),l=1,2,…,L,X(l)代表組l中與因變量有關(guān)的X列的子矩陣,β(l)是組的系數(shù)向量,pl是第l組中包含的變量個(gè)數(shù)。這種方法利用了‖β(l)‖2在β(l)=0處不可微的性質(zhì),將該組從模型中剔除。其主要思想是篩選出對(duì)因變量有影響的特征組,同時(shí)通過選擇合適的參數(shù)λ調(diào)整組別個(gè)數(shù),λ值越大,對(duì)各分組作用的懲罰越大,則模型中保留的組數(shù)越少。雖然分組lasso回歸能夠?qū)崿F(xiàn)對(duì)組別的篩選,但是只能篩選出模型中整個(gè)組內(nèi)的變量回歸系數(shù)β(l)=0的特征組,這一缺點(diǎn)限制了其應(yīng)用。當(dāng)每個(gè)特征組內(nèi)只包含一個(gè)自變量時(shí),則該方法即為傳統(tǒng)的lasso方法。
6.稀疏組lasso回歸
實(shí)際研究中不僅僅需要實(shí)現(xiàn)組別的稀疏化,同時(shí)還需要實(shí)現(xiàn)組內(nèi)變量的稀疏化,例如,研究者識(shí)別感興趣基因通路,同時(shí)對(duì)該條基因通路中的關(guān)鍵基因進(jìn)行篩選。因此Simon(2010)提出將lasso回歸和分組lasso回歸相結(jié)合,引進(jìn)了稀疏組lasso(sparse-group lasso)回歸[19],其表達(dá)式為
(10)
稀疏組lasso回歸的方法與彈性網(wǎng)方法相似,不同點(diǎn)在于該種方法是利用在懲罰回歸系數(shù)為0時(shí)不可微的性質(zhì),將稀疏為0的組別從模型中去除,實(shí)現(xiàn)組間稀疏化,而彈性網(wǎng)方法則保留了所有的分組。Simon(2013)將lasso回歸、分組lasso回歸和稀疏組lasso回歸的方法應(yīng)用于乳腺癌患者的基因表達(dá)數(shù)據(jù)中,并比較了三種方法的篩選效果。結(jié)果表明稀疏組lasso回歸的變量篩選性能優(yōu)于lasso回歸和分組lasso回歸:稀疏組lasso回歸的分類正確率達(dá)到70%,而分組lasso回歸和lasso回歸的分類正確率分別為60%和53%。由于在癌癥數(shù)據(jù)中添加分組的信息對(duì)于分類非常有意義,同時(shí)分組信息可以幫助更加深入的了解生物學(xué)機(jī)制,因此對(duì)于癌癥數(shù)據(jù)的分析,稀疏組lasso回歸的方法有很大的優(yōu)勢(shì)[19]。
傳統(tǒng)的變量篩選方法一般均可應(yīng)用于單一組學(xué)數(shù)據(jù)的變量篩選,如基于回歸、基于機(jī)器學(xué)習(xí)和基于網(wǎng)絡(luò)的方法等,目前應(yīng)用于多組學(xué)融合數(shù)據(jù)變量篩選方法相對(duì)較少,而基于正則化的變量篩選方法可以實(shí)現(xiàn)多組學(xué)數(shù)據(jù)的融合和變量篩選。
1.稀疏廣義典型相關(guān)分析
典型相關(guān)分析(canonical correlation analysis,CCA)是用于研究兩組變量之間關(guān)系的常用方法。Tenenhaus于2011年提出了正則化的廣義典型相關(guān)分析(regularized generalized canoncial analysis,RGCCA)方法,該方法可用于分析三個(gè)或者更多的變量集合間的關(guān)系[20]。RGCCA是一種基于主成分分析的方法,用于研究多個(gè)數(shù)據(jù)集中變量之間的關(guān)系。RGCCA成分的性質(zhì)和解釋性受每組變量之間有用性和相關(guān)性的影響。RGCCA主要基于使多個(gè)數(shù)據(jù)樣本中生成新的綜合變量的相關(guān)程度最大化的思想進(jìn)行求解。
實(shí)際中,在每組變量中識(shí)別出在組間關(guān)系中起顯著作用的變量子集非常重要,因此Tenenhaus在2014年提出了稀疏的廣義典型相關(guān)分析(sparse generalized canonical correlation analysis,SGCCA),這種方法通過對(duì)外部權(quán)重向量加上L1懲罰,在同一方法中結(jié)合了RGCCA和L1懲罰項(xiàng)[20]。Tenenhaus將SGCCA方法應(yīng)用到兒科神經(jīng)膠質(zhì)瘤數(shù)據(jù),結(jié)果表明:與RGCCA比較,SGCCA方法能夠篩選出在組間相關(guān)作用中具有更小差異的變量組合[21]。
2.稀疏偏最小二乘回歸
偏最小二乘(PLS)的變量篩選方法已經(jīng)成功地應(yīng)用于代謝組學(xué)數(shù)據(jù)。其主要原理是分別在自變量和因變量中提取出成分,使各自提取出的成分盡可能多的解釋各數(shù)據(jù)的變異信息,同時(shí)使提取成分的相關(guān)程度達(dá)到最大。
Cao(2008)在此基礎(chǔ)上提出了稀疏偏最小二乘回歸(sparse least squares regression,SPLS)的變量篩選方法,該方法能夠同時(shí)實(shí)現(xiàn)數(shù)據(jù)的整合和變量篩選,其主要思想是在PLS的基礎(chǔ)上,通過Q2值作為評(píng)價(jià)指標(biāo)對(duì)構(gòu)建模型的成分?jǐn)?shù)量進(jìn)行選擇,同時(shí)對(duì)每個(gè)成分加上lasso懲罰,實(shí)現(xiàn)變量篩選。研究表明:該方法應(yīng)用于高維數(shù)據(jù)集分析時(shí),與PLS相比具有更高的穩(wěn)定性,能夠更好地進(jìn)行變量篩選[22]。
3.結(jié)構(gòu)組稀疏算法
多尺度數(shù)據(jù)分析的關(guān)鍵性問題是,數(shù)據(jù)結(jié)構(gòu)異質(zhì)性整合和特征變量篩選的穩(wěn)定性?;诮Y(jié)構(gòu)組稀疏算法(structure grouping sparsity,SGS)的多尺度數(shù)據(jù)變量篩選方法的目的是根據(jù)實(shí)際數(shù)據(jù)給出一個(gè)可以解釋和預(yù)測(cè)的模型。其主要思想是,根據(jù)實(shí)際數(shù)據(jù)建立應(yīng)變量Y與自變量X=(X1,X2,…,Xm)關(guān)系的廣義線性模型(可拓展至非線性),實(shí)現(xiàn)對(duì)不同來源的異質(zhì)性數(shù)據(jù)在不同水平上進(jìn)行組間和組內(nèi)特征篩選。其表達(dá)式為:
(11)
基于正則化回歸的變量篩選方法,克服了傳統(tǒng)變量篩選方法的不足,且隨著研究的深入,需要不斷的更新和發(fā)展。該方法的發(fā)展一直圍繞著擬合較好的模型應(yīng)該具有預(yù)測(cè)準(zhǔn)確度高、模型的可解釋性強(qiáng)的特點(diǎn);本著模型本身具有優(yōu)良的參數(shù)估計(jì)性質(zhì),即無偏性、有效性、相合性和漸進(jìn)正態(tài)性?;谡齽t化的變量篩選方法不僅能夠應(yīng)用于單一組學(xué)數(shù)據(jù)的變量篩選,也能夠應(yīng)用于多組學(xué)數(shù)據(jù)的融合和變量篩選。然而,這種方法的懲罰項(xiàng)選擇及其統(tǒng)計(jì)性質(zhì),以及參數(shù)求解等問題都有待進(jìn)一步研究。展望未來,高維組學(xué)的數(shù)據(jù)研究將實(shí)現(xiàn)跨組學(xué)的超高維變量篩選,從而更全面的研究疾病的發(fā)生機(jī)制,因此這類方法將會(huì)具有較好的前景。
[1]趙奕林,朱真峰,周清雷.適用于大規(guī)模高維多類別數(shù)據(jù)分類的并行非線性最小二乘分類器.小型微型計(jì)算機(jī)系統(tǒng),2014,3:579-583.
[2]Daniel PB,Pierluigi C.Introduction to the theory of complexity.Prentice Hall.ISBN 0-13-915380-2,1994.
[3]Zou H,Hastie T.Regularization and variable selection via the elastic net.J.R.Statist.Soc.B,2005,67:301-320.
[4]閔濤,葛寧國,黃娟,等.正則參數(shù)求解的微分進(jìn)化算法.應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué)學(xué)報(bào),2010,24(2):23-27.
[5]Hoerl A,Kennar R.Ridge regression.In Encyclopedia of Statistical Sciences,1998,8:129-136.
[6]Tibshirani R.Regression shrinkage and selection via the lasso.J.R.Statist.Soc.B,1996,58:267-288.
[7]Efron B,Hastie T,Johnstone I,et al.Least angle regression.Ann.Statist.,2004,32:407-499.
[8]Zou H.The adaptive lasso and its oracle properties.Journal of the American Statistical Association,2006,101:1418-1429.
[9]Li ZT,Mikko J,Sillanpaa.Overview of lasso-related penalized regression methods for quantitative trait mapping and genomic selection.Theor Appl Genet,2012,125:419-435.
[10]Scheetz T,Kim K,et al.Regulation of gene expression in the mammalian eye and its revevance to eye disease.Proc.Natl.Acad.Sci,2006,103:14429-14434.
[11]Patrick B,Jian H.Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors.Stat Comput,2015,25:173-187.
[12]Zou H,Hastie T.Regression shrinkage and selection via the elastic net,with application to microarrays,2003,1-26.
[13]Zou H,Zhang H.On the adaptive elastic-net with a diverging number of parameters.Ann.Statist.,2009,37:1733-1751.
[14]Hastie,Tibshirani R,Friedman JH.The elements of statistical learning.Springer,New York,2009.
[15]Zou H,Hastie T,Tibshirani R.On the “degrees of freedom” of the lasso.Ann Stat,2007,35:2173-2192.
[16]Akaike H.New look at the statistical model identification.IEEE T Autom Contr,1974,19:716-723.
[17]Chen J,Chen Z.Extended Bayesian information criteria for model selection with large model spaces.Biometrika,2008,95:759-771.
[18]Yuan M,Lin Y.Model selection and estimation in regression with grouped variables.Journal of the Royal Statistical Society,Series B,2007,68(1):49-67.
[19]Simon N,Friedman J,et al.A Sparse-Group lasso.Journal of computational and Graphical Statistics,2013,22:231-245.
[20]Tenenhaus A,Tenenhuas M.Regularized generalized canonical analysis.Psychometrika,2011,76:257-284.
[21]Tenenhaus.Variable selection for generalized canonical correlation analysis.Biostatistics,2014:1-15.
[22]Le Cao KA,Rossouw D,et al.A sparse PLS for variable selection when integrating omics data.Stat Appl Genet Mol Biol,2008,7(1):1-32.
(責(zé)任編輯:郭海強(qiáng))
國家自然科學(xué)基金資助(81573256,81473072);黑龍江省博士后資助經(jīng)費(fèi)(LBH-Z14174)
李康,E-mail:likang@ems.hrbmu.edu.cn