多水平零膨脹計數模型在骨關節疼痛影響因素中的應用研究*

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

層次或組群結構(hierarchical or clustered structure)數據是多中心臨床試驗、社區規范化管理、流行病學復雜抽樣調查、家系資料縱向研究及重復測量數據分析中經常遇到的問題。不同地區(中心或社區)觀察個體的反應或研究時間結局不僅受個體特征的影響，有些還不可避免地受到其周圍環境的制約，形成組群效應。若采用傳統的分析方法常會低估效應的估計誤差，對影響因素效應會作出不合理或錯誤的推斷。本文針對層次或組群結構計數資料，闡明多水平計數模型原理與方法，以山西省農村部分社區居民骨關節疼痛部位數的影響因素為例，完成SAS軟件實現，為醫藥衛生監測中組群結構數據模型分析提供理論依據。

模型原理與方法

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

1.多水平Poisson回歸模型

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

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

水平1模型記作

其兩水平模型表示為:

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

水平1隨機回歸系數β0j和β1j相對應的是兩個水平2方程見式(3)。

式(4)中，括號內為復合殘差結構即模型的隨機成分:u0j和u1j為水平2殘差，其服從二元正態分布;eij為水平1殘差，其服從Poisson分布;u1jzij為組群與水平1變量zij之間的交互作用。

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

針對實際調查計數數據中出現的零過多現象，Lambert(1992)〔3〕首次提出了零膨脹模型(zero-inflated model)。它是把事件數的發生看成兩種可能的過程，第一種是假定事件發生取值只能為0;第二種是假定對應事件發生數取值為0或正的事件數，且服從計數發生過程。

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

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

式中，wij'和xij'為模型logistic和Poisson部分的協變量;γ'和 β'為協變量 wij和 xij的參數向量;φj和 υj為模型隨機效應，分別服從

3.多水平計數回歸模型估計方法及其軟件應用

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

農村地區居民骨關節疼痛部位數多水平零膨脹計數模型分析

課題小組由統一培訓的醫療衛生技術人員組成，于2010年1月至3月抽取山西省陽城縣9個村，偏關縣14個村居民進行風濕病入戶調查。最終納入滿足入選標準要求的分析對象6737人，包括男性3412人，女性3325人。基于該調查資料，以調查點(村莊)為水平2單位，個體為水平1單位，擬合兩水平計數回歸模型。模型結局變量取每個個體檢查的全身(包括頭、頸、肩等20個部位)骨關節疼痛部位數。解釋變量包括地區、性別、年齡、婚姻狀況、高血壓與否。地區(1=陽城，2=偏關)作為水平2變量，其他4個為水平1變量，并對年齡進行均值中心化處理。6737名居民中5171人(76.8%)未曾有過骨關節疼痛，525人(7.8%)疼痛部位數為1。應變量取值最小為0，最大為14，人均0.55個，結果見表1。

表1 骨關節疼痛部位數頻數分布

據分層模型分析理論，逐個引入可能作為隨機斜率的變量，擬合不同地區影響骨關節疼痛部位數因素的隨機截距和隨機斜率計數模型，建模過程表明性別變量隨機斜率有統計學意義。由資料擬合的Poisson、多水平Poisson和多水平ZIP模型參數估計值及分析結果見表2。

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

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

表2 Poisson、多水平Poisson和多水平ZIP模型估計系數與標準誤

結 論

社區抽樣層次結構數據研究中，較低層次單位常嵌套于較高層次單位中。針對該類數據的非獨立性，可采用非線性多水平模型來擬合結局變量為計數數據的資料。文中介紹的多水平Poisson和多水平ZIP模型是解決組群結構效應計數數據模型分析的方法。

實例可見兩地區村民骨關節疼痛部位數為0者占76.8%，零頻數較多，且資料源于10個村莊，故考慮建立多水平零膨脹模型。Vuong檢驗進一步表明多水平ZIP模型比多水平Poisson回歸模型擬合效果更好，提示集群效應計數資料存在額外零時采用多水平ZIP模型分析是更恰當的選擇。

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

實際問題分析過程中，許多人將計數性質的結局變量轉換為二分類變量進行簡單的分析，將導致丟失計數數據蘊含的部分信息，分析結果不可靠。計數回歸模型應變量要求服從二項、Poisson、負二項等分布，但其前提條件是每次事件發生具有相同權重，而實際問題中這一假設條件很難滿足。例如實例中骨關節疼痛共包括20個部位，每個部位疼痛的輕重程度可能是不同的，即每次事件發生的權重是不同的。有關這一問題的解決方法詳見另文。

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

2．方積乾，陸盈．現代醫學統計學．北京:人民衛生出版社，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．王濟川，謝海義，姜寶法．多層統計分析模型:方法與應用．北京:高等教育出版社，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．