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

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

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

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

{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為研究對象規(guī)定的退出時間，此時可將反應(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ù)型隨機變量，則分布函數(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維隨機效應(yīng)矩陣?？梢娮兿禂?shù)模型中不僅位置參數(shù)隨時間而變化，且離散參數(shù)也可以是一個時間依賴參數(shù)，即βi和φi隨u的變化而變化。

1．參數(shù)估計

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

2．貝葉斯懲罰樣條

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

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

3．敏感性分析

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

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

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

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

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

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

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

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

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

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

討 論

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

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

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

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．