北京協(xié)和醫(yī)學(xué)院中國醫(yī)學(xué)科學(xué)院國家心血管病中心阜外醫(yī)院醫(yī)學(xué)統(tǒng)計(jì)部（100037） 嚴(yán)若華 李 衛(wèi)

?

Cox回歸模型比例風(fēng)險(xiǎn)假定的檢驗(yàn)方法研究

北京協(xié)和醫(yī)學(xué)院中國醫(yī)學(xué)科學(xué)院國家心血管病中心阜外醫(yī)院醫(yī)學(xué)統(tǒng)計(jì)部（100037） 嚴(yán)若華 李 衛(wèi)△

生存分析（survival analysis）是一系列處理“事件發(fā)生時(shí)間”變量的統(tǒng)計(jì)方法總稱［1］，常用于研究疾病的發(fā)生、轉(zhuǎn)歸、痊愈和死亡。Cox比例風(fēng)險(xiǎn)回歸模型（Cox proportional hazards regression model）［2］作為生存分析中最重要的多因素分析方法之一，被廣泛應(yīng)用于臨床隨訪資料的危險(xiǎn)因素篩選及預(yù)測。

Cox回歸模型將風(fēng)險(xiǎn)函數(shù)表達(dá)為基準(zhǔn)風(fēng)險(xiǎn)函數(shù)與相應(yīng)協(xié)變量函數(shù)的乘積，通過描述不同人群在不同時(shí)刻的風(fēng)險(xiǎn)，來探索各危險(xiǎn)因素對生存的影響。其基本形式為：

通過以上公式可以看出，該模型的參數(shù)估計(jì)不依賴于基準(zhǔn)風(fēng)險(xiǎn)函數(shù)的分布類型，即對于不同個(gè)體X1和X2，其風(fēng)險(xiǎn)比

與基準(zhǔn)風(fēng)險(xiǎn)函數(shù)無關(guān)，且不隨時(shí)間t發(fā)生變化。這就是Cox回歸模型最基本的比例風(fēng)險(xiǎn)（proportional hazards，PH）假定。此外，Cox回歸模型還要求滿足對數(shù)線性假定，即協(xié)變量與對數(shù)風(fēng)險(xiǎn)函數(shù)間呈線性［3］。

Cox回歸模型雖然使生存數(shù)據(jù)中的多因素分析成為可能，然而，由于它依賴于嚴(yán)格的假定條件，若資料無法滿足，則會(huì)較大程度地影響結(jié)果的解讀，甚至導(dǎo)致錯(cuò)誤的結(jié)論。因此，在統(tǒng)計(jì)分析前，對PH假定的檢驗(yàn)是重要且必要的。目前Cox回歸模型存在一定濫用現(xiàn)象，多數(shù)研究者在使用該方法時(shí)忽略了對PH假定的檢驗(yàn)，影響了結(jié)果的真實(shí)性和可靠性。本文希望能夠?qū)ΜF(xiàn)有的PH假定的檢驗(yàn)方法進(jìn)行歸納總結(jié)，以提示讀者合理使用Cox回歸模型，選擇恰當(dāng)?shù)姆椒z驗(yàn)數(shù)據(jù)的適用性，建立穩(wěn)定有效的模型。

比例風(fēng)險(xiǎn)（PH）假定的檢驗(yàn)方法

目前，檢驗(yàn)Cox回歸模型PH假定的方法主要有圖示法［4-5］和假設(shè)檢驗(yàn)法［6］兩種。圖示法包括：（1）Cox & K-M比較法，（2）累積風(fēng)險(xiǎn)函數(shù)法，（3）Schoenfeld殘差圖法；假設(shè)檢驗(yàn)法包括：（1）時(shí)協(xié)變量法，（2）線性相關(guān)檢驗(yàn)法，（3）加權(quán)殘差Score法；（4） Omnibus檢驗(yàn)法。本文將對上述方法逐一做簡要介紹并附優(yōu)劣分析。

1.圖示法

（1）Cox & K-M比較法

即比較Cox回歸模型與其他非參數(shù)方法（如K-M曲線）估計(jì)的生存曲線的形態(tài)差異，若二者趨勢基本一致，且無交叉，則認(rèn)為數(shù)據(jù)滿足PH假定。該方法最早由Cox［2］本人提出，后經(jīng)Harrel和Lee［7］等人拓展，可用于計(jì)數(shù)資料、等級資料和計(jì)量資料的分析，其中計(jì)量資料需先將變量離散化，再比較各組的Cox和K-M曲線結(jié)果。

