呂晶晶 陳金寶 吳天敏 賀鈺錕 駱雨桐 劉 丹 侯雅文 陳 征△
固定點處k組生存率比較的非參數(shù)統(tǒng)計方法研究*
呂晶晶1陳金寶1吳天敏1賀鈺錕1駱雨桐1劉 丹1侯雅文2陳 征1△
目的 當僅需比較某固定時間點上組間生存率的差異,或者不滿足組間風險率成比例假設時,如生存曲線存在交叉,log-rank檢驗不再適用,且現(xiàn)有固定點檢驗法僅限于兩組間,故本文發(fā)展固定點處多組間生存率的比較方法。方法 首先提出多組間五種固定點檢驗法(經(jīng)典法、對數(shù)轉(zhuǎn)換法、雙對數(shù)轉(zhuǎn)換法、反正弦平方根轉(zhuǎn)換法及邏輯轉(zhuǎn)換法),并通過Monte Carlo模擬評價五種方法在不同情形下的一類錯誤和檢驗效能。最后對滿足和不滿足風險率成比例假設兩個實例用上述方法進行分析。結(jié)果 綜合Monte Carlo模擬得到的一類錯誤及檢驗效能結(jié)果,以經(jīng)典法和反正弦平方根轉(zhuǎn)換法最為激進,對數(shù)轉(zhuǎn)換法略保守,邏輯轉(zhuǎn)換法最為保守,而雙對數(shù)轉(zhuǎn)換法最為穩(wěn)健。結(jié)論 在進行多組間生存率比較時,當僅想比較多組間某固定點處生存率差異或者組間不滿足成比例假設,可使用上述五種固定點檢驗法,其中優(yōu)先建議使用雙對數(shù)轉(zhuǎn)換法。
固定點 多組間比較 生存率 Kaplan-Meier估計 Monte Carlo模擬
臨床隨訪研究中,多組間生存率差異的比較是重要研究內(nèi)容之一,當組間生存曲線交叉時,通常結(jié)果顯示多組在整個隨訪期上是無統(tǒng)計學意義的,Liu[1]和Li[2]指出由于交叉而不滿足風險率成比例假設,導致log-rank檢驗的檢驗效能降低,此時會掩飾有部分時間段或者某個有重要意義的時間點上有差異的事實。另外,相對于整個隨訪期內(nèi)生存率的比較,有時候研究者感興趣的是某個時間點上生存率的差異,如在第5年多組療效的差異。例如在一項探究三類種族對淋巴細胞白血病預后影響的研究中,生存曲線大約在第5、10年出現(xiàn)交叉(圖1)。經(jīng)Grambsch-Therneau檢驗[3]發(fā)現(xiàn)三組間風險率不成比例(χ2=7.030,P=0.008),因此log-rank檢驗發(fā)現(xiàn)無統(tǒng)計學差異(χ2=3.500,P=0.177)的結(jié)果不可靠。除此之外,多條生存曲線存在交叉在臨床研究中頻繁出現(xiàn),如Oosten[4]和Kawada[5]分別指出Llor等[6]和Benito-León等[7]的研究中由于三條生存曲線存在交叉,log-rank檢驗或Cox回歸模型并不適用。另外,Iacobelli等[8]指出生存曲線發(fā)生交叉時比較整體上的差異并不是最重要的,可能更感興趣的是某部分區(qū)域或者某些固定時間點上組間的差異問題。本文僅針對某些固定時間點上組間的差異進行探討研究。
陳金寶等[9]介紹了基于對生存率進行不同形式轉(zhuǎn)換,構(gòu)造兩組間固定點上生存率差異的多種檢驗法,但卻無法直接比較多組(3組及以上)的情況。Klein等[10-11]雖然提到多組固定點生存率比較的問題,卻沒有嘗試多種方法的轉(zhuǎn)換校正。因此,本文將針對多組間固定點處生存率差異的比較及其多種轉(zhuǎn)換方法展開研究。
假設在第k組(k=1,2,…,p)中,nk為該組樣本量,tki表示該組第i個個體事件發(fā)生的時間,其中tk1<tk2<…<tkn,dki和Yki分別表示該組tki上的事件數(shù)和風險人數(shù)。Kaplan-Meier估計及其方差分別為^Sk(t)=
圖1 不同種族下淋巴細胞白血病患者的生存曲線
在兩條生存曲線某固定點的比較中[10],原假設為在某固定點t上,兩組生存率相等,即S1(t)=S2(t)。設對第k組生存率的轉(zhuǎn)換函數(shù)為φ,則生存率轉(zhuǎn)換后形式為φ(^Sk(t)),對應方差為V[φ(^Sk(t))],檢驗法基本構(gòu)造形式為:
在原假設成立的前提下,檢驗統(tǒng)計量Z服從自由度為1的χ2分布。為了將其推廣到適應于多組固定時刻點上生存率的差異檢驗,此時的原假設H0為在某固定時間點t,多組生存率都相等,即S1(t)=S2(t)=…=Sp(t)。首先將公式(1)的分子部分擴充為一個含(p-1)個元素的向量,記為A,表示如下:
然后將公式(1)的分母部分擴充為一個(p-1)×(p-1)的方差協(xié)方差矩陣,記為∑。其中,∑的主對角線上元素為^V[φ(^Sk(t))]+^V[φ(^Sp(t))],k=1,2,…,p-1。非對角線上元素為^V[φ(^Sp(t))]。由此得到多組固定點上比較的統(tǒng)計量基本形式為:
在H0成立的前提下,檢驗統(tǒng)計量X(公式(3))服從自由度為(p-1)的χ2分布。
由于生存率不服從正態(tài)分布,并基于兩組比較的研究[9-10],本文提出5種多組間固定點上生存率比較的檢驗方法,均構(gòu)造不同的向量A(公式(2))以及對應的方差協(xié)方差陣∑,主對角線元素和非主對角線元素構(gòu)造不同,最后代入公式(3)得到最終檢驗統(tǒng)計量X。
1.經(jīng)典法(naive)
經(jīng)典法中,向量A(公式(2))里的φ變換為φ(^Sk(t))=^Sk(t),則
∑矩陣主對角線上元素為^V[^Sk(t)]+^V[^Sp(t)],非對角線上元素為^V[^Sp(t)]。在H0成立的前提下,最終統(tǒng)計量服從自由度為(p-1)的χ2分布,下列檢驗法也全部服從自由度為(p-1)的χ2分布。
2.對數(shù)轉(zhuǎn)換法(log)
3.雙對數(shù)轉(zhuǎn)換法(cloglog)
4.反正弦平方根轉(zhuǎn)換方法(arcsin)
為了評價五種多組固定點檢驗法的性能,采用Monte-Carlo模擬研究各檢驗法的檢驗效能和一類錯誤,其中在檢驗效能方面,三組的生存時間均由指數(shù)分布產(chǎn)生(Exp(λk)),為控制生存曲線開口大小,選擇不同參數(shù)(λ1=0.693,λ2=0.462,λ3=0.315)使得時間點為1時的生存率分別約為0.50,0.63,0.73(如圖2),即相對于第一組而言,第二組HR=1.5,第三組HR=2.2;刪失時間C均由指數(shù)分布產(chǎn)生。在一類錯誤方面,三組的生存時間均由參數(shù)為0.2的指數(shù)分布(Exp(0.2))產(chǎn)生,刪失時間C由服從于U(0,a)、U(0,b)和U(0,c)的均勻分布產(chǎn)生。記錄時間為t=min(T,C),δ=1[T≤C],通過改變刪失時間分布參數(shù),可使得每組的平均刪失率相同且約為0、15%、30%、50%。通過預模擬發(fā)現(xiàn),合并樣本時間25%分位數(shù)近似為2,進而比較三組在時刻點為2時的生存情況??紤]樣本均衡(n1,n2,n3均為30、60或100)和不均衡(n1=30,n2=n3=60;n1=n2=30,n3=60;n1=30,n2=60,n3=100)的情形,每一種參數(shù)組合下模擬10000次,顯著水平α=0.05。
圖2 檢驗效能模擬情形
表1展示的是一類錯誤模擬結(jié)果。naive轉(zhuǎn)換法相比其他檢驗有更高的一類錯誤,arcsin轉(zhuǎn)換法次之,且都高于檢驗水準0.05,易得出有差異的結(jié)果,其余三種檢驗法的一類錯誤均較小,顯得稍微保守。隨著樣本量增加,五種檢驗法一類錯誤均越接近0.05,但在樣本量不均衡時,log轉(zhuǎn)換法的一類錯誤相對地出現(xiàn)較大上升,特別是在樣本量為(30,60,100)時。
表2是五種方法檢驗效能的模擬結(jié)果。隨樣本量的增加和刪失率的減少,五種檢驗法檢驗效能均呈上升趨勢。naive轉(zhuǎn)換法和arcsin轉(zhuǎn)換法有較高的檢驗效能,cloglog轉(zhuǎn)換法的檢驗效能比其余兩種檢驗法的高。
綜合一類錯誤和檢驗效能,以naive法和arcsin轉(zhuǎn)換法最易得有統(tǒng)計學意義的結(jié)論,其余三種轉(zhuǎn)換法顯得相對保守,其中l(wèi)og轉(zhuǎn)換法一類錯誤高于其余兩組,但是檢驗效能卻是最低的,顯得最不易發(fā)現(xiàn)差異的存在,cloglog轉(zhuǎn)換法顯得最為穩(wěn)健。
表1 五種檢驗方法一類錯誤模擬結(jié)果
表2 五種檢驗方法檢驗效能模擬結(jié)果
續(xù)表2
本文分別提供滿足和不滿足風險率成比例假設的兩個實例分析,用于進行多組生存曲線的固定點檢驗法的驗證。
【例1】一項關于探究不同種族對淋巴細胞白血病預后影響的研究,即圖1對應的實例。研究起點為初診確認為淋巴細胞白血病,終點事件是患者發(fā)生死亡,其余為右刪失。3類人種的平均生存時間分別約為5年、7年和8年,刪失率分別約為25%、23%和26%。三組的log-rank檢驗結(jié)果(χ2=3.500,P=0.177)并不可靠,進一步對第3、5、10、15、20年進行固定點上生存率檢驗,結(jié)果顯示除第5年和第10年無統(tǒng)計學差異外,其余時間點上均有。通過模擬發(fā)現(xiàn),cloglog轉(zhuǎn)換法所得到的結(jié)果最為穩(wěn)健,故固定點檢驗結(jié)果以cloglog轉(zhuǎn)換法為準。進一步選取有統(tǒng)計學差異的固定時間點,利用Bonferroni法進行多重比較,發(fā)現(xiàn)在第3年上黑種人和美國印第安人患者的預后情況(P=0.012)、第15年(P=0.033)和第20年(P=0.033)上黑種人和亞洲或太平洋島民患者的預后情況具有統(tǒng)計學差異。
表3 例1的分析結(jié)果
【例2】一項關于三種療法對白血病患者預后影響的研究,共納入66名白血病患者,按所接受的不同療法分為3個組別(圖3):療法1含25人,療法2含19人,療法3含22人,刪失率分別約為5%、10%和11%,中位生存時間分別約為29天、75天和40天。三組滿足風險率成比例假設(χ2=0.349,P=0.555),log-rank檢驗(χ2=5.600,P=0.061)顯示不同療法下的白血病患者生存率無統(tǒng)計學差異。但由圖3觀察發(fā)現(xiàn)三條曲線中后期開口較大,因此進一步對第25、50、75、100、150、200天進行固定點上生存率比較檢驗,發(fā)現(xiàn)在第75天、100天和150天存在差異。進一步選取有統(tǒng)計學差異的固定點,利用Bonferroni法進行多重比較,發(fā)現(xiàn)在第75天(P=0.050)、第100天(P=0.044)和第150天(P=0.034)上療法1和療法2下患者的預后情況具有統(tǒng)計學差異。
圖3 不同療法下白血病患者的生存曲線
表4 例2的分析結(jié)果
多組間生存率比較是臨床隨訪研究中最重要的研究內(nèi)容之一,其中l(wèi)og-rank檢驗是整體差異檢驗的經(jīng)典檢驗方法之一。但在進行多組間生存率比較時,僅對生存曲線上某固定點處的生存率差異感興趣時,可使用本文所提出的固定點處多組間生存率比較的檢驗方法。同時,由本文模擬檢驗結(jié)果得到,在五種多組間的固定點檢驗方法中,以naive轉(zhuǎn)換法和arcsin轉(zhuǎn)換法較激進,log轉(zhuǎn)換法和logit轉(zhuǎn)換法較保守,以cloglog轉(zhuǎn)換法最為穩(wěn)健,建議使用cloglog轉(zhuǎn)換法。
[1]Liu K,Qiu P,Sheng J.Comparing two crossing hazard rates by Cox proportional hazardsmodelling.Stat Med,2007,26(2):375-391.
[2]Li H,Han D,Hou Y,et al.Statistical inference methods for two crossing survival curves:a comparison ofmethods.PLoSOne,2015,10(1):e116774.
[3]Grambsch PM,Therneau TM.Proportional hazards tests and diagnostics based on weighted residuals.Biometrika,1994,81(3):515-526.
[4]van Oosten DC.Re:Efficacy of anti-inflammatory or antibiotic treatment in patients with non-complicated acute bronchitis and discoloured sputum:randomised placebo controlled trial.BMJ,2014.(http://www.bm j.com/content/347/bm j.f5762/rr/680601).
[5]Kawada T.Long sleep duration in elders without dementia increases risk of dementia mortality(NEDICES).Neurology,2015,85(4):388.
[6]Llor C,Moragas A,Bayona C,et al.Efficacy of anti-inflammatory or antibiotic treatment in patientswith non-complicated acute bronchitis and discoloured sputum:randomised placebo controlled trial.BMJ,2013,347:f5762.
[7]Benito-Leon J,Louis ED,Villarejo-Galende A,et al.Long sleep duration in elderswithout dementia increases risk of dementiamortality(NEDICES).Neurology,2014,83(17):1530-1537.
[8]Iacobelli S,EBMT Statistical Committee.Suggestions on the use of statisticalmethodologies in studies of the European Group for Blood and Marrow Transplantation.Bone Marrow transplantation,2013,48(1),S1-37.
[9]陳金寶,邱李斌,王北琪,等.固定點處組間生存率比較的統(tǒng)計檢驗法.中華流行病學雜志,2015,36(2):186-188.
[10]Klein JP,Logan B,Harhoff M,et al.Analyzing survival curves at a fixed point in time.Stat Med,2007,26(24):4505-4519.
[11]Klein JP,Moeschberger ML,et al.Survival Analysis:Techniques for cencored and truncated data.2th ed.New York:Springer,2003.234-238.
(責任編輯:鄧 妍)
M ethods of Com paring M ultiple Survival Rates at a Fixed Time Point
Lv Jingjing,Chen Jinbao,Wu Tianm in,et al(Department of Biostatistics,School of Public Health,Southern Medical University(510515),Guangzhou)
Objective In comparing multiple survival curves at a fixed pointin time,log-rank test is inapplicable.Besides,its power would be worse in crossing survival curves because of notmeet the proportional hazard assumption.Hence,we use themethod of comparing survival rates at fixed point to dealw ith them.However,thismethod can only use for two groups.In view of the above,we considered the comparison of multiple survival curves at fixed point.M ethods We first proposed 5 methods to comparemultiple survival curves at fixed point(naive,log,cloglog,arcsin,logit).Monte Carlo simulationswere carried out to evaluate the type Ierror and power of thesemethods.Finally,we used two examples for analysis by using abovemethods.Results Comprehensive results of type Ierror and power,naive and arcsin were themost radical ways;log and logitwere more conservative;and cloglog was themost robust.Conclusion In comparison ofmultiple survival rates,someone can choose our methods of comparingmultiple survival curves at fixed-pointwhen these survival curves do notmeet the proportional hazard assumption or only interested at fixed-point in time.And we suggested cloglogmethod.
Fixed point;Comparison ofmultiple survival curves;Kaplan-Meier estimation;Monte Carlo simulations
國家自然科學基金(81673268),廣東省自然科學基金(2017-1714050008015)
1.南方醫(yī)科大學公共衛(wèi)生學院生物統(tǒng)計學系(510515)
2.暨南大學經(jīng)濟學院統(tǒng)計學系
△通信作者:陳征,E-mail:zchen@smu.edu.cn