半?yún)?shù)與非參數(shù)可加風(fēng)險模型對達菲和中藥治療H1N1流感效果的比較

曹志強王 楊李 衛(wèi)

半?yún)?shù)與非參數(shù)可加風(fēng)險模型對達菲和中藥治療H1N1流感效果的比較

曹志強1，2王 楊1△李 衛(wèi)1

目的介紹半?yún)?shù)與非參數(shù)可加風(fēng)險模型在事件發(fā)生時間分析（time-to-event analysis）方面的應(yīng)用，并比較它們結(jié)果的異同。方法利用鞅方法得到半?yún)?shù)可加風(fēng)險模型的參數(shù)估計，運用核密度估計非參數(shù)可加風(fēng)險模型，用R軟件進行相應(yīng)的計算及統(tǒng)計分析。結(jié)果半?yún)?shù)可加風(fēng)險模型中的系數(shù)相當(dāng)于非參數(shù)可加風(fēng)險模型可變系數(shù)的平均。相對于對照組，達菲組、中藥組、達菲加中藥組分別治療53、56、29個發(fā)熱病人就能多一個發(fā)熱消退的病例。兩個模型的生存概率圖都表明，達菲、中藥、達菲加中藥相對不吃藥都能有效治療H1N1流感，中藥與達菲的療效相當(dāng)，中西藥結(jié)合的療效最好。結(jié)論在研究者需要對風(fēng)險的絕對變化進行刻畫或Cox比例風(fēng)險假定不滿足時，可加風(fēng)險模型或許是一個合適的替代。

半?yún)?shù)模型 非參數(shù)模型 H1N1流感

事件發(fā)生時間分析（time-to-event analysis）常用來研究風(fēng)險因子與臨床事件發(fā)生的聯(lián)系，可選擇的模型包括參數(shù)、半?yún)?shù)以及非參數(shù)模型。在參數(shù)模型中，常用的是加速失效模型，半?yún)?shù)模型主要是Cox模型，非參數(shù)模型中用的最多的是KM估計。

一個很有名但應(yīng)用還不廣的用來分析風(fēng)險函數(shù)與協(xié)變量關(guān)系的是非參數(shù)可加風(fēng)險模型［1］，1989年由Aalen提出，基本形式如下：

其中α0（t）是基線函數(shù)，Z（t）＝［Z1（t），…，Zp（t）］T是協(xié)變量矩陣，αj（t）是第j個協(xié)變量的回歸系數(shù)，對模型（1）進一步的探討見文獻［2］。受模型（1）的啟發(fā)，D.Y.Lin和Z.Ying在1994年提出了下面的半?yún)?shù)可加風(fēng)險模型［3］：

其中β0（t）是基線函數(shù)，βj為不隨時間變化的回歸系數(shù)，Zj（t）是第j個協(xié)變量在t時刻的值，模型（2）進一步的擴展見文獻［4］。不管是半?yún)?shù)還是非參數(shù)可加風(fēng)險模型，回歸系數(shù)度量的都是絕對風(fēng)險，這點與Cox模型中的系數(shù)度量相對風(fēng)險不同。因此，如果研究者感興趣的是風(fēng)險的絕對變化或Cox比例風(fēng)險假定不滿足時，可加風(fēng)險模型或許是一個合適的替代。

資料與方法

1.資料

一項以治療H1N1流感為目的的研究［5］，比較奧司他韋（達菲）和傳統(tǒng)中藥湯（麻杏石甘湯和銀翹散加減方）的治療效果，410例確診為輕癥H1N1流感的成年患者被隨機非盲分成4組：對照組、達菲組、中藥組、達菲加中藥組。目標(biāo)變量time為從入組治療到結(jié)束的時間；status指發(fā)熱是否消退；age是患者的年齡；g2、g3、g4是三個啞變量，分別指患者服用的是達菲、中藥湯、達菲加中藥；fb48h是發(fā)病至入組時間是否大于48小時的二值變量；s2、s3、s4是三個中心啞變量。為了解決結(jié)點問題，將time每個值加上［0，1］之間的隨機數(shù)。

2.方法

考慮半?yún)?shù)和非參數(shù)可加風(fēng)險模型，分析達菲和傳統(tǒng)中藥治療H1N1流感的效果。對于半?yún)?shù)可加風(fēng)險模型，即模型（2），當(dāng)協(xié)變量與時間獨立時，其形式與Cox模型類似。它們都有一個未知的基線函數(shù)，回歸系數(shù)均為常數(shù)，不過Cox模型是定義在乘積尺度上，而模型（2）是定義在加性尺度上。利用計數(shù)過程中的鞅方法，借鑒Cox模型估計的偏極大似然法，可以得到模型（2）中回歸系數(shù)和標(biāo)準(zhǔn)差的估計，估計的構(gòu)造及大樣本性質(zhì)見文獻［3］。

