區(qū)間刪失數(shù)據(jù)的3種統(tǒng)計(jì)模型分析及其SAS實(shí)現(xiàn)

張倩倩， 王純杰， 佟知真， 李純凈

(長春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院， 吉林 長春 130012)

?

區(qū)間刪失數(shù)據(jù)的3種統(tǒng)計(jì)模型分析及其SAS實(shí)現(xiàn)

張倩倩， 王純杰*， 佟知真， 李純凈

(長春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院， 吉林 長春 130012)

借助SAS9.4中PROC ICLIFETEST、PROC ICPHREG過程步編寫宏程序,同步實(shí)現(xiàn)了區(qū)間刪失數(shù)據(jù)的生存函數(shù)估計(jì)、廣義Log-Rank檢驗(yàn)和PH類型回歸模型的統(tǒng)計(jì)推斷。結(jié)合回溯研究中368個樣本HIV-1感染時間的區(qū)間刪失數(shù)據(jù)給出實(shí)證分析。

區(qū)間刪失； ICLIFETEST； 廣義對數(shù)秩檢驗(yàn)； ICPHREG； 宏程序

0 引 言

生存分析是對試驗(yàn)或調(diào)查得到的人或生物的生存時間數(shù)據(jù)進(jìn)行推斷,在醫(yī)學(xué)實(shí)踐中有著廣泛應(yīng)用。一般稱給定事件的出現(xiàn)時間為生存時間[1]，分析生存時間數(shù)據(jù)通常意味著解決3個問題：估計(jì)生存函數(shù)，比較處理組或者生存函數(shù)，評估協(xié)變量的影響或者依靠生存時間的解釋變量。區(qū)間刪失數(shù)據(jù)是生存時間中越來越常見的一種數(shù)據(jù)，在過去幾十年里，出現(xiàn)了許多分析區(qū)間刪失數(shù)據(jù)的統(tǒng)計(jì)方法。Turnbull[2]找到了類似右刪失數(shù)據(jù)下的Kaplan-Meier估計(jì)的自相合算法來獲得生存函數(shù)估計(jì)；王弄升[3]2012年利用SAS軟件中宏程序%EMICM給出區(qū)間刪失數(shù)據(jù)生存函數(shù)的估計(jì)；Sun[4]等把Log-Rank檢驗(yàn)推廣到區(qū)間刪失數(shù)據(jù)中，提出廣義對數(shù)秩檢驗(yàn)；Finkelstein D M[5]給出區(qū)間刪失數(shù)據(jù)的COX回歸模型。但是基于SAS軟件還沒有完整的程序可以同步實(shí)現(xiàn)區(qū)間刪失數(shù)據(jù)3種統(tǒng)計(jì)分析任務(wù)。因此，文中借助SAS9.4中PROC ICLIFETEST[6]、PROC ICPHREG過程步編寫宏程序,實(shí)現(xiàn)了區(qū)間刪失數(shù)據(jù)的生存函數(shù)估計(jì)、廣義Log-Rank檢驗(yàn)和PH比例風(fēng)險類型的回歸模型統(tǒng)計(jì)推斷。

1 模型及方法介紹

1.1 區(qū)間刪失的數(shù)據(jù)類型

設(shè)T為非負(fù)的隨機(jī)變量，代表研究中個體的生存時間，對于區(qū)間刪失數(shù)據(jù)，只能知道T落在某個區(qū)間內(nèi)，即T∈(L,R],在這里L(fēng)≤R。區(qū)間刪失數(shù)據(jù)可分為Ⅰ型區(qū)間刪失與Ⅱ型區(qū)間刪失。Ⅰ型區(qū)間刪失數(shù)據(jù)可以表示{C,δ=I(T≤C)},C代表觀測時間，I(·)是示性函數(shù)。Ⅱ型區(qū)間刪失數(shù)據(jù)是包括有限區(qū)間的區(qū)間刪失數(shù)據(jù)(L,R)，假設(shè)每個個體觀測兩次，L,R是兩個隨機(jī)變量，L≤R以概率1成立。

1.2 區(qū)間刪失數(shù)據(jù)的非參數(shù)估計(jì)

1.2.1 EMICM算法[1]