（2）累積風(fēng)險(xiǎn)函數(shù)法

即觀察：①累積風(fēng)險(xiǎn)函數(shù)H（t|X）對于時(shí)間t的趨勢圖，若比例恒定；②累積風(fēng)險(xiǎn)函數(shù)H（t| X）對于基準(zhǔn)累積風(fēng)險(xiǎn)函數(shù)H0（t）的趨勢圖，若斜率恒定；③對數(shù)累積風(fēng)險(xiǎn)函數(shù)lnH（t|X）對于時(shí)間t的趨勢圖，若相互平行；④對數(shù)累積風(fēng)險(xiǎn)函數(shù)差異lnH（t|X1）-lnH（t|X2）對于時(shí)間t的趨勢圖，若恒為常數(shù)。以上四者均可認(rèn)為數(shù)據(jù)滿足PH假定。

對于Cox回歸模型，下列公式等價(jià)：

因此當(dāng)PH假定成立時(shí)，累積風(fēng)險(xiǎn)函數(shù)圖中不同組間的比值H（t|X1）/ H（t| X2）=exp（（X1-X2）β）［11］、對數(shù)累積風(fēng)險(xiǎn)函數(shù)圖中不同組間的差異lnH（t|X1）-lnH（t|X2）=（X1-X2）β［12］應(yīng)為常數(shù)，即顯示為互成比例或平行［13-15］。

此外，對于二值變量（X =0，1），Anderson［16］等提出可以采用H1（t）對H0（t）圖，即H-H圖來檢驗(yàn)PH假定。由于H1（t）=H（t|X =1）=H0（t）eβ，因此H-H圖應(yīng)為一條直線，截距為0，斜率為eβ。此方法同樣適用于連續(xù)變量和等級變量［17-18］。

上述方法雖然實(shí)現(xiàn)了以一條曲線代替兩條曲線來評價(jià)數(shù)據(jù)的適用性，但值得注意的是，在不滿足PH假定的情況下，＾β并不能靠lnH1（t）/ H0（t）進(jìn)行簡單估計(jì)，因?yàn)榇藭r(shí)lnH1（t）/ H0（t）≠lnh1（t）/ h0（t）。故當(dāng)該值隨時(shí)間波動(dòng)較大時(shí)，只能得出數(shù)據(jù)不服從PH假定的結(jié)論，不能以該值的形狀推斷β［5］。

（3）Schoenfeld殘差圖法

基于上述方法均與時(shí)間變量相關(guān)，Schoenfeld［19］定義了一個(gè)不依賴于t的偏殘差，以檢驗(yàn)Cox回歸模型的PH假定。

令Ri為ti時(shí)刻（即第i個(gè)個(gè)體發(fā)生事件時(shí)）的風(fēng)險(xiǎn)集，則該個(gè)體的Schoenfeld殘差可表示為ri=XI-E （Xi|Ri），其中

為給定風(fēng)險(xiǎn)集下的條件期望［20］，可由最大偏似然估計(jì)（maximum partial likelihood estimate）得到的代入計(jì)算。由于為的解［2］，可以證明，當(dāng)PH假定滿足時(shí)故ri對ti圖應(yīng)在0附近波動(dòng)。

此后，Grambsch和Therneau［21］對Scheonfeld殘差進(jìn)行了標(biāo)度調(diào)整，得到加權(quán)Scheonfeld殘差rid·Sβ，其中d為資料中的所有事件數(shù)。ri*的尺度與＾β一致，即ri*對ti圖在＾β附近波動(dòng)，說明數(shù)據(jù)滿足PH假定。同時(shí)，上述兩位學(xué)者在另一篇文獻(xiàn)［22］中提到Score殘差，作為Scheonfeld殘差的擴(kuò)展，還可用于更復(fù)雜的情況。

