多水平零膨脹計(jì)數(shù)模型在骨關(guān)節(jié)疼痛影響因素中的應(yīng)用研究*

田嬌妮 劉桂芬△ 張軍鋒，2 羅天娥

層次或組群結(jié)構(gòu)(hierarchical or clustered structure)數(shù)據(jù)是多中心臨床試驗(yàn)、社區(qū)規(guī)范化管理、流行病學(xué)復(fù)雜抽樣調(diào)查、家系資料縱向研究及重復(fù)測(cè)量數(shù)據(jù)分析中經(jīng)常遇到的問(wèn)題。不同地區(qū)(中心或社區(qū))觀察個(gè)體的反應(yīng)或研究時(shí)間結(jié)局不僅受個(gè)體特征的影響，有些還不可避免地受到其周?chē)h(huán)境的制約，形成組群效應(yīng)。若采用傳統(tǒng)的分析方法常會(huì)低估效應(yīng)的估計(jì)誤差，對(duì)影響因素效應(yīng)會(huì)作出不合理或錯(cuò)誤的推斷。本文針對(duì)層次或組群結(jié)構(gòu)計(jì)數(shù)資料，闡明多水平計(jì)數(shù)模型原理與方法，以山西省農(nóng)村部分社區(qū)居民骨關(guān)節(jié)疼痛部位數(shù)的影響因素為例，完成SAS軟件實(shí)現(xiàn)，為醫(yī)藥衛(wèi)生監(jiān)測(cè)中組群結(jié)構(gòu)數(shù)據(jù)模型分析提供理論依據(jù)。

模型原理與方法

多水平模型〔1〕處理具有層次結(jié)構(gòu)特征數(shù)據(jù)時(shí)不僅考慮了資料誤差的層次性，且可將傳統(tǒng)模型的誤差隨機(jī)項(xiàng)分解到與數(shù)據(jù)層次結(jié)構(gòu)相應(yīng)的水平上。通過(guò)估計(jì)個(gè)體水平的誤差，考慮解釋變量對(duì)方差的影響，既充分利用了個(gè)體水平內(nèi)的聚集信息、獲取回歸系數(shù)的有效估計(jì)，又可純化個(gè)體的隨機(jī)誤差〔2〕。它不僅能在層次結(jié)構(gòu)中正確處理模型參數(shù)估計(jì)的問(wèn)題，而且還可同時(shí)分析微觀與宏觀變量效應(yīng)及跨層的交互作用。

1.多水平Poisson回歸模型

Poisson回歸是計(jì)數(shù)回歸的基礎(chǔ)模型，屬?gòu)V義線性模型指數(shù)分布族，其模型連接函數(shù)多采用對(duì)數(shù)函數(shù)。Poisson回歸要求事件的發(fā)生獨(dú)立;事件發(fā)生數(shù)的條件均值等于方差，即滿足等離散(equi-dispersion)特征。設(shè)yij為第j個(gè)群第i個(gè)個(gè)體的觀察值，多水平Poisson模型表示為:

帶有兩個(gè)水平1(解釋變量xij、zij)和1個(gè)水平2(解釋變量wij)的兩水平Poisson回歸模型表示為:

水平1模型記作

其兩水平模型表示為:

式(1)中，yij為第j個(gè)水平2單位中第i個(gè)個(gè)體水平1結(jié)局測(cè)量值，i=1，2，…，n(n 是樣本總量)，j=1，2，…，j(j是水平2的觀察單位數(shù))。式(2)中，α1為水平1模型固定斜率;β0j為水平1模型隨機(jī)截距(j表示水平1截距跨水平2單位的變化);β1j為水平1模型隨機(jī)斜率(j表示水平1斜率跨水平2單位的變化)。

水平1隨機(jī)回歸系數(shù)β0j和β1j相對(duì)應(yīng)的是兩個(gè)水平2方程見(jiàn)式(3)。

式(4)中，括號(hào)內(nèi)為復(fù)合殘差結(jié)構(gòu)即模型的隨機(jī)成分:u0j和u1j為水平2殘差，其服從二元正態(tài)分布;eij為水平1殘差，其服從Poisson分布;u1jzij為組群與水平1變量zij之間的交互作用。

