駱煥園 劉慧芳 吳艷喬 何成奇 楊 珉,4△

·應(yīng)用研究·

1.四川大學(xué)華西公共衛(wèi)生學(xué)院(610041)

2.四川省醫(yī)學(xué)科學(xué)院·四川省人民醫(yī)院康復(fù)醫(yī)學(xué)科

3.四川大學(xué)華西醫(yī)院康復(fù)醫(yī)學(xué)科

4.英國(guó)諾丁漢大學(xué)醫(yī)學(xué)院

#駱煥園，劉慧芳為共同第一作者

△通信作者：楊珉，E-mail:yangmin2013@scu.edu.cn；何成奇，E-mail：hxkfhcq@126.com

根據(jù)重復(fù)觀察數(shù)據(jù)分析臨床病人轉(zhuǎn)歸的統(tǒng)計(jì)學(xué)問(wèn)題和實(shí)踐

駱煥園1#劉慧芳2,3#吳艷喬1何成奇3△楊 珉1,4△

目的探討疾病轉(zhuǎn)歸中重復(fù)測(cè)量數(shù)據(jù)存在的數(shù)據(jù)、統(tǒng)計(jì)分析問(wèn)題及策略。方法以蘆山震區(qū)傷員疾病轉(zhuǎn)歸影響因素研究為例,以Barthel指數(shù)為病人結(jié)局測(cè)量值,采用多水平多模型嘗試法、直尺逼近法和生存分析模型等解決疾病臨床變化、轉(zhuǎn)折平臺(tái)期確定及影響因素等臨床研究問(wèn)題。結(jié)果獲得了克服天花板效應(yīng)的不同傷員Barthel指數(shù)隨時(shí)間變化的多水平模型,確定了疾病轉(zhuǎn)歸事件的時(shí)點(diǎn),探討了疾病轉(zhuǎn)歸的影響因素。結(jié)論分析疾病轉(zhuǎn)歸的重復(fù)測(cè)量數(shù)據(jù)有一系列的統(tǒng)計(jì)問(wèn)題和分析策略。需進(jìn)一步發(fā)展處理天花板效應(yīng)的隨機(jī)效應(yīng)模型。

疾病轉(zhuǎn)歸 重復(fù)測(cè)量數(shù)據(jù) 多水平模型 統(tǒng)計(jì)學(xué)問(wèn)題 天花板效應(yīng)

疾病的轉(zhuǎn)歸需通過(guò)重要指標(biāo)的變化情況來(lái)表征,從而了解病情變化,調(diào)整治療方案[1-3]。比如在比較兩種治療老年股骨粗隆間骨折的方法療效時(shí),需分別在術(shù)后隨訪12～30個(gè)月,觀察Harris髖關(guān)節(jié)功能評(píng)分等[4]。這種重復(fù)測(cè)量數(shù)據(jù)具有以下特征:(1)初次測(cè)量時(shí)間、病情及基線資料存在差異;(2)多次測(cè)量時(shí)間間隔和次數(shù)不一致;(3)多次測(cè)量數(shù)據(jù)間存在相關(guān)性。目前,多水平或隨機(jī)系數(shù)模型是處理這類數(shù)據(jù)的有效方法[5-8]。而在數(shù)據(jù)處理及分析過(guò)程中,測(cè)量指標(biāo)的“天花板”效應(yīng)也需依據(jù)實(shí)際情況予以考慮[9],但目前醫(yī)學(xué)研究中少有探討如何處理此問(wèn)題。

本文通過(guò)臨床實(shí)例,探討不平衡重復(fù)測(cè)量資料中涉及的數(shù)據(jù)問(wèn)題(如天花板效應(yīng))和統(tǒng)計(jì)方法問(wèn)題(如選擇時(shí)間變量的形式)等,為今后相關(guān)研究提供參考。

資 料