然而，Scheonfeld殘差圖中的散點(diǎn)趨勢難以評價(jià)，尤其對于二值變量，圖中僅有兩條近乎平行的趨勢線r1=1-E（1|R1）和r0=-E（0|R0），無法了解其波動(dòng)情況［5］。LOWESS（locally-weighted scatterplot smoothing）平滑法［23-24］及其置信區(qū)間的估計(jì)［21］可能會(huì)對殘差圖的可讀性起到較大的幫助，并且將實(shí)際的PH檢驗(yàn)結(jié)果從隨機(jī)趨勢中分離出來。

2.假設(shè)檢驗(yàn)法

（1）時(shí)協(xié)變量（time-by-covariate interactions）法

Cox在建立回歸模型時(shí)，曾提出一個(gè)構(gòu)造時(shí)協(xié)變量的方法來檢驗(yàn)PH假定［2］，即在模型中加入一個(gè)交互作用項(xiàng)X·f（t），其中f（t）為時(shí)間函數(shù)，使模型擴(kuò)展為：

此時(shí)，

對PH假定的檢驗(yàn)即可轉(zhuǎn)化為對γ=0的檢驗(yàn)。

此外，該方法對時(shí)間函數(shù)的選擇亦是十分關(guān)鍵的。線性函數(shù)（f（t）=a + bt）、對數(shù)函數(shù)（f（t）=a + blnt）、指數(shù)函數(shù)（f（t）=a + bexpt）等單調(diào)函數(shù)在以往的文獻(xiàn)中較為常見［2，10，26］，因?yàn)楣潭ǖ膮?shù)a、b便于計(jì)算，且隨時(shí)間單調(diào)增/減的交互項(xiàng)也更易于解讀。對非單調(diào)時(shí)間函數(shù)的研究相對少見，Stablein等［27］曾考慮使用二次函數(shù)（f（t）=at + bt2）進(jìn)行建模，Hess［5］提出也可采用分段線性函數(shù)擬合時(shí)間交互項(xiàng)。20世紀(jì)90年代初期，非參數(shù)模型逐漸被提出并使用，一些研究使用樣條函數(shù)估計(jì)時(shí)協(xié)變量［28-30］，無需確定時(shí)間函數(shù)的具體形式，可以避免模型的選擇錯(cuò)誤導(dǎo)致的結(jié)果偏差，提高檢驗(yàn)效率。但如何選擇合適的節(jié)點(diǎn)位置和數(shù)量將成為新的問題。

（2）線性相關(guān)檢驗(yàn)法

該方法起源于Schoenfeld［19］提出的偏殘差概念，后經(jīng)Harrel和Lee［7］改進(jìn)，逐漸發(fā)展為一種簡便的PH假定檢驗(yàn)方法。其原理主要是：Schoenfeld殘差不依賴于時(shí)間變量，因此Schoenfeld殘差與生存時(shí)間的秩次線性無關(guān)。此時(shí)只需檢驗(yàn)相關(guān)系數(shù)ρ=0（檢驗(yàn)統(tǒng)計(jì)量，即可證明數(shù)據(jù)滿足PH假定。

線性相關(guān)檢驗(yàn)法也同樣適用于其他形式的殘差，如鞅殘差（martingale residual）［22，31］等。定義Mi=δi-H0（ti）exp（Xβ）為第i個(gè)個(gè)體的鞅殘差，則Mi可認(rèn)為是該個(gè)體觀察到的事件數(shù)與Cox回歸模型下期望的事件數(shù)之差，即資料中存在但未被預(yù)測到的部分。在PH假定成立的條件下，鞅殘差同樣具有與時(shí)間t無關(guān)的特性，因此可用于相關(guān)性檢驗(yàn)；此外，鞅殘差還可用于對數(shù)線性假定的識別［32］。

（3）加權(quán)殘差Score法

此方法類似于上述兩種方法的綜合。令β（t）=β+ f（t）γ為時(shí)間相關(guān)的參數(shù)，即

則可以證明［21］，上述方程的Schoenfeld殘差，其期望約為Eri≈Sβ（ti）f（ti）γ，對β（t）=β的檢驗(yàn)等價(jià)于對Schoenfeld殘差的廣義最小二乘檢驗(yàn)。

其中

由此可以構(gòu)造一個(gè)自由度為p的漸進(jìn)χ2統(tǒng)計(jì)量

若能夠得出β（t）與時(shí)間無關(guān)，則表明數(shù)據(jù)滿足PH假定。

（4）Omnibus檢驗(yàn)法