2.多水平零膨脹Poisson回歸模型

針對(duì)實(shí)際調(diào)查計(jì)數(shù)數(shù)據(jù)中出現(xiàn)的零過(guò)多現(xiàn)象，Lambert(1992)〔3〕首次提出了零膨脹模型(zero-inflated model)。它是把事件數(shù)的發(fā)生看成兩種可能的過(guò)程，第一種是假定事件發(fā)生取值只能為0;第二種是假定對(duì)應(yīng)事件發(fā)生數(shù)取值為0或正的事件數(shù)，且服從計(jì)數(shù)發(fā)生過(guò)程。

若計(jì)數(shù)模型構(gòu)建中既考慮組群效應(yīng)，又要討論零頻數(shù)過(guò)多問(wèn)題，就需要用多水平(隨機(jī)效應(yīng))零膨脹Poisson回歸模型(multilevel zero-inflated poisson model，ML_ZIP)〔4，5〕來(lái)解釋組群內(nèi)觀察對(duì)象的非獨(dú)立性。假定yij為第j個(gè)群第i個(gè)個(gè)體的觀察值，多水平零膨脹模型可表示為:

式中，f(·)為Poisson分布的概率密度函數(shù);I(yij)為指示變量，當(dāng) yij=0 時(shí)，I(yij)=1，當(dāng) yij＞0 時(shí)，I(yij)=0;πij為第j個(gè)群第i個(gè)個(gè)體事件發(fā)生數(shù)取值為0的概率。當(dāng)πij=0時(shí)，多水平零膨脹模型退化為多水平Poisson回歸。多水平零膨脹計(jì)數(shù)模型有兩部分組成，連接函數(shù)分別記作:

式中，wij'和xij'為模型logistic和Poisson部分的協(xié)變量;γ'和 β'為協(xié)變量 wij和 xij的參數(shù)向量;φj和 υj為模型隨機(jī)效應(yīng)，分別服從

3.多水平計(jì)數(shù)回歸模型估計(jì)方法及其軟件應(yīng)用

計(jì)數(shù)模型屬非線性模型，積分似然函數(shù)多采用近似估計(jì)。常借助SAS、MLWin和stata等軟件進(jìn)行多水平計(jì)數(shù)模型構(gòu)建。MLWin軟件運(yùn)用Taylor展開(kāi)對(duì)似然函數(shù)進(jìn)行線性化處理，利用迭代廣義最小二乘或限制迭代廣義最小二乘法獲得模型的無(wú)偏估計(jì)。SAS軟件LIMIXED過(guò)程雖可完成隨機(jī)效應(yīng)計(jì)數(shù)回歸模型估計(jì)，但它不提供隨機(jī)效應(yīng)的統(tǒng)計(jì)檢驗(yàn)。而NLMIXED過(guò)程可采用數(shù)值積分近似估計(jì)邊際似然函數(shù)〔6〕，Pinheiro和Bates證明適應(yīng)性高斯求積法結(jié)果更可靠〔7〕，所產(chǎn)生的邊際積分似然函數(shù)常用二元準(zhǔn)牛頓算法進(jìn)行極大化。數(shù)值積分近似估計(jì)的優(yōu)點(diǎn)在于，它既可基于原分析數(shù)據(jù)估計(jì)似然函數(shù)，其-2LL也可用于模型的似然比檢驗(yàn);同時(shí)也能提供隨機(jī)效應(yīng)的顯著性檢驗(yàn)。本資料通過(guò)SAS9.2 NLMIXED過(guò)程實(shí)現(xiàn)。

農(nóng)村地區(qū)居民骨關(guān)節(jié)疼痛部位數(shù)多水平零膨脹計(jì)數(shù)模型分析