對模型（2）擬合效果的檢驗，可以定義類似于Cox模型中的Cox-Snell殘差。比如，在觀測時刻ti，假設(shè)協(xié)變量為事件發(fā)生的示性函數(shù)為δi，那么殘差ri的估計為，其中若模型（2）正確，那么ri是來自標(biāo)準(zhǔn)指數(shù)分布的樣本［6］。因此，我們用Nelson-Aalen估計，計算ti時刻的累積風(fēng)險函數(shù)（ti），如果（ri，（ti））近似為一條斜率為1的直線，那么模型（2）正確。

模型（2）假定協(xié)變量對風(fēng)險函數(shù)的影響是常數(shù)，實際上協(xié)變量的影響也可能與時間相關(guān)。從形式上來看，模型（2）是模型（1）的特殊情況，但它們參數(shù)估計的方法完全不同。在模型（1）中，直接估計回歸系數(shù)是困難的，因而轉(zhuǎn)向估計與其等價的累積回歸系數(shù)，第j個協(xié)變量的累積回歸系數(shù)定義如下：

采用非參數(shù)最小二乘法，將累積回歸系數(shù)矩陣A（t）＝（A0（t），A1（t），…，Ap（t））T估計出來，估計的大樣本性質(zhì)見文獻［1］。為了得到回歸系數(shù)αj（t）的估計，需要選擇合適的核函數(shù)K（x）以及窗寬b，并且利用下面的核密度估計公式：

至于模型（1）擬合效果的檢驗，也可以定義類似的Cox-Snell殘差。R軟件中的timereg、ahaz包有專門分析非參數(shù)和半?yún)?shù)可加風(fēng)險模型的程序，本文所有的統(tǒng)計分析結(jié)果在R 3.0.1中得到。

結(jié) 果

雖然Cox模型是分析該類問題的首選模型，但是對Cox模型的協(xié)變量進行檢驗時發(fā)現(xiàn)，變量g3（P＝0.006）、s3（P＝0.015）并不滿足比例風(fēng)險假定。文獻［7］指出，當(dāng)Cox模型的假設(shè)不滿足時，可以嘗試用可加風(fēng)險模型進行分析。

表1是用半?yún)?shù)與非參數(shù)可加風(fēng)險模型擬合檢驗的結(jié)果。在半?yún)?shù)模型中，除了age、fb48h、s3，其他變量都有統(tǒng)計學(xué)意義。控制協(xié)變量的前提下，平均來說，相對于不吃藥，服用達菲的治療為1／0.0189≈53人，說明用達菲每治療53例H1N1流感病人就能多一個發(fā)熱消退的病例。類似的，相對不吃藥，服用中藥湯每治療56（1／0.0177）個病人就會多一個發(fā)熱消退的病例，服用達菲加中藥每治療29（1／0.0342）個病人就能多一個發(fā)熱消退的病例。通過Wald檢驗發(fā)現(xiàn)，達菲與中藥湯系數(shù)之間的差別不顯著，達菲加中藥相對達菲或中藥湯的系數(shù)差別均顯著。從表1不難得知，相對對照組，三個試驗組對治療H1N1流感均有效，達菲的藥效與中藥湯相當(dāng)，但達菲加中藥比單純達菲或中藥湯的藥效都要好。

表1 半?yún)?shù)與非參數(shù)可加風(fēng)險模型的擬合檢驗

非參數(shù)可加風(fēng)險模型的系數(shù)隨著時間變化，這里是用TST統(tǒng)計量及相應(yīng)的P值檢驗變量有無統(tǒng)計學(xué)意義，TST的構(gòu)造以及大樣本性質(zhì)見文獻［1］。從表1還可以得到，在變量有無統(tǒng)計學(xué)意義方面，兩個模型的結(jié)果是一樣的。圖1至圖3是非參數(shù)可加風(fēng)險模型中三組藥的累積系數(shù)圖以及95%的置信區(qū)間，直線是該變量在半?yún)?shù)可加風(fēng)險模型中的累積回歸系數(shù)。

圖1 達菲的累積系數(shù)圖及置信區(qū)間

圖2 中藥湯的累積系數(shù)圖及置信區(qū)間

圖3 達菲加中藥的累積系數(shù)圖及置信區(qū)間