一些研究者［7，14］認(rèn)為，檢驗(yàn)Cox回歸模型的PH假定，還可以通過分割生存時(shí)間，即預(yù)先將時(shí)間按順序劃分為幾個(gè)互不相交的區(qū)間，在每個(gè)區(qū)間內(nèi)分別擬合Cox模型，以比較不同區(qū)間內(nèi)的回歸系數(shù)是否一致。

假設(shè)將時(shí)間劃分為q個(gè)區(qū)間（τ0，τ1），…，（τq-1，τq），其中τ0=0，τq=∞，令

則分段的Cox回歸模型可表達(dá)為

該方法最初由Moreau等［33-34］提出，后經(jīng)O′Quigley和Pessione［35］拓展，成為時(shí)協(xié)變量法的一個(gè)特例［15，23］。此法要求每一區(qū)間內(nèi)的事件數(shù)均衡可比，且不能太少。針對如何確定區(qū)間分割點(diǎn)的位置和數(shù)量，曾被廣泛探討，一些作者［25，35］認(rèn)為選擇事件發(fā)生時(shí)間的分位點(diǎn)是合適的，但該考慮受到刪失機(jī)制的影響，若生存數(shù)據(jù)不滿足隨機(jī)刪失，則難以給出令人信服的結(jié)論；同時(shí)區(qū)間的數(shù)量也依賴于樣本量的大小和生存時(shí)間的分布。有作者提出，可以采用事后分析來確定分割點(diǎn)［5］，但這樣的做法必然會(huì)導(dǎo)致統(tǒng)計(jì)效率的下降，影響結(jié)果的可信度。

討 論

比例風(fēng)險(xiǎn)假定作為Cox回歸模型非常重要的應(yīng)用條件之一，應(yīng)成為模型建立的基本前提。目前，國際上已有一些研究在進(jìn)行生存分析時(shí)提到了PH假定的檢驗(yàn)［36-38］，然而，絕大多數(shù)使用Cox回歸模型的文章中，仍然缺乏對該假定的考慮。截至2015年7月，在Pubmed中搜索“Cox regression”，共析出文獻(xiàn)24167篇，但其中同時(shí)提到“proportional hazards assumption”的文獻(xiàn)僅有29篇，占0.12%；即便放寬搜索條件，在提到“Cox”的121342篇文獻(xiàn)中，也僅有413篇（0.34%）提到了“assumption”，可見在使用Cox回歸模型時(shí)，大多數(shù)作者對其假定的檢驗(yàn)并不關(guān)心。另一方面，在29篇提及比例風(fēng)險(xiǎn)假定的文章中，探討方法學(xué)的文章達(dá)16篇之多，僅13篇臨床研究在建立Cox回歸模型時(shí)考慮了PH假定的適用性，其中8篇得到不符合PH假定的結(jié)論，并依此對模型進(jìn)行了調(diào)整，4篇認(rèn)為變量符合PH假定，還有1篇對PH假定進(jìn)行了提及，但未進(jìn)行檢驗(yàn)。研究者對于PH假定的重視程度還有待提高，目前尚無公認(rèn)的檢驗(yàn)方法，也給PH假定的推廣造成了困難。

本文對目前常見的PH假定檢驗(yàn)方法進(jìn)行了總結(jié)，其中不乏一些簡單直觀的方法，便于使用。本文希望通過描述，讓更多人對PH假定有更深入的理解，使之日后納入Cox回歸模型的評價(jià)成為可能。

本文將檢驗(yàn)PH假定的方法分為圖示法和假設(shè)檢驗(yàn)法兩種。其中圖示法具有清晰、實(shí)用等特點(diǎn)，無需復(fù)雜的計(jì)算，當(dāng)發(fā)現(xiàn)數(shù)據(jù)可能不滿足PH假定時(shí)，還可以指導(dǎo)恰當(dāng)?shù)臋z驗(yàn)統(tǒng)計(jì)量的選擇。圖示法除了上文中提到的三種外，尚有UCP（updated covariate percentile）圖［25］、Aalen圖［39-41］、Arjas圖［42］等，由于計(jì)算過程相對繁瑣，且無相應(yīng)的軟件模塊可以直接繪制，故在此不做贅述。