課題小組由統(tǒng)一培訓(xùn)的醫(yī)療衛(wèi)生技術(shù)人員組成，于2010年1月至3月抽取山西省陽(yáng)城縣9個(gè)村，偏關(guān)縣14個(gè)村居民進(jìn)行風(fēng)濕病入戶調(diào)查。最終納入滿足入選標(biāo)準(zhǔn)要求的分析對(duì)象6737人，包括男性3412人，女性3325人?；谠撜{(diào)查資料，以調(diào)查點(diǎn)(村莊)為水平2單位，個(gè)體為水平1單位，擬合兩水平計(jì)數(shù)回歸模型。模型結(jié)局變量取每個(gè)個(gè)體檢查的全身(包括頭、頸、肩等20個(gè)部位)骨關(guān)節(jié)疼痛部位數(shù)。解釋變量包括地區(qū)、性別、年齡、婚姻狀況、高血壓與否。地區(qū)(1=陽(yáng)城，2=偏關(guān))作為水平2變量，其他4個(gè)為水平1變量，并對(duì)年齡進(jìn)行均值中心化處理。6737名居民中5171人(76.8%)未曾有過(guò)骨關(guān)節(jié)疼痛，525人(7.8%)疼痛部位數(shù)為1。應(yīng)變量取值最小為0，最大為14，人均0.55個(gè)，結(jié)果見(jiàn)表1。

表1 骨關(guān)節(jié)疼痛部位數(shù)頻數(shù)分布

據(jù)分層模型分析理論，逐個(gè)引入可能作為隨機(jī)斜率的變量，擬合不同地區(qū)影響骨關(guān)節(jié)疼痛部位數(shù)因素的隨機(jī)截距和隨機(jī)斜率計(jì)數(shù)模型，建模過(guò)程表明性別變量隨機(jī)斜率有統(tǒng)計(jì)學(xué)意義。由資料擬合的Poisson、多水平Poisson和多水平ZIP模型參數(shù)估計(jì)值及分析結(jié)果見(jiàn)表2。

影響關(guān)節(jié)疼痛部位數(shù)因素的回歸擬合結(jié)果表明，三個(gè)模型中Poisson回歸的-2LL、BIC及AIC值最大，似然比檢驗(yàn)結(jié)果可見(jiàn)多水平Poisson回歸模型優(yōu)于一般Poisson回歸(χ2=856，v=3，P ＜0.001)。由于多水平Poisson和多水平ZIP屬非嵌套模型，故采用Vuong檢驗(yàn)進(jìn)行模型的比較與選擇，結(jié)果表明V＞1.96，可認(rèn)為多水平ZIP擬合關(guān)節(jié)疼痛部位數(shù)影響結(jié)果最優(yōu)。

多水平ZIP回歸部分的隨機(jī)項(xiàng)φj和vj的方差有統(tǒng)計(jì)學(xué)意義，表明不同村莊居民關(guān)節(jié)是否疼痛以及疼痛部位數(shù)差別有統(tǒng)計(jì)學(xué)意義。logistic回歸估計(jì)參數(shù)表明年齡、婚姻狀況和高血壓是影響關(guān)節(jié)疼痛與否的主要因素，提示年齡越大骨關(guān)節(jié)疼痛的風(fēng)險(xiǎn)越大;離異的居民骨關(guān)節(jié)疼痛風(fēng)險(xiǎn)是在婚居民的1.71倍;患高血壓的居民是未患高血壓居民的1.68倍。而居民生活的地區(qū)、年齡、地區(qū)和性別交互作用是影響關(guān)節(jié)疼痛部位數(shù)的主要因素，實(shí)例中陽(yáng)城地區(qū)居民骨關(guān)節(jié)疼痛數(shù)是偏關(guān)地區(qū)的1.65倍，年齡越大關(guān)節(jié)疼痛部位數(shù)平均越多;生活環(huán)境地區(qū)與性別之間的跨層交互作用為正值，表明居民平均關(guān)節(jié)疼痛部位數(shù)之性別差以偏關(guān)地區(qū)較多。例如，偏關(guān)地區(qū)女性居民關(guān)節(jié)疼痛數(shù)的期望值是男性的exp(0.418-0.280)=1.15倍。

表2 Poisson、多水平Poisson和多水平ZIP模型估計(jì)系數(shù)與標(biāo)準(zhǔn)誤

結(jié) 論

社區(qū)抽樣層次結(jié)構(gòu)數(shù)據(jù)研究中，較低層次單位常嵌套于較高層次單位中。針對(duì)該類數(shù)據(jù)的非獨(dú)立性，可采用非線性多水平模型來(lái)擬合結(jié)局變量為計(jì)數(shù)數(shù)據(jù)的資料。文中介紹的多水平Poisson和多水平ZIP模型是解決組群結(jié)構(gòu)效應(yīng)計(jì)數(shù)數(shù)據(jù)模型分析的方法。

實(shí)例可見(jiàn)兩地區(qū)村民骨關(guān)節(jié)疼痛部位數(shù)為0者占76.8%，零頻數(shù)較多，且資料源于10個(gè)村莊，故考慮建立多水平零膨脹模型。Vuong檢驗(yàn)進(jìn)一步表明多水平ZIP模型比多水平Poisson回歸模型擬合效果更好，提示集群效應(yīng)計(jì)數(shù)資料存在額外零時(shí)采用多水平ZIP模型分析是更恰當(dāng)?shù)倪x擇。