資料來(lái)源于四川某教學(xué)醫(yī)院一項(xiàng)地震傷員康復(fù)干預(yù)研究項(xiàng)目。研究對(duì)象為四川蘆山震區(qū)306例傷員,主要研究目的為分析康復(fù)干預(yù)在地震傷員轉(zhuǎn)歸中的作用及影響轉(zhuǎn)歸的因素。以Barthel指數(shù)(Barthel index,BI)作為反映傷員轉(zhuǎn)歸情況的主要指標(biāo),年齡、性別、主診斷類別、受傷被困時(shí)間、受傷至接受康復(fù)治療時(shí)間、康復(fù)干預(yù)項(xiàng)目數(shù)、總住院時(shí)間為協(xié)變量。其中,BI為目前臨床常用的一種評(píng)定日常生活活動(dòng)能力的指標(biāo),以獨(dú)立完成10項(xiàng)任務(wù)(如進(jìn)食等)程度來(lái)反映傷員身體狀況[10],最低0分,最高100分。傷員在不同時(shí)間點(diǎn)入院,均接受基線BI評(píng)定,根據(jù)評(píng)分情況開始干預(yù)后間隔不等時(shí)間評(píng)定其BI直至出院前的末次評(píng)定。每位傷員評(píng)定次數(shù)從1次到11次不等。入院時(shí)評(píng)分高的傷員干預(yù)時(shí)間一般較短,干預(yù)評(píng)估次數(shù)較少,反之亦多。

本例關(guān)心三個(gè)問(wèn)題:第一,BI隨時(shí)間變化規(guī)律如何？據(jù)此推測(cè)BI變化是否有平臺(tái)期？若有則第二,BI到達(dá)平臺(tái)期需多少時(shí)間？哪些協(xié)變量對(duì)此有影響？第三,可否預(yù)測(cè)不同傷員所需時(shí)間？本文介紹我們回答這些問(wèn)題時(shí)在統(tǒng)計(jì)方法選擇上的一些嘗試和思考。

圖1 306例傷員個(gè)體的BI隨時(shí)間變化趨勢(shì)圖

傷員數(shù)據(jù)特征見(jiàn)表1,多數(shù)傷員測(cè)量次數(shù)在1～6之間,測(cè)量次數(shù)少的傷員BI基線值較高、病情較輕,提示測(cè)量次數(shù)因傷員轉(zhuǎn)歸情況而異。隨時(shí)間推移(測(cè)量次數(shù)增加),每次測(cè)量的BI均數(shù)逐漸增加,增加速度不等,傷員大致趨于康復(fù)。由圖1可見(jiàn),BI大致呈拋物線狀隨時(shí)間推移而增加,增加速度先快后慢,個(gè)體間基線值及增加速度有很大差異。部分個(gè)體BI在入院后很短時(shí)間內(nèi)快速上升并持續(xù)停留在最大值100分處。

統(tǒng)計(jì)分析方法

1.BI變化規(guī)律分析

BI為不平衡重復(fù)測(cè)量資料,具有兩個(gè)層次結(jié)構(gòu)(兩個(gè)水平),水平1是重復(fù)測(cè)量值,用i指示;水平2是個(gè)體,用j指示。由于多水平模型能夠接受模型中的截距及多個(gè)協(xié)變量的系數(shù)(斜率)估計(jì)可為隨機(jī)變量的事實(shí),并設(shè)定特定統(tǒng)計(jì)量來(lái)估計(jì)截距及斜率估計(jì)值的隨機(jī)效應(yīng)及分布[5],故采用兩水平模型來(lái)反映不同傷員的BI有不同的基線值和不同的變化速度的事實(shí)。

以tij表示BI測(cè)量時(shí)間,首次測(cè)量時(shí)記為第0天,則模型一般形式為

yij=f(tij)+uj(tij)+e0ij

(1)

f(tij)是時(shí)間變量tij的任意函數(shù),uj(tij)是與tij有關(guān)的水平2隨機(jī)效應(yīng),e0ij是水平1的模型殘差。