SAS軟件計(jì)算非參數(shù)生存函數(shù)的方法是EMICM法。該算法是一個把自相合與ICM算法簡單結(jié)合的混合算法。考慮一個包含n個獨(dú)立個體帶有生存函數(shù)為S(t)的失效時間的研究。令Ti代表個體i(i=1,2，…,n)的生存時間。假設(shè)在Ti上的區(qū)間刪失數(shù)據(jù)被觀測如下：

O={(Li,Ri];i=1,2，…,n}

則區(qū)間刪失數(shù)據(jù)的似然函數(shù)為：

NPMLE算法就是要最大化此似然函數(shù)。在自相合算法中對數(shù)似然為：

Wellner,Zhan[7]指出，當(dāng)NPMLE存在且唯一并且對數(shù)似然函數(shù)連續(xù)可微的時候,EMICM算法收斂到NPMLE。應(yīng)用EMICM算法需要選擇一個收斂準(zhǔn)則。這里選擇基于Robertson[8]1988年提出的Fenchel優(yōu)化條件。這種準(zhǔn)則下如果滿足

1.2.2 廣義對數(shù)秩檢驗(yàn)[1]

j=1,2，…,m

類似的對于第l組，l=1,2，…,p+1，且j=1,2，…,m時，在sj處失效數(shù)和在風(fēng)險數(shù)的估計(jì)為：

為了檢驗(yàn)H0,給出統(tǒng)計(jì)量

此檢驗(yàn)統(tǒng)計(jì)量是服從自由度為P的卡方分布。

1.2.3 PH比例風(fēng)險模型[1]

則似然函數(shù)的對數(shù)形式為:

利用極大似然函數(shù)的得分方程估計(jì)β，基準(zhǔn)生存函數(shù)S0可由連續(xù)的階梯函數(shù)來估計(jì)。

2 實(shí)證分析

2.1 區(qū)間刪失數(shù)據(jù)及變量情況

數(shù)據(jù)來源于美國和歐洲的16個研究中心的第5中心[9]，此中心采用回溯研究的方法,主要研究帶有血友病的病人感染HIV-1的風(fēng)險。血友病人的治療血液制品來自于成千上萬個捐贈者的血漿制成的Ⅷ型凝血因子和Ⅳ型因子，所以這些病人都存在感染HIV-1的風(fēng)險。文中對病人的HIV-1感染時間只得到了區(qū)間刪失數(shù)據(jù)，且病人依據(jù)他們每年獲得血液制品的平均劑量被分配到不同的組。來自第5中心的368個病人的觀測數(shù)據(jù)，不考慮進(jìn)入研究的病人HIV-1抗體狀態(tài)，不接受或接受低劑量(1～20 000U)的Ⅷ型凝血因子(兩組病人的組別記為NF和LDF)見表1。

兩組病人的數(shù)量分別為236和13,這些數(shù)據(jù)以季度為單位，0代表研究開始時間(1978年1月1日)。

表1 HIV-1感染數(shù)據(jù)

續(xù)表1

2.2 SAS程序

對于上述血友病患者HIV-1感染的區(qū)間刪失數(shù)據(jù)，首先要分別建立兩組的SAS數(shù)據(jù)集，然后用ICLIFETEST過程步估計(jì)出兩個受試組HIV-1感染的時間生存函數(shù)，并檢驗(yàn)兩組生存曲線的區(qū)別，最后用ICPHREG過程步建立PH比例風(fēng)險回歸模型。區(qū)間刪失數(shù)據(jù)3種生存統(tǒng)計(jì)分析可由SAS宏程序%ICLIFETEST-PHREG同步實(shí)現(xiàn)。

2.2.1 創(chuàng)建NF組的數(shù)據(jù)集

data NF;

input lTime rTime @@;

Stage=0;

datalines;

55. 55. 56. 54. 53. 57.

56. 56. 54. 56. 54. 55.

...;

run;

2.2.2 創(chuàng)建LDF組的數(shù)據(jù)集

data LDF;

input lTime rTime @@;

Stage=1;

datalines;

720 920 025 57.

2326 821 2026 2527

...;

run;

2.2.3 估計(jì)兩組病人的生存概率及繪制生存曲線

data zq;

set LDF NF;

run;

proc iclifetest data=zq plot=s(cl) impute(seed=1234);/*ICLIFETEST過程步對應(yīng)數(shù)據(jù)，畫出帶95%置信帶的生存曲線圖*/

