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

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

·應(yīng)用研究·

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

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

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

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

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

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

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

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

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

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

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

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

資 料

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

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

圖1 306例傷員個體的BI隨時間變化趨勢圖

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

統(tǒng)計分析方法

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

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

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

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

(1)

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

因BI值最大為100(即天花板效應(yīng)),故f(tij)固定效應(yīng)預(yù)測值不應(yīng)大于100。對于具有天花板效應(yīng)的數(shù)據(jù),這是衡量最優(yōu)模型時最關(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ù)趨勢的三角矩陣

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

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

2.BI變化的平臺期確定

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

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

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

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

(6)

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

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

結(jié) 果

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

1.BI時間變化模型比較

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

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

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

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

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

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

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

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

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

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

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

討 論

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

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

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

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

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

[2] 王玲,徐霞,李娟,等.霧化吸入痰熱清治療慢性阻塞性肺疾病急性加重期的護理觀察.護士進修雜志,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)計模型．北京:北京大學(xué)醫(yī)學(xué)出版社,2007:1-61.

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

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

[8] 李新,包紅.多水平模型在交叉設(shè)計臨床試驗數(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ān)生命質(zhì)量及其對衛(wèi)生服務(wù)利用影響的研究.北京協(xié)和醫(yī)學(xué)院,2010.

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

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

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