因BI值最大為100(即天花板效應(yīng)),故f(tij)固定效應(yīng)預(yù)測(cè)值不應(yīng)大于100。對(duì)于具有天花板效應(yīng)的數(shù)據(jù),這是衡量最優(yōu)模型時(shí)最關(guān)注的限制條件之一。又根據(jù)圖1顯示的拋物線狀變化規(guī)律,故考慮擬合以下4種具有類似圖像的函數(shù)。

(2)

f(tij)=β0+β1tij+β2/(tij+1)

(3)

f(tij)=β0+β1tij+β2log(tij+1)

(4)

f(tij)=β0+β1log(tij+1)

(5)

表1 傷員康復(fù)趨勢(shì)的三角矩陣

水平2隨機(jī)效應(yīng)包含上述模型中截距估計(jì)值β0的隨機(jī)效應(yīng)u0j和斜率估計(jì)值β1的隨機(jī)效應(yīng)u1j。假設(shè)它們是服從均數(shù)為0,方差為Ω的聯(lián)合正態(tài)分布的隨機(jī)變量。

個(gè)體BI均數(shù)由β0+u0j估計(jì),個(gè)體BI隨時(shí)間變化的斜率由β1+u1j估計(jì)。

2.BI變化的平臺(tái)期確定

從式(1)～(5)中找出最優(yōu)模型,以估計(jì)BI預(yù)測(cè)值均數(shù)變化曲線,再?gòu)脑撉€上找出變化速度開始減緩的時(shí)間點(diǎn)tc,定義為平臺(tái)期。由于到達(dá)該點(diǎn)后多數(shù)傷員BI增速變緩,指示康復(fù)干預(yù)效果增速變緩,故將該點(diǎn)作為康復(fù)干預(yù)的“轉(zhuǎn)折事件”點(diǎn)。假設(shè)確定平臺(tái)期BI值為90,傷員A入院后15天BI值達(dá)到90,則該傷員有觀察到的干預(yù)轉(zhuǎn)折事件和到達(dá)時(shí)間t。

3.治療干預(yù)情況到達(dá)平臺(tái)期時(shí)間和影響因素分析

若BI首次測(cè)量值大于平臺(tái)期BI,則傷員到達(dá)平臺(tái)期事件所需時(shí)間t視為缺失值;若末次測(cè)量值小于平臺(tái)期BI,則t為截尾數(shù)據(jù);若某次測(cè)量值等于平臺(tái)期BI,則t為該次測(cè)量時(shí)間(完全數(shù)據(jù));否則采用內(nèi)插法計(jì)算t(完全數(shù)據(jù))。用Cox比例風(fēng)險(xiǎn)回歸模型分析協(xié)變量對(duì)BI到達(dá)平臺(tái)期所需時(shí)間的影響,模型表達(dá)如下:

hi(t)=h0(t)exp(β1xi1+β2xi2+…+βmxim)

(6)

其中,hi(t)為第i例傷員在時(shí)刻t到達(dá)平臺(tái)期的概率函數(shù),協(xié)變量X=(xi1,xi2,…,xim)′為可能影響t的m個(gè)因素,β=(β1,β2,…,βm)′為回歸系數(shù),h0(t)為所有協(xié)變量取值為0時(shí)到達(dá)平臺(tái)期的概率函數(shù)。

數(shù)據(jù)整理用SAS9.2,多水平模型用MLwiN 2.30,Cox模型用SPSS 21軟件。

結(jié) 果

測(cè)量次數(shù)太少的傷員BI隨時(shí)間變化規(guī)律不明顯,故本文只分析測(cè)量次數(shù)大于等于3的155例傷員數(shù)據(jù)。

1.BI時(shí)間變化模型比較