相比于圖示法，很大程度上依賴于研究者的主觀判斷，假設(shè)檢驗(yàn)法可以更加客觀準(zhǔn)確地給出令人信服的結(jié)論。本文探討了四種較為常見的假設(shè)檢驗(yàn)方法，Nicholas［6］的模擬研究表明，時(shí)協(xié)變量法、線性相關(guān)檢驗(yàn)法和加權(quán)殘差Score法在檢驗(yàn)PH假定時(shí)均有較高的準(zhǔn)確率；其他可行的檢驗(yàn)方法還包括：廣義矩檢驗(yàn)［43］、線性秩檢驗(yàn)［44］、Score檢驗(yàn)［22，45］。假設(shè)檢驗(yàn)法雖然更為正式可靠，但檢驗(yàn)結(jié)果受到樣本量的影響，當(dāng)樣本量較小時(shí)，檢驗(yàn)的靈敏度降低，樣本量過大時(shí)，又可能因?yàn)楦怕试蚓芙^原本成立的假設(shè)。因此圖示法在已有的研究中更受歡迎，一般只需判斷既定模型近似滿足PH假定即可。各方法的優(yōu)劣分析見表1。

由表1可知，各方法特點(diǎn)不同，適應(yīng)于不同條件下的應(yīng)用。在探索性分析階段，建議采用Cox & K-M比較法或累積風(fēng)險(xiǎn)函數(shù)法進(jìn)行簡單的比例風(fēng)險(xiǎn)描述，若符合PH假定，即可直接建模；當(dāng)發(fā)現(xiàn)數(shù)據(jù)偏離PH假定時(shí)，再考慮計(jì)算Schoenfeld殘差進(jìn)行繪圖或檢驗(yàn)，以便為改進(jìn)Cox回歸模型提供幫助。

表1　比例風(fēng)險(xiǎn)檢驗(yàn)方法的優(yōu)劣分析

參考文獻(xiàn)

［1］Langova K.Survival analysis for clinical studies.Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub，2008，152：303-307.

［2］Cox DR.Regression models and life-tables（with discussion）.Journal of the Royal Statistical Society，Series B：Methodology，1972，34：187-220.

［3］Anderson PK.Statistical models based on counting processes.Berlin：Springer-Verlag，1993.

［4］Hess KR.Graphical methods for assessing violations of the proportional hazards assumption in Cox regression.Statistics in medicine，1995，14：1707-1723.

［5］余紅梅，何大衛(wèi).檢查Cox模型比例風(fēng)險(xiǎn)假定的幾種圖示法.中國衛(wèi)生統(tǒng)計(jì)，2000，17：215-218.

［6］Nicholas HN.An empirical comparison of statistical tests for assessing the proportional hazards assumption of Cox′s model.Statistics in Medicine，1997，16：611-626.

［7］Harrel FE，Lee KL.Verifying assumptions of the Cox proportional hazards mode1.Proceedings of the Eleventh Annual SAS Users Group International Conference，1986：823-828.

［8］Breslow NE.Covariance analysis of censored survival data.Biometrics，1974，30：88-89.

［9］Link CL.Confidence intervals for the survival function using Cox′s proportional-hazard model with covariates.Biometrics，1984，40：601-609.

［10］Kalbfleisch JD，Prentice RL.The Statistical analysis of failure time data.New York：Wiley，1980.

［11］Muenz LR.Comparing survival distributions：a review for nonstatisticians II.Cancer Investigation，1983，1：537-545.

［12］Schumacher M.Two-sample tests of Cramer-von Mises-and Kolmogorov-Smirnov-type for randomly censored data.International Statistical Review，1984，52：263-281.

［13］Lagakos SW.The graphical evaluation of explanatory variables in proportional hazard regression models.Biometrika，1981，68：93-98.

［14］Kay R.Goodness of fit methods for the proportional hazards regression model：a review.Revue Epidemiol Sante Publique，1984，32：185-198.

［15］Gray RJ.Some diagnostic methods for Cox regression models through hazard smoothing.Biometrics，1990，46：93-102.

［16］Andersen PK.Testing goodness of fit of Cox′s regression and life model.Biornetrics，1982，38：67-77.

［17］Wei LJ.Testing goodness of fit for proportional hazards model with censored observations.Journal of the American Statistical Association，1984，79：649-652.

［18］Gill RD，Schumacher M.A simple test of the proportional hazards assumption.Biometrika，1987，74：289-300.

