縱向監(jiān)測(cè)連續(xù)非隨機(jī)缺失數(shù)據(jù)變系數(shù)模型及其應(yīng)用*

山西醫(yī)科大學(xué)公共衛(wèi)生學(xué)院衛(wèi)生統(tǒng)計(jì)教研室(030001) 徐麗紅 劉志永 劉桂芬 羅天娥

縱向研究是指同一觀察對(duì)象在不同時(shí)間點(diǎn)連續(xù)監(jiān)測(cè)的一種方法。通常情況下，觀測(cè)時(shí)間多是相對(duì)固定的(如隔月檢查、每周測(cè)量并記錄等);但在醫(yī)學(xué)監(jiān)測(cè)過(guò)程中，常會(huì)由于休假、病人有事等實(shí)際問(wèn)題，在原計(jì)劃規(guī)定的時(shí)間點(diǎn)未能獲取到監(jiān)測(cè)值。若忽視時(shí)間點(diǎn)不統(tǒng)一，或收集資料的時(shí)間是在一段時(shí)間內(nèi)獲取的具體問(wèn)題，采用常規(guī)的重復(fù)測(cè)量資料的方差分析，刪除有缺失數(shù)據(jù)的觀察單位，有可能引致估計(jì)偏差加大。變系數(shù)模型就是將其均值參數(shù)和方差成分看作是隨觀測(cè)時(shí)間變化的函數(shù)，利用模式混合模型的基本原理，通過(guò)貝葉斯懲罰樣條構(gòu)建變系數(shù)模型，從而獲得更精確的參數(shù)估計(jì)。該方法不但允許均值參數(shù)隨時(shí)間變化，而且允許模型中方差成分隨缺失時(shí)間變化，實(shí)際應(yīng)用中具有較大的靈活性。

變系數(shù)模型的基本原理

假定監(jiān)測(cè)資料{Yi(t)}和{Xi(t)}分別表示第i個(gè)體在t時(shí)刻的反應(yīng)變量與協(xié)變量值，且Y為連續(xù)型隨機(jī)變量，則有:

{Yi(t)|Xi(t)}∈N(μi(t)，Vi(s，t))，

式中 μi(t)=Xi(t)β，Vi(s，t)=cov(Yi(s)，Yi(t)|Xi(s)，Xi(t))，(s≤t)。令 θ為模型參數(shù)，則{Yi(t)|Xi(t)}∈F(θ;Xi(t))，設(shè)Ui為研究對(duì)象規(guī)定的退出時(shí)間，此時(shí)可將反應(yīng)變量的聯(lián)合密度函數(shù)分解為:

f(Yi(t)，Xi(t)，Ui)=f(Yi(t)|Xi(t)，Ui)f(Ui|Xi(t))(1)模型中{Yi(t)|Xi(t)，Ui}∈F{θ(U);Xi(t)}，若反應(yīng)變量為連續(xù)型隨機(jī)變量，則分布函數(shù)記作:

{Yi|Xi，Zi，bi，Ui}∈N［Xiβ(u)+Zibi，Ri{φ(u)}］(2)

{bi|Ui，δi}∈N(0，Gi{φ(u)})，式中Xi:ni×p維固定效應(yīng)矩陣，Zi:ni×q維隨機(jī)效應(yīng)矩陣?？梢?jiàn)變系數(shù)模型中不僅位置參數(shù)隨時(shí)間而變化，且離散參數(shù)也可以是一個(gè)時(shí)間依賴參數(shù)，即βi和φi隨u的變化而變化。

1．參數(shù)估計(jì)

令隨訪時(shí)間分布為f(u|x，π)，Θ為參數(shù)集合，則第i個(gè)體的似然函數(shù)記作

2．貝葉斯懲罰樣條

懲罰樣條函數(shù)是非參數(shù)法的一種，它在參數(shù)設(shè)定方面具有很大靈活性。貝葉斯法可利用樣本信息結(jié)合先驗(yàn)分布得到后驗(yàn)，利用后驗(yàn)分布進(jìn)行推斷。貝葉斯懲罰樣條結(jié)合了以上兩方面優(yōu)點(diǎn)。θ(·)的三次樣條函數(shù)記作:

式中 φ =(η0，η1，γ1，γ2，…，γk)T，γk為平滑參數(shù)，而 vk為節(jié)點(diǎn)，它是U的第k/(k+1)百分位數(shù)。為保證曲線的靈活性應(yīng)盡量選擇較多節(jié)點(diǎn)，但過(guò)多的節(jié)點(diǎn)數(shù)也可能會(huì)造成過(guò)度擬合，Ruppert認(rèn)為采用5-20個(gè)節(jié)點(diǎn)效果較好〔1〕。

3．敏感性分析