以實(shí)際測(cè)量開始時(shí)間和最長(zhǎng)時(shí)間(0天,464天)為X值估計(jì)得到表2中固定效應(yīng)(即不同時(shí)點(diǎn)的BI均數(shù))預(yù)測(cè)值范圍。此范圍不應(yīng)該大于100,這是當(dāng)受到天花板效應(yīng)的限制、選擇最優(yōu)模型時(shí)最重要的考量之一。可知只有式(5)的BI最大預(yù)測(cè)值在實(shí)際取值范圍內(nèi)。另外,擬合優(yōu)度評(píng)價(jià)的常用指標(biāo)可以作為選擇模型的參考,又由表2可知各個(gè)模型的-2Log Likelihood和AIC差異均不大,故在此不著重考慮。

由圖2可進(jìn)一步看出式(5)預(yù)測(cè)效果與實(shí)際情況最為吻合。又式(5)標(biāo)準(zhǔn)化殘差正態(tài)分?jǐn)?shù)圖顯示水平1和水平2的隨機(jī)效應(yīng)均近似正態(tài)分布,提示擬合式(5)模型合理。故采用式(5)作為擬合BI變化趨勢(shì)的模型。

表2 兩水平隨機(jī)系數(shù)模型擬合優(yōu)度

圖2 模型(2)～(5)估計(jì)的BI均數(shù)隨時(shí)間變化曲線與實(shí)際均數(shù)散點(diǎn)圖對(duì)比

模型(5)各參數(shù)估計(jì)值如下:

BIij～N(XB,Ω)

BIij=β0ijcons+β1jlog(tij+1)

β0ij=26.612(1.478)+u0j+e0ij

β1j=11.764(0.377)+u1j

[e0ij]～N(0,Ωe):Ωe=[147.587(9.256)]

-2loglikelihood(IGLSDeviance)=6761.381(808of808casesinuse)。

此例定義BI平臺(tái)期是隨時(shí)間變化速度開始減緩的時(shí)間點(diǎn),我們用直尺逼近法確定此點(diǎn)。在圖2中模型(5)擬合圖加上網(wǎng)格背景,用直尺逼近曲線變緩且重合較多的線段(圖中線段標(biāo)示),得傷員平臺(tái)期BI約為46(圖中圓點(diǎn)標(biāo)示),進(jìn)而得各傷員到達(dá)平臺(tái)期所需時(shí)間t。剔除t缺失的數(shù)據(jù)后傷員為138例,其中126例(91%)t為完全數(shù)據(jù),12例(9%)為截尾數(shù)據(jù)。

2.影響傷員治療干預(yù)情況到達(dá)平臺(tái)期的因素

Cox模型單因素分析中有統(tǒng)計(jì)學(xué)意義的協(xié)變量有性別、受傷至接受康復(fù)治療時(shí)間、主診斷類別、總住院時(shí)間?？祻?fù)干預(yù)項(xiàng)目數(shù)在單因素分析中無(wú)統(tǒng)計(jì)學(xué)意義,但作為本例主要目的指標(biāo),仍納入多因素分析。向前逐步法結(jié)果見(jiàn)表3。

結(jié)果顯示主診斷類別及住院時(shí)間相同時(shí),男性、接受多項(xiàng)康復(fù)干預(yù)的傷員BI到達(dá)平臺(tái)期的概率比女性和接受較少干預(yù)者大;性別及康復(fù)干預(yù)項(xiàng)目數(shù)相同時(shí),住院時(shí)間越長(zhǎng)和有脊柱骨折的傷員BI到達(dá)平臺(tái)期的概率比參照組小。受傷被困時(shí)間及受傷至接受康復(fù)治療時(shí)間的聯(lián)合檢驗(yàn)顯示時(shí)間越長(zhǎng),BI到達(dá)平臺(tái)期的概率降低。用ROC曲線評(píng)價(jià)模型預(yù)測(cè)效果,以模型預(yù)測(cè)的傷員到達(dá)平臺(tái)期的概率為檢驗(yàn)變量,實(shí)際是否到達(dá)平臺(tái)期為狀態(tài)變量,計(jì)算ROC曲線下面積為0.694(Z=1.9596,P=0.034)。