strata stage;/*識別定義分組的變量*/

time (lTime,rTime);/*區(qū)間刪失的左右觀測時間*/

run;

2.2.4 兩組病人數(shù)據(jù)的廣義Log-Rank檢驗(yàn)

proc iclifetest data=zq impute(seed=1234);

time (lTime,rTime);

test stage;/*檢驗(yàn)的組別*/

run;

2.2.5 生存時間與協(xié)變量的PH比例風(fēng)險回歸

proc icphreg data=zq;

class Stage / desc;

model (lTime,rTime) = Stage / basehaz=pch(intervals=(10));/*協(xié)變量Stage與區(qū)間刪失數(shù)據(jù)建立PH比例風(fēng)險模型*/

hazardratio Stage;/*求風(fēng)險比例*/

run;

上述5個程序是分別建立數(shù)據(jù)集，做區(qū)間刪失數(shù)據(jù)的3種統(tǒng)計(jì)模型分析的SAS程序，對于臨床試驗(yàn)得到的這類區(qū)間刪失數(shù)據(jù)，如果利用SAS程序做分析，需要進(jìn)行上述復(fù)雜的5個程序的運(yùn)行，文中給出一個宏程序，同步實(shí)現(xiàn)上述模型分析。

2.2.6 區(qū)間刪失數(shù)據(jù)同步實(shí)現(xiàn)3種統(tǒng)計(jì)任務(wù)的SAS宏程序%ICLIFETEST-PHREG

%macro iclifetest-phreg(data,var1,var2);

proc iclifetest data=&data plot=s(cl) impute(seed=1234);

strata &var1;

time (lTime,rTime);

run;

proc iclifetest data=&data impute(seed=1234);

time(lTime,rTime);

test &var1;

run;

proc icphreg data=&data;

class &var2/ desc;

model (lTime,rTime) = &var2/ basehaz=pch(intervals=(10));

hazardratio &var1;

run;

%mend;

%iclifetest-phreg(zq,stage,stage)

宏程序%ICLIFETEST-PHREG可以對含二分變量和其他協(xié)變量的區(qū)間刪失數(shù)據(jù)同步實(shí)現(xiàn)生存函數(shù)的估計(jì)、廣義的Log-Rank檢驗(yàn)及PH比例風(fēng)險類型的回歸分析。其語法結(jié)構(gòu)是：data：含二分變量和其他協(xié)變量的區(qū)間刪失數(shù)據(jù)集；Var1：某一個二分協(xié)變量；Var2：做回歸分析的所有的協(xié)變量集合。

2.3 結(jié)果分析

兩個不同組的非參數(shù)生存估計(jì)與生存函數(shù)曲線的宏程序結(jié)果分別如圖1和表2所示。

圖1 NF組與LDF組估計(jì)的生存函數(shù)

組別區(qū)間失效概率生存概率標(biāo)準(zhǔn)誤差NF4180.00550.99450.005419210.01950.98050.010922260.07420.92580.018727270.08560.91440.019928310.11500.88500.021132570.12340.87660.0214LDF0140.00001.00000.000015180.13900.86100.036819200.35350.64650.047121230.36480.63520.047224250.43730.56270.047426260.50620.49380.046027290.54560.45440.043930570.56060.43940.0432

圖1中相應(yīng)的NF組的生存概率要大于LDF組的生存概率，即不接受Ⅷ型凝血因子的血友病實(shí)驗(yàn)組的生存概率要高于接受低劑量Ⅷ型凝血因子的實(shí)驗(yàn)組。表2中NF組的失效概率要小于LDF組的失效概率。

兩組生存概率是否相等的檢驗(yàn)結(jié)果及回歸的顯著性檢驗(yàn)結(jié)果分別見表3和表4。

表3 檢驗(yàn)兩組生存概率是否相等

表4 回歸的顯著性檢驗(yàn)

在顯著性水平為0.05的情況下，表3和表4的p值均小于0.000 1，即拒絕原假設(shè)，LDF組和NF組的生存概率有明顯區(qū)別。表4說明該回歸模型參數(shù)是顯著的。分組的風(fēng)險比估計(jì)見表5。

表5 分組的風(fēng)險比估計(jì)