［19］Schoenfeld D.Partial residuals for the proportional hazards regression model.Biornetrika，1982，69：239-241.

［20］Cox DR.Partial likehood.Biometrika，1975，62：269-276.

［21］Grambsch PM，Therneau TM.Proportional hazard tests and diagnostics based on weighted residuals.Biometrika，1994，81：515-526.

［22］Therneau TM，Grambsch PM.Martingale-based residuals for survival models.Biomatrika，1990，77：147-160.

［23］Pettitt AN，Daud IB.Investigating time dependence in Cox′s proportional hazards model.Applied Statistics，1990，39：313-329.

［24］Cleveland WS.Robust locally weighted regression and smoothing scatterplots.Journal of the American Statistical Association，1979，74：829-836.

［25］Thaler HT.Nonparametric estimation of the hazard ratio.Journal of the American Statistical Association，1984，79：290-293.

［26］Anderson JA，Senthilselvan A.A two-step regression model for hazard functions.Applied Statistics，1982，31：44-51.

［27］Stablein DM，Carter WH，Novak JW.Analysis of survival data with nonproportional hazard functions.Controlled Clinical Trials，1981，2：149-159.

［28］Gray RJ.Flexible methods for analyzing survival data using splines，with applications to breast cancer prognosis.Journal of the American Statistical Association，1992，87：942-951.

［29］Hess KR.Assessing time-by-covariate interactions in proportional hazards regression model using cubic spline functions.Statistics in Medicine，1994，13：1045-1062.

［30］Heinzl H，Kaider A，Zlabriger G.Assessing interactions binary timedependent covariates with time in Cox proportional hazards regression models using cubic spline functions.Statistics in medicine.1996，15：2589-2601.

［31］Kay R.Proportional hazards regression models and the analysis of censored survival data.Applied Statistics，1977，26：227-237.

［32］Marzec L，Marzec P.Generalized martingale-residual processes for goodness-of-fit inference in Cox′s type regression models.The Annals of Statistics，1997，25：683-714.

［33］Moreau T，O′Quigley J，Mesbah M.A global goodness-of-fit statistic for the proportional hazards model.Applied Statistics，1985，34：212-218.

［34］Moreau T，O′Quigley J，Lellouch J.On D Schoenfeld′s approach for testing the proportional hazards assumption.Biometrics，1986，73：513-515.

［35］O′Quigley J，Pessione F.Score tests for homogeneity of regression effect in the proportional hazards model.Biometrics，1989，45：135-144.

［36］Santori G，F(xiàn)ontana I，Bertocchi M，et al.Application and validation of Cox regression models in a single-center series of double kidney transplantation.Transplantation Proceedings.2010，42：1098-1103.

［37］Abadi A，Saadat S，Yavari P，et al.Comparison of Aalen′s additive and Cox proportional hazards models for breast cancer survival：analysis of population-based data from British Columbia，Canada.Asian Pac J Cancer Prev，2011，12：3113-3116.

［38］Abadi A，Yavari P，Dehghani-Arani M，et al.Cox models survival analysis based on breast cancer treatments.Iran J Cancer Prev，2014，7：124-129.

［39］Aalen OO.Interaction between life history events.Nonparametric analysis for prospective and retrospective data in the presence of censoring.Scandinavian Journal of Statistics，1980，7：161-171.

［40］Mau J.On a graphical method for the detection of time-dependent effects of covariates in survival data.Applied Statistics，1986，35：245-255.

［41］Henderson R，Milner A.Aalen plots under proportional hazards.Applied Statistics，1991，40：401-409.

［42］Arjas E.A graphical method for assessing goodness of fit in Cox′s proportional hazards mdoel.Journal of the American Statistical Association，1988，83：204-212.

［43］Horowitz JL，Neulnann GR.A generalized moments specification test of the proportional hazards model.Journal of the American Statistical Association，1992，87：234-240.

［44］Chappell R.A note on linear rank test and Gill and Schumacher′s tests of proportionality.Biometrika，1992，79：l99-201.

［45］Lin DY，Wei LJ，Ying Z.Checking the Cox model with cumulative sums of martingale-based residuals.Biometrika，1993，80：557-572.

（責(zé)任編輯：郭海強(qiáng)）

通信作者：△李衛(wèi)，E-mail：liwei@ mrbc-nccd.com