張?zhí)鹛?,劉紅偉 ,劉媛媛 ,李長平 ,2,胡良平
(1.天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院,天津 300070;2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會,北京 100029;3.軍事科學(xué)院研究生院,北京 100850
本期科研方法專題的第三篇文章介紹了生存資料參數(shù)回歸模型有關(guān)的基礎(chǔ)知識,包括構(gòu)建三個(gè)常見的生存資料參數(shù)回歸模型的基本原理、基于圖示法判斷生存時(shí)間服從何種概率分布的方法、基于最大似然估計(jì)法求解參數(shù)回歸模型中的參數(shù)和兩個(gè)參數(shù)回歸模型擬合優(yōu)度的比較。本文通過一個(gè)實(shí)例,利用SAS軟件的LIFEREG過程介紹生存資料參數(shù)回歸模型的SAS實(shí)現(xiàn)方法。
【例1】本文所用數(shù)據(jù)是對149例糖尿病患者17年追蹤調(diào)查數(shù)據(jù),包括生存結(jié)局、生存時(shí)間、隨訪開始時(shí)年齡、體重指數(shù)等[1]。由于存在刪失數(shù)據(jù),此資料為生存資料。Diabetic數(shù)據(jù)集中生存時(shí)間(t,年)、隨訪開始時(shí)年齡(age1,歲)、體重指數(shù)(BMI)、診斷出糖尿病時(shí)的年齡(age0,歲)、收縮壓(SBP,mmHg)和舒張壓(DBP,mmHg)是定量自變量;吸煙狀況(smk,0表示不吸煙,1表示曾吸煙,2表示吸煙)、心電圖讀數(shù)(ECG,1表示正常,2表示臨界,3表示異常)、是否有冠心?。–HD,0表示無,1表示有)和結(jié)局(status,0表示截尾,1表示死亡)都是定性變量。利用以下SAS數(shù)據(jù)步程序,創(chuàng)建名為Diabetic的數(shù)據(jù)集:
【SAS程序說明】因篇幅所限,在CARDS語句后的數(shù)據(jù)僅列出前5個(gè)和后5個(gè)觀測。因?yàn)槭湛s壓和舒張壓有一定關(guān)聯(lián)[2],分析時(shí)取加權(quán)平均血壓(MBP),即令 MBP=SBP*(1/3)+DBP*(2/3),權(quán)重系數(shù)分別為1/3和2/3。
對數(shù)生存圖、對數(shù)-對數(shù)生存圖和對數(shù)失效比生存圖分別見圖1、圖2和圖3。
圖1 對數(shù)生存圖
圖2 對數(shù)-對數(shù)生存圖
圖3 對數(shù)失效比生存圖
圖1中的線圖有些彎曲,圖2在整體上呈線性趨勢,與圖3接近。圖示法表明,基于現(xiàn)有數(shù)據(jù),該生存資料同時(shí)適合Weibull分布模型和Log-logistic分布模型,本文選擇Weibull分布模型進(jìn)行擬合(圖示法雖然簡單直觀,但不夠精確,若深入探討該數(shù)據(jù)是否更加適合Weibull分布模型,可參考其他方法[3])。
因?yàn)閃eibull分布回歸模型嵌套于廣義Gamma分布回歸模型[4],為選擇最優(yōu)的模型,本文將分別擬合Weibull分布回歸模型和廣義Gamma分布回歸模型,根據(jù)似然比檢驗(yàn)結(jié)果來確定最優(yōu)模型。利用以下SAS過程步程序,構(gòu)建Weibull分布回歸模型和廣義Gamma分布回歸模型。
【SAS程序說明】采用LIFEREG過程分別擬合Weibull分布回歸模型和廣義Gamma分布回歸模型,class語句列出離散型自變量,model語句左邊為生存時(shí)間和生存結(jié)局變量(括號內(nèi)為截尾值的標(biāo)記),右邊為預(yù)測變量(即自變量),包括離散和連續(xù)自變量。
Weibull分布回歸模型輸出結(jié)果:
廣義Gamma分布回歸模型輸出結(jié)果:
【SAS結(jié)果說明】Fit Statistics表是基于t為響應(yīng)變量的最大似然估計(jì)得到的統(tǒng)計(jì)量,可以用來比較不同協(xié)變量的模型。Analysis of Maximum Likelihood Parameter Estimates表給出了參數(shù)的估計(jì),包括預(yù)測變量的回歸系數(shù)和參數(shù)分布中參數(shù)的估計(jì)值。其中“Scale”代表“尺度參數(shù)”,“Shape”代表“形狀參數(shù)”。
因Weibull分布回歸模型包含2個(gè)參數(shù),廣義Gamma分布回歸模型包含3個(gè)參數(shù),則似然比擬合優(yōu)度檢驗(yàn)[5]的自由度為 1,χ2統(tǒng)計(jì)量為 8.249(=85.899-77.650)(注意:;顯然,8.249>6.635),P<0.01,即兩個(gè)分布擬合效果之間差異有統(tǒng)計(jì)學(xué)意義,故采用-2logL值最小的分布,即廣義Gamma分布(-2logL值為77.650)。廣義Gamma分布回歸模型輸出結(jié)果中,只有age1、MBP、CHD、smk和ECG五個(gè)自變量有統(tǒng)計(jì)學(xué)意義。故僅保留這五個(gè)自變量,重新采用廣義Gamma分布回歸模型來進(jìn)行分析。結(jié)果如下:
廣義Gamma分布回歸模型的分析結(jié)果表明,隨訪開始時(shí)年齡的回歸系數(shù)為負(fù)值(-0.0408),說明隨訪年齡越大,生存時(shí)間越短,死亡風(fēng)險(xiǎn)越大;smk變量0水平與2水平比較時(shí),對應(yīng)的P值小于0.05,其回歸系數(shù)為正值(0.4048),說明不吸煙者生存時(shí)間長于吸煙者;ECG變量1水平與3水平比較時(shí),對應(yīng)的P值小于0.001,其回歸系數(shù)為正值(1.0879),說明心電圖正常者生存時(shí)間長于心電圖異常者。
在現(xiàn)有的統(tǒng)計(jì)軟件中,進(jìn)行生存資料參數(shù)回歸模型建模時(shí)尚不能篩選自變量,只能依據(jù)計(jì)算結(jié)果,采用手工方法刪除沒有統(tǒng)計(jì)學(xué)意義的自變量,這在一定程度上影響了最優(yōu)回歸模型的產(chǎn)生。當(dāng)自變量中含有兩個(gè)以上定量自變量時(shí),根據(jù)統(tǒng)計(jì)分析的經(jīng)驗(yàn),盡可能產(chǎn)生出一些派生自變量(例如平方項(xiàng)、交叉乘積項(xiàng)等)參與構(gòu)建回歸模型,可能有助于找到擬合優(yōu)度更好的回歸模型。生存資料參數(shù)回歸模型的建模策略與logistic回歸模型建模策略類似,對最后的結(jié)果而言,保留在回歸模型中的自變量,既要考慮其應(yīng)具有統(tǒng)計(jì)學(xué)意義也要考慮其實(shí)際意義,即回歸系數(shù)正負(fù)號的含義在專業(yè)上是可以得到合理解釋的。
本文通過一個(gè)實(shí)例,比較詳細(xì)地介紹了基于SAS軟件實(shí)現(xiàn)參數(shù)回歸模型的擬合方法;通過圖示法大致判斷生存時(shí)間可能服從何種分布類型;最后,還針對所擬合的兩個(gè)參數(shù)回歸模型,進(jìn)行了擬合優(yōu)度檢驗(yàn),從而可以確定最優(yōu)回歸模型。需注意的是,當(dāng)基于圖示法確定生存時(shí)間符合某種參數(shù)分布且其對應(yīng)的模型屬于嵌套模型(特指其參數(shù)取某些特定值時(shí),原先復(fù)雜的模型就退化成一個(gè)相對簡單的模型,例如,當(dāng)威布爾分布模型中的形狀參數(shù)γ=1時(shí),它就退化成指數(shù)分布模型了)時(shí),則需要采用似然比檢驗(yàn)確定擬合效果最優(yōu)的模型;否則,需要根據(jù)所擬合的兩個(gè)模型各自的參數(shù)數(shù)目和其他擬合統(tǒng)計(jì)量的數(shù)值,綜合考慮后再選定相對較優(yōu)的模型。