從圖1至圖3可以直觀地看出三組藥在各個時段對風(fēng)險函數(shù)的累積影響。當(dāng)協(xié)變量的累積系數(shù)過了某個時刻近似平穩(wěn)，則該變量不再對風(fēng)險函數(shù)有影響，即該風(fēng)險因子對治療H1N1流感不再有效。從三個圖還可以看出，半?yún)?shù)可加風(fēng)險模型中三組藥的累積系數(shù)在非參數(shù)可加風(fēng)險模型中的95%置信區(qū)間內(nèi)，實際上，其他幾個變量也是如此。文獻［8］指出，如果半?yún)?shù)可加風(fēng)險模型中協(xié)變量的累積系數(shù)在非參數(shù)可加風(fēng)險模型95%置信區(qū)間內(nèi)，可以將非參數(shù)可加風(fēng)險模型中的可變系數(shù)視為固定。

由累積回歸系數(shù)，利用核密度估計可以得到回歸系數(shù)。在核密度估計中，文獻［6，9］表明核函數(shù)對估計結(jié)果影響不大，本文選文獻中廣泛應(yīng)用的Epanechnikov核，窗寬的選擇是用有偏交叉驗證法確定。非參數(shù)可加風(fēng)險模型累積系數(shù)的核密度估計見圖4。

圖4 非參數(shù)可加風(fēng)險模型的核密度估計

圖4的核密度估計展示了非參數(shù)模型中系數(shù)隨時間變化的情況，其直觀地說明了在控制協(xié)變量的前提下，相對不吃藥，三組藥在各時刻增加風(fēng)險的大小。在前40小時，達菲、達菲加中藥前的系數(shù)隨時間變化的趨勢類似，不過在前28小時達菲加中藥的系數(shù)大一些。中藥湯的系數(shù)與它們很不相同，尤其是12～20小時期間，達菲、達菲加中藥的系數(shù)處于較低水平，中藥湯的卻處于較高水平。40小時后，三組藥的系數(shù)波動很大，可能是由于樣本量減少的緣故。從圖4的縱坐標(biāo)和表1中的結(jié)果發(fā)現(xiàn)，半?yún)?shù)模型中的系數(shù)是非參數(shù)模型中可變系數(shù)的平均。

圖5 半?yún)?shù)可加風(fēng)險模型的生存概率圖

圖6 非參數(shù)可加風(fēng)險模型的生存概率圖

圖5和圖6是半?yún)?shù)和非參數(shù)可加風(fēng)險模型在其他變量的調(diào)整下，對照組、達菲組、中藥組、達菲加中藥組的生存概率圖。兩模型中三組藥的療效都比不吃藥的好，達菲加中藥的療效比單純達菲或中藥湯的都要好。在半?yún)?shù)模型中，中藥湯與達菲的療效幾乎沒有差別。而非參數(shù)模型中，在前30小時中藥湯的療效稍好于達菲，30小時后達菲的療效稍好于中藥湯。

圖7 達菲加中藥在三個模型中的生存曲線

圖7是達菲加中藥在KM估計、半?yún)?shù)和非參數(shù)可加風(fēng)險模型中的生存曲線。從圖中可以看出，半?yún)?shù)和非參數(shù)可加風(fēng)險模型的生存曲線與KM估計的差不多。文獻［7］指出可以利用模型估計的生存曲線與KM估計的相比，曲線與KM估計的越相近，表明模型擬合的越好。因此，圖7顯示不論是半?yún)?shù)還是非參數(shù)可加風(fēng)險模型，都較好地擬合了這組數(shù)據(jù)。

圖8和圖9是半?yún)?shù)與非參數(shù)可加風(fēng)險模型的Cox-Snell殘差圖。相對來說，半?yún)?shù)可加風(fēng)險模型的擬合程度比非參數(shù)可加風(fēng)險模型要好。

圖8 非參數(shù)可加風(fēng)險模型的cox-snell殘差圖

圖9 半?yún)?shù)可加風(fēng)險模型的cox-snell殘差圖

討 論

本文用半?yún)?shù)和非參數(shù)可加風(fēng)險模型研究了達菲和中藥治療H1N1流感的效果，并從參數(shù)估計、變量有無統(tǒng)計學(xué)意義、在其他變量調(diào)整下各實驗組的生存概率、模型擬合程度等方面進行比較。在本文的例子中，兩個模型中各變量有無統(tǒng)計學(xué)意義是類似的，半?yún)?shù)模型中的系數(shù)相當(dāng)于非參數(shù)模型中累積系數(shù)核密度估計后的平均。相對于不吃藥，達菲、中藥湯、達菲加中藥分別治療53、56、29個發(fā)熱病人就能多一個發(fā)熱消退的病例。生存概率圖則表示相對不吃藥，三組藥都能有效治療甲流，達菲與中藥湯的療效相當(dāng)，達菲加中藥湯的療效最好。Cox-Snell殘差圖表明，半?yún)?shù)可加風(fēng)險模型比非參數(shù)模型的整體擬合效果要好。