表5結(jié)果顯示，患有血友病的病人接受低劑量的Ⅷ型凝血因子組感染HIV-1病毒的風(fēng)險概率要高于不接受組病人感染HIV-1病毒的風(fēng)險，且是后者的7.360倍。

3 結(jié) 語

基于區(qū)間刪失數(shù)據(jù)生存函數(shù)的估計(jì)算法除了EMICM算法還有自相合、ICM、Turnbull算法等。廣義對數(shù)秩檢驗(yàn)也成為檢驗(yàn)生存函數(shù)是否相等的有用工具。實(shí)例說明，NF與LDF治療組的生存函數(shù)有顯著差異。回歸模型檢驗(yàn)出分組協(xié)變量對生存時間存在有效的影響。文中給出的SAS宏程序可以針對帶有協(xié)變量的區(qū)間刪失數(shù)據(jù)同步實(shí)現(xiàn)生存函數(shù)的估計(jì)、比較處理組及協(xié)變量對生存時間的影響這3種統(tǒng)計(jì)推斷，對臨床工作人員整理、分析實(shí)驗(yàn)結(jié)果有很大幫助。

[1] Sun J. The statistical analysis of interval-censored failure time data[M]. [S.l.]： Springer Science & Business Media,Inc,2006.

[2] Turnbull B W. The empirical distribution function with arbitrarily grouped, censored and truncated data[J]. Journal of the Royal Statistical Society,1976,38(3):290-295.

[3] 王弄升,張晉聽,駱福添.含有區(qū)間刪失數(shù)據(jù)時的生存函數(shù)估計(jì)及其SAS實(shí)現(xiàn)[J].中國醫(yī)院統(tǒng)計(jì),2012,19(1):1-4.

[4] Zhao Q, Sun J. Generalized log-rank test for mixed interval-censored failure time data [J]. Statistics in Medicine,2004,23(10):1621-1629.

[5] Finkelstein D M. A proportional hazards model for interval-censored failure time data [J]. Biometrics,1987,42(4):845-54.

[6] Guo C, Ying S, Gordon J. Analyzing interval-censored data with the ICLIFETEST procedure [J]. SAS Global Forum,2014,279:327-345.

[7] Jon A Wellner, Yihui Zhan. A hybrid algorithm for computation of the nonparametric maximum likelihood estimator from censored data [J]. Journal of the American Statistical Association,1997,439(92):945-959.

[8] Robertson T, Wright F T, Dykstra R L. Order restricted statistical inference [J]. Journal of the American Statistical Association,1990,85(410):111-112.

[9] Kroner B L, Rosenberg P S, Aledort L M, et al. HIV-1 infection incidence among persons with hemophilia in the United States and western Europe,1978-1990. Multicenter Hemophilia Cohort Study[J]. Journal of Acquired Immune Deficiency Syndromes,1994,7(3):279-86.

SAS-based study on three statistical models for interval censored data analysis

ZHANG Qianqian, WANG Chunjie*, TONG Zhizhen, LI Chunjing

(School of Basic Sciences, Changchun University of Technology, Changchun 130012， China)

Applying the procedures such as PROC ICLIFETEST and PROC ICPHREG in SAS9.4, we realize the statistic deduction for the survival function, Generalized log-rank test and PH regression for the interval censored data by means of macro programming. Interval censored data during the infection time from 368 samples are analyzed in the prospectively study.

interval censored; ICLIFETEST; generalized Log-rank test; ICPHREG; macro program.

2017-01-22

國家自然科學(xué)基金(青年基金)資助項(xiàng)目(11301037); 國家自然科學(xué)基金資助項(xiàng)目(11571051，11671054); 吉林省教育廳“十三五”規(guī)劃項(xiàng)目(2016317)

張倩倩(1991-)，女，漢族，山東德州人，長春工業(yè)大學(xué)碩士研究生，主要從事生存分析方向研究,E-mail:zhangqqxiao@163.com.*通訊作者：王純杰(1978-)，女，漢族，遼寧燈塔人，長春工業(yè)大學(xué)副教授，博士，主要從事生存分析方向研究,E-mail:cjwang2014@126.com.

10.15923/j.cnki.cn22-1382/t.2017.2.02

O 212.3

A

1674-1374(2017)02-0111-06