表3 Cox模型參數(shù)估計(jì)及檢驗(yàn)結(jié)果

*主診斷類別:二分類變量,包括脊柱骨折及其他診斷類別,參考水平是其他診斷類別。

討 論

本文數(shù)據(jù)為接受康復(fù)干預(yù)的蘆山震區(qū)傷員轉(zhuǎn)歸情況的重復(fù)觀察值,BI在傷員水平有聚集性、觀察時(shí)間點(diǎn)數(shù)目不等、部分?jǐn)?shù)據(jù)缺失、基線值及變化速度不等,多水平模型能夠處理具有這種特征的數(shù)據(jù),通過(guò)隨機(jī)效應(yīng)擬合出針對(duì)每一個(gè)個(gè)體的每一個(gè)模型,得到不止一個(gè)模型,在圖像上則體現(xiàn)為具有個(gè)體數(shù)相同數(shù)量的、截距和斜率不同的不平行回歸線[5-8]。多水平模型用隨機(jī)截距反映傷員臨床干預(yù)效果觀察指標(biāo)BI的不同基線情況,用時(shí)間變量t的固定效應(yīng)反映BI隨時(shí)間推移增加的幅度,用t的隨機(jī)效應(yīng)反映不同傷員BI增加幅度的差異,較客觀地反映了圖1描述的BI隨時(shí)間變化的實(shí)際規(guī)律。

本文結(jié)局指標(biāo)BI存在上限100且部分傷員BI堆積在上限處(上刪截?cái)?shù)據(jù)),出現(xiàn)天花板效應(yīng),而一般的模型未考慮此現(xiàn)象導(dǎo)致擬合結(jié)果產(chǎn)生誤差且預(yù)測(cè)值可能超出限制。本文嘗試在多個(gè)模型間做選擇以獲得預(yù)測(cè)值在觀察期內(nèi)不超出實(shí)測(cè)最大值100分的模型,從而可進(jìn)行每位傷員的平臺(tái)期估計(jì)、到達(dá)平臺(tái)期時(shí)間和影響因素分析。目前處理天花板效應(yīng)的方法有Tobin提出的Tobit模型[9,11-12],可通過(guò)Stata軟件的GLLAMM命令實(shí)現(xiàn)。下一步研究需探討能處理隨機(jī)效應(yīng)的Tobit模型,并與本文中簡(jiǎn)單的2水平隨機(jī)效應(yīng)模型比較擬合效果。另外,也可考慮非線性的多水平模型如2水平logistic曲線模型、logistic全局化回歸模型[13]、或光滑樣條非參數(shù)回歸模型是否可提高擬合優(yōu)度[14]。

本文對(duì)平臺(tái)期的定義是BI隨時(shí)間變化速度開始減緩的時(shí)間點(diǎn),而臨床尚沒(méi)有明確定義該時(shí)間點(diǎn),故需根據(jù)模型確定。由于平臺(tái)期并非BI到達(dá)峰值的時(shí)間點(diǎn)或BI曲線的凹凸分界點(diǎn),故不能采用BI趨勢(shì)模型導(dǎo)數(shù)為零或不存在的點(diǎn)即極值點(diǎn),也不能采用BI趨勢(shì)模型二階導(dǎo)數(shù)為零或不存在的點(diǎn)即拐點(diǎn)。考慮是否可采用BI曲線某點(diǎn)切線傾斜角為45度處的時(shí)間點(diǎn),但查閱文獻(xiàn)未找到理論依據(jù)。故本文通過(guò)在最優(yōu)模型估計(jì)的BI均數(shù)隨時(shí)間變化趨勢(shì)圖上,用直尺逼近曲線變緩且重合較多的線段,尋找BI變化速度開始減緩的時(shí)間點(diǎn)作為平臺(tái)點(diǎn)。該法僅為嘗試,后續(xù)研究可驗(yàn)證討論或?qū)ふ腋‘?dāng)?shù)姆椒ā?/p>