正如預(yù)期的一樣，不管是半?yún)?shù)還是非參數(shù)可加風(fēng)險模型，其回歸系數(shù)估計值都比Cox模型中的小很多。這不奇怪，因為可加風(fēng)險模型與Cox模型針對不同的問題。Cox模型估計的相對風(fēng)險在理解事件的關(guān)聯(lián)程度方面非常有用，某種藥對疾病的絕對風(fēng)險影響可能較小，但相對風(fēng)險仍然可以很大。如果感興趣的是因治療而實際獲益的病例數(shù)目時，本文方法則更有優(yōu)勢，因為絕對風(fēng)險在公共衛(wèi)生規(guī)劃和干預(yù)方面發(fā)揮著重要作用［10］。

若感興趣的是協(xié)變量的平均效果，可以用半?yún)?shù)可加風(fēng)險模型。如果想檢驗協(xié)變量相對于對照組的影響是否隨時間變化，非參數(shù)可加風(fēng)險模型是一個不錯的選擇。更多的文獻則建議對于同一組數(shù)據(jù)，既要用Cox模型也需要用可加風(fēng)險模型擬合，因為它們反映了風(fēng)險函數(shù)與協(xié)變量之間兩種不同的關(guān)系。文獻［11］指出，即使協(xié)變量滿足比例風(fēng)險，Cox模型有時也會得出誤導(dǎo)性的結(jié)果。

1.Aalen OO.A linear regression model for the analysis of life times.Statistics in medicine，1989，8（8）：907-925.

2.Aalen OO.Further results on the non-parametric linear regression model in survival analysis.Statistics in medicine，1993，12（17）：1569-1588.

3.Lin DY，Ying Z.Semiparametric analysis of the additive risk model.Biometrika，1994，81（1）：61-71.

4.Lin DY，Ying Z.Semiparametric analysis of general additive-multiplicative hazard models for counting processes.The Annals of Statistics，1995，23（5）：1712-1734.

5.Wang C，Cao B，Liu QQ，et al.Oseltamivir compared with the Chinese traditional therapy maxingshigan-yinqiaosan in the treatment of H1N1 influenza：a randomized trial.Annals of Internal Medicine，2011，155（4）：217-25.

6.Klein，Moeschberger JM.Survival Analysis：Techniques for Censored and Truncated Data.Springer，1997.

7.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 Pacific Journal of Cancer Prevention，2011，12：3113-3116.

8.Mckeague IW，Sasieni PD.A partly parametric additive risk model.Biometrika，1994，81（3）：501-514.

9.Aydemir，Biller.Kernel smoothing of Aalen′s linear regression model. Sonderforschungsbereich，386，1997，101.

10.Xie X，Strickier HD，Xue X.Additive Hazard Regression Models：An Application to the Natural History of Human Papillomavirus.Computational and mathematical methods in medicine，2013，796270.

11.Amato DA.A generalized Kaplan-Meier estimator for heterogeneous populations Communications in Statistics-Theory and Methods，1988，17：263-286.

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

Comparison of Semi-parametric and Nonparametric Additive Hazards Models：with Tamiflu and Traditional Chinese Medicine to Treat H1N1 Flu Effect for Example

Cao Zhiqiang，Wang Yang，Li Wei（Chinese Academy of Medical Sciences，Beijing Union Medical College，National Center for Cardiovascular Diseases，F(xiàn)u Wai Hospital，State Key Laboratory of Cardiovascular Disease（100037），Beijing）

ObjectiveTo introduce semi-parametric and nonparametric additive hazards models and its application in time-to-event analysis，and compare similarities and differences of the results.MethodsUsing the method of martingale to get semi-parametric additive risk model parameter estimation，using kernel density technique to estimate the nonparametric additive risk model，using R software to do corresponding calculation and statistical analysis.ResultsSemi-parametric additive hazard model coefficient is equal to the nonparametric additive hazard model of variable coefficient of the average.Relative to the control group，Tamiflu，traditional Chinese medicine，tamiflu plus traditional Chinese medicine treat53，56，29 patients respectively can be one more example of a fever subsided.The survival probability plots in two models show that，relatively to control group，tamiflu，traditional Chinese medicine，tamiflu plus traditional Chinese medicine all of three are effective in treating H1N1 influenza.Traditional Chinese medicine′s efficacy is similar to tamiflu，the efficacy of Chinese and Western medicine combined is the best.ConclusionWhen the absolute change in risk is of primary interest or when the proportional hazard assumption for the Cox proportional hazard model is violated，additive hazard model may be more appropriate.

Semi-parametric model；Nonparametric model；H1N1 influenza

1.中國醫(yī)學(xué)科學(xué)院，北京協(xié)和醫(yī)學(xué)院，國家心血管病中心，阜外心血管病醫(yī)院，心血管疾病國家重點實驗室（100037）

2.北京師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院（100875）

△通信作者：王楊，Email：wangyang＠m(xù)rbc-nccd.com