假定包含缺失信息的觀測(cè)值與未包含缺失信息的觀測(cè)值，建模時(shí)都具有相同的參數(shù)估計(jì)結(jié)果〔2〕。為了驗(yàn)證在適當(dāng)延長(zhǎng)觀測(cè)時(shí)間后是否仍能得到相似的結(jié)論，實(shí)際問(wèn)題分析中有必要進(jìn)行敏感性分析，因?yàn)樵摷俣o(wú)法采用已觀測(cè)到的數(shù)據(jù)進(jìn)行驗(yàn)證。若在模型后增加敏感項(xiàng)ω(u)(tij－u)+，令模型中無(wú)法證明的部分為 ω(u)，ω(u)=a{max(U) －u}/range(U)，則稱 a為敏感性參數(shù)。

高血壓社區(qū)規(guī)范化管理的變系數(shù)模型實(shí)現(xiàn)

1．隨機(jī)抽取全國(guó)高血壓社區(qū)規(guī)范化管理項(xiàng)目山西分中心負(fù)責(zé)的太原市四個(gè)社區(qū)衛(wèi)生服務(wù)中心1 275名高血壓患者的監(jiān)測(cè)資料。根據(jù)危險(xiǎn)因素量化分層標(biāo)準(zhǔn)，將社區(qū)規(guī)范化管理的高血壓患者分為1級(jí)、2級(jí)和3級(jí)，分別實(shí)施不同的管理程序。1級(jí)和2級(jí)高血壓患者分別在基線調(diào)查后每隔三個(gè)月隨訪一次，一年內(nèi)隨訪四次;3級(jí)高血壓患者隔月隨訪一次，一年內(nèi)隨訪六次。本文選取3級(jí)高血壓患者301名作為研究對(duì)象，共有264名3級(jí)高血壓患者完成了六次隨訪，失訪率為12.3%。由于本次調(diào)查采取患者主動(dòng)前來(lái)接受調(diào)查的方法，許多研究對(duì)象未能按照研究方案規(guī)定的兩個(gè)月一次接受調(diào)查，表現(xiàn)為不規(guī)則的觀測(cè)時(shí)間和連續(xù)的缺失時(shí)間。分別以收縮壓或舒張壓(mmHg)作為反應(yīng)變量，年齡(歲)、性別、高血壓病程(年)和社區(qū)中心作為自變量，構(gòu)建不同參數(shù)隨觀測(cè)時(shí)間變化的變系數(shù)模型。

采用 SAS9.2整理分析數(shù)據(jù)集，通過(guò) Win-BUGS14.0建立變系數(shù)模型，并通過(guò) R軟件中的R2WINBUGS包來(lái)調(diào)用WinBUGS程序，通過(guò) Gibbs抽樣進(jìn)行MCMC迭代，迭代進(jìn)行15000次，退火5000次。根據(jù)迭代序列圖可見(jiàn)所有參數(shù)的迭代歷史軌跡呈現(xiàn)“毛毛蟲(chóng)”狀，概率密度均呈現(xiàn)中間高，兩邊低的倒鐘形分布，且自相關(guān)圖逐漸趨近于0，說(shuō)明參數(shù)值達(dá)到收斂。最終分析結(jié)果輸出到R軟件中。

圖1(1)～(7)分別給出了截距、年齡、性別、高血壓病程、社區(qū)二、社區(qū)三和社區(qū)四的回歸系數(shù)隨時(shí)間的變化情況。由圖可見(jiàn)，無(wú)論收縮壓還是舒張壓，年齡參數(shù)隨觀測(cè)時(shí)間延長(zhǎng)呈下降趨勢(shì)，表明年齡較大的高血壓患者血壓控制更理想，且舒張壓值下降的幅度比收縮壓值大;收縮壓的性別參數(shù)隨時(shí)間呈上升趨勢(shì)，而舒張壓的性別參數(shù)隨訪期內(nèi)基本不變;收縮壓的高血壓病程參數(shù)基本保持為0，而舒張壓則上升且為正值，表明高血壓病程長(zhǎng)的患者舒張壓控制不理想，而對(duì)收縮壓的影響則較小;中心2和中心3不論對(duì)收縮壓還是舒張壓，在整個(gè)觀測(cè)時(shí)間內(nèi)基本不變，而中心3則呈現(xiàn)下降趨勢(shì)，說(shuō)明中心3的血壓控制較好。圖1(8)～(10)分別給出隨機(jī)系數(shù)截距、隨機(jī)系數(shù)斜率以及個(gè)體間誤差隨觀測(cè)時(shí)間的變化情況。隨機(jī)系數(shù)截距和斜率隨時(shí)間延長(zhǎng)基本為一個(gè)常數(shù)，收縮壓的個(gè)體間誤差呈現(xiàn)下降趨勢(shì)，而舒張壓則略呈上升趨勢(shì)。