綜上,本文資料為臨床觀察性記錄資料,測(cè)量次數(shù)少的傷員比例較多、測(cè)量間隔不統(tǒng)一、觀察數(shù)據(jù)存在缺失值,導(dǎo)致樣本量減少、分析方法復(fù)雜從而研究結(jié)論可靠性降低。這提示我們?cè)谝院笱芯恐袘?yīng)注意研究設(shè)計(jì)和嚴(yán)格的質(zhì)量控制以獲得完整有代表性的數(shù)據(jù)。本文對(duì)解決數(shù)據(jù)及統(tǒng)計(jì)學(xué)問(wèn)題的方法進(jìn)行了嘗試和思考,供其他研究者參考。

[1] 潘培育.綜合性護(hù)理干預(yù)對(duì)慢性阻塞性肺疾病患者生活質(zhì)量的影響.中國(guó)社區(qū)醫(yī)師(醫(yī)學(xué)專業(yè)),2012,14(17):333-335.

[2] 王玲,徐霞,李娟,等.霧化吸入痰熱清治療慢性阻塞性肺疾病急性加重期的護(hù)理觀察.護(hù)士進(jìn)修雜志,2013,28(5):472-473.

[3] Wunderlich MT,Hanhoff T,Goertler M,et al.Release of brain-type and heart-type fatty acid-binding proteins in serum after acute ischaemic stroke.J Neurol,2005,252(6):718-724.

[4] 董佩龍,唐曉波,王健,等.股骨近端防旋髓內(nèi)釘與骨水泥柄人工股骨頭置換治療老年股骨粗隆間骨折的療效分析.中華關(guān)節(jié)外科雜志(電子版),2016,10(1):37-42.

[5] 楊珉,李曉松．醫(yī)學(xué)和公共衛(wèi)生研究常用多水平統(tǒng)計(jì)模型．北京:北京大學(xué)醫(yī)學(xué)出版社,2007:1-61.

[6] 黃坤,倪宗瓚,程薇波.混合線性模型在臨床試驗(yàn)中重復(fù)測(cè)量資料的應(yīng)用.現(xiàn)代預(yù)防醫(yī)學(xué),2005,32(11):1584-1584.

[7] 黃高明,周穎川,梁秋萍.臨床研究中重復(fù)測(cè)量資料的統(tǒng)計(jì)分析方法.廣西醫(yī)科大學(xué)學(xué)報(bào),2000,17(2):264-266.

[8] 李新,包紅.多水平模型在交叉設(shè)計(jì)臨床試驗(yàn)數(shù)據(jù)分析中的應(yīng)用.數(shù)理醫(yī)藥學(xué)雜志,2013,26( 1):100-103

[9] Twisk J,Rijmen F.Longitudinal tobit regression:a new approach to analyze outcome variables with floor or ceiling effects.J Clin Epidemiol,2009,62(9):953-958.

[10] 南登昆主編.康復(fù)醫(yī)學(xué).北京:人民衛(wèi)生出版社,2004:202.

[11] Tobin J.Estimationof Relationships for Limited Dependent Variables.Econometrica,1958,26(1):24-36.

[12] 王煜.中國(guó)居民健康相關(guān)生命質(zhì)量及其對(duì)衛(wèi)生服務(wù)利用影響的研究.北京協(xié)和醫(yī)學(xué)院,2010.

[13] 李雪迎,王琳,朱賽楠,等.全局化模型在描述臨床治療過(guò)程中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2007,24(4):341-344.

[14] 陳生長(zhǎng),徐勇勇,夏結(jié)來(lái).光滑樣條非參數(shù)回歸方法及醫(yī)學(xué)應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),1999,16(6):23-26.

(責(zé)任編輯：張 悅)