計(jì)數(shù)數(shù)據(jù)模型很多，除本文提及的多水平Poisson和多水平ZIP模型外，Greene和Long〔8〕認(rèn)為，多水平零膨脹負(fù)二項(xiàng)回歸模型(ZINB)有時(shí)比ZIP模型更具優(yōu)勢(shì)。若將文中公式(5)中替換為負(fù)二項(xiàng)分布，也可得到多水平ZINB模型，同理probit模型也可代替logistic模型作為模型中的二值概率函數(shù)。進(jìn)一步說(shuō)明實(shí)際問(wèn)題研究可按文中介紹的原理和軟件實(shí)現(xiàn)，完成具體的多水平計(jì)數(shù)模型構(gòu)建。

實(shí)際問(wèn)題分析過(guò)程中，許多人將計(jì)數(shù)性質(zhì)的結(jié)局變量轉(zhuǎn)換為二分類變量進(jìn)行簡(jiǎn)單的分析，將導(dǎo)致丟失計(jì)數(shù)數(shù)據(jù)蘊(yùn)含的部分信息，分析結(jié)果不可靠。計(jì)數(shù)回歸模型應(yīng)變量要求服從二項(xiàng)、Poisson、負(fù)二項(xiàng)等分布，但其前提條件是每次事件發(fā)生具有相同權(quán)重，而實(shí)際問(wèn)題中這一假設(shè)條件很難滿足。例如實(shí)例中骨關(guān)節(jié)疼痛共包括20個(gè)部位，每個(gè)部位疼痛的輕重程度可能是不同的，即每次事件發(fā)生的權(quán)重是不同的。有關(guān)這一問(wèn)題的解決方法詳見(jiàn)另文。

1．Goldstein H．Multilevel models 3nd．London:Amold，2003．

2．方積乾，陸盈．現(xiàn)代醫(yī)學(xué)統(tǒng)計(jì)學(xué)．北京:人民衛(wèi)生出版社，2002．

3．Lambert D．Zero-inflated Poisson regression with an application to defects in manufacturing．Technometrics，1992，1(34):1-14．

4．Siddiqui O．Modeling clustered count and survival data with an application to a school-based smoking prevention study［thesis］．Chicago:Center for health statistics，University of Illinois at Chicago，1996．

5．Wang K，Yau K，Lee A．A zero-inflated Poisson mixed model to analyze diagnosis related groups with majority of same-day hospital stays．Computer Methods and Program in biomedicine，2002，68:195-203．

6．王濟(jì)川，謝海義，姜寶法．多層統(tǒng)計(jì)分析模型:方法與應(yīng)用．北京:高等教育出版社，2008．

7．Pinheiro J，Bates D．Approximations to the Log-likelihood Function in the Nonlinear Mixed-effects Model．Journal of Computational and Graphical Statistics，1995，1(4):12-35．