圖1 固定效應(yīng)參數(shù)與隨機(jī)效應(yīng)參數(shù)隨觀測(cè)時(shí)間(天)變化圖

表1 分別以收縮壓和舒張壓為反應(yīng)變量的參數(shù)估計(jì)結(jié)果

表1為由變系數(shù)模型得到的收縮壓和舒張壓的各個(gè)參數(shù)估計(jì)值及其可信區(qū)間。其中收縮壓的參數(shù)，性別的可信區(qū)間未包括0，表明女性的收縮壓低于男性;舒張壓的參數(shù)，性別和中心3的可信區(qū)間未包括0，表明女性舒張壓低于男性，社區(qū)3的舒張壓均值低于社區(qū)1。

分別設(shè)定敏感性參數(shù)a為5和10，比較各個(gè)社區(qū)分性別的血壓估計(jì)值變化情況。結(jié)果見(jiàn)表2?？梢钥闯鋈N條件下的估計(jì)結(jié)果均較為接近。

表2 不同敏感性參數(shù)的分析結(jié)果

討 論

隨訪研究中經(jīng)常會(huì)遇到監(jiān)測(cè)時(shí)間不規(guī)則問(wèn)題。完全數(shù)據(jù)分析中，個(gè)體內(nèi)相關(guān)系數(shù)的變化一般僅影響方差成分，而對(duì)回歸參數(shù)的影響較小。當(dāng)數(shù)據(jù)中包含信息缺失時(shí)，個(gè)體內(nèi)相關(guān)系數(shù)的變化也可能會(huì)引致回歸系數(shù)的變化〔3〕。因此本文建模時(shí)也同時(shí)建立了方差成分隨時(shí)間變化的平滑函數(shù)，以得到更為精確的分析結(jié)果。

由于貝葉斯分析結(jié)果對(duì)先驗(yàn)選取非常敏感，若錯(cuò)誤地選擇先驗(yàn)，可能引致錯(cuò)誤的參數(shù)估計(jì)結(jié)果。因此，在參數(shù)分布未知的情況下，應(yīng)選擇較為保守的無(wú)信息先驗(yàn)。Crainiceanu在平滑參數(shù)估計(jì)中將方差成分的先驗(yàn)分布設(shè)定為Gamma(A，B)，并指出當(dāng)A≤0.01且B≤‖b‖2/2時(shí)，先驗(yàn)分布與后驗(yàn)分布之間相互獨(dú)立〔5〕。當(dāng)不滿足上述條件時(shí)，先驗(yàn)分布就會(huì)對(duì)后驗(yàn)分布有影響。因此本文在先驗(yàn)選取上對(duì)位置參數(shù)設(shè)定了無(wú)信息先驗(yàn)，而將方差成分的先驗(yàn)分布設(shè)定為 Gamma(10－6，10－6)。

WinBUGS軟件是進(jìn)行貝葉斯分析最常用的工具。該軟件通過(guò)Gibbs抽樣和MCMC方法來(lái)進(jìn)行復(fù)雜統(tǒng)計(jì)模型推斷，其優(yōu)點(diǎn)是運(yùn)算速度快，結(jié)果準(zhǔn)確。但它對(duì)數(shù)據(jù)輸入格式要求嚴(yán)格，在數(shù)據(jù)庫(kù)操作方面還無(wú)法與目前常用的SAS和R軟件相提并論，且由于其儲(chǔ)存空間有限，當(dāng)數(shù)據(jù)量較大時(shí)常常會(huì)出現(xiàn)超過(guò)記憶容量的情況，實(shí)際應(yīng)用中受到較大限制。Sibylle〔5〕認(rèn)為R軟件推出的R2WinBUGS包允許在R環(huán)境下運(yùn)行WinBUGS程序，可解決此問(wèn)題。本文通過(guò)四個(gè)社區(qū)高血壓規(guī)范化管理實(shí)例，進(jìn)一步證實(shí)并很好地解決了這一問(wèn)題。

1．Diggle P，Heagerty P，Liang KY，et al．Analysis of Longitudinal Data．New York:Oxford University Press，2002．

2．Daniels M，Hogan，J．Missing data in longitudinal studies:strategies for bayesian modeling and sensitivity analysis．Monographs on Statistics and Applied Probability，2008，101．New York:Chapman ＆ Hall．

3．Li S，Joseph WH．Varying-coefficient models for longitudinal processes with continuous-time informative dropout．Biostatistics，2010，11(1):93-110．

4．Ciprian MC，David R，Wand MP．Bayesian analysis for penalized spline regression using WinBUGS．Johns Hopkins University，Dept．of Biostatistics Working Papers，2007:40-74．

5．Sibylle S，Uwe L，Andrew G．R2WinBUGS:A package for running Win-BUGS from R．Journal of Statistical Software，2005，12(3):1-16．