復雜抽樣數據多水平模型分析方法及其應用

于石成廖加強于 妺郭 瑩肖革新△金承剛馮國雙胡躍華馬林茂

復雜抽樣數據多水平模型分析方法及其應用

于石成1廖加強2于 妺1郭 瑩1肖革新1△金承剛3馮國雙1胡躍華1馬林茂1

目的本文通過抽樣調查實例，闡述多階段抽樣、不等抽樣概率和事后分層特性不同產生的復雜抽樣數據，其應用多水平模型分析的原理和方法。方法對我國某省行為危險因素抽樣調查的數據，應用未加權和加權的隨機截距logistic回歸模型分析了某些因素與跌倒性傷害的關系。結果實際分析包括50個區縣（PSU），250個鄉鎮街道（2水平），12086個體（1水平）。未加權估計結果顯示：對跌倒性傷害有統計學影響的變量是健康狀況中等和差、未被雇傭和未婚，年齡為負相關，即年齡越大，發生跌倒性傷害的危險性越小；復雜抽樣2水平logistic回歸分析顯示：對跌倒性傷害有統計學影響的變量與未加權的結果基本一致，但未婚失去了統計學意義。體重指數、性別和受教育程度與跌倒性傷害的發生沒有統計學聯系。結論與未加權的結果比，加權分析對跌倒性傷害有統計學影響的變量基本一致，但加權復雜抽樣PMLE估計的標準誤偏大，結果更保守；對性別的分析發現，加權后的結果符合目前對跌倒性傷害發生機制的認識，因此納入權重的多水平分析方法對該資料可能更合理。

復雜抽樣 多水平模型 多階段抽樣 隨機效應logistic回歸

目前主流統計分析軟件，如SAS，SPSS，Stata，MPlus和SUDAN，都將抽樣權重納入統計分析過程，除可進行復雜抽樣數據的描述性統計分析外，還可進行復雜抽樣數據的多元線性回歸、logistic回歸、Poisson回歸和Cox回歸等，使得復雜抽樣數據的統計推斷方法越來越多地在數據分析中得到應用［1-4］。復雜抽樣數據大多具有層次結構即多水平，其特點是反應變量的分布在個體間不具獨立性，存在地理距離內、行政區劃內或特定空間范圍內的聚集性［5］。多水平模型在醫學領域已有多年的應用，在處理層次結構數據上已發揮了重要的作用［6-8］；復雜抽樣數據的分析既要考慮抽樣權重，又要兼顧數據的層次結構，一般是將權重納入廣義線性混合模型（generalized linearm ixed models，GLMMs）來處理這類數據［9-11］。當GLMMs將抽樣權重納入模型后，使GLMMs能處理復雜抽樣數據，解決了復雜抽樣數據多水平模型統計分析和計算問題。鑒于目前我國全國性的流行病學抽樣調查多采用多階段抽樣設計，并且數據具有層次結構，因此復雜抽樣數據的多水平模型分析已有明顯的應用價值。

多水平復雜抽樣數據的線性模型，Pfeffermann應用偽最大似然估計（pseudo-maximum-likelihood estimation）算法來估計模型的參數，效果很好。但對廣義線性混合模型（generalized linearm ixed models），認為較好的模型估計方法是全偽最大似然估計（full pseudo-maximum-likelihood estimation），它應用自適應積分法（adaptive quadrature）估計模型參數，標準誤的估計采用泰勒線性化三明治估計量［10］。目前國外復雜抽樣數據多水平模型理論和應用研究已有二十幾年的歷史，我國在這方面的理論研究鮮有報告，可能由于數據本身和應用的限制，在醫學研究領域國內還沒有應用復雜抽樣多水平模型來處理流行病學抽樣調查數據。本文介紹了復雜抽樣數據多水平模型分析的原理和方法，并用STATA軟件對我國某省行為危險因素調查的數據擬合復雜抽樣數據多水平模型，并與未考慮復雜抽樣的多水平模型結果進行了比較和解釋，以闡述具有層次結構的復雜抽樣數據應用復雜抽樣數據多水平模型的合理性。

模型基本原理

廣義線性混合模型參數估計構造一個常邊際對數似然函數（usualmarginal log likelihood function），其公式如下：

復雜抽樣數據的廣義線性混合模型參數估計構造一個加權對數偽似然函數（log pseudo-likelihood function w ith weights），其公式如下：

這里，wj＝1/πj，wi｜j＝1/πi｜j，i和j分別代表1水平個體和2水平的層。一般情況下，復雜抽樣數據多水平模型采用牛頓-拉夫遜最大算法（New ton-Raphson maximum algorithm）使加權對數偽似然函數達到最大值，即PMLE（pseudo-maximum-likelihood estimation）參數估計值；其標準誤估計采用泰勒線性化三明治估計量（sandw ich estimator）。

復雜抽樣數據廣義線性混合模型用上式（2）做PMLE估計時，須考慮2水平權重；且不能直接使用1水平的個體權重，這樣可能影響參數估計的準確性［11］，解決方法是對1水平個體權重進行權重的縮放（scaling of weights）。Longford等［12］1995年提出的權重縮放方法1如下：

1998年Pfeffermann等［13］提出了權重縮放方法2如下：

應用實例及結果解釋

我國某省在2010年進行了行為危險因素的調查，該省共有150個區縣，1900個鄉鎮，人口7000多萬。根據計算本調查所需樣本量為15000人，使用了三階段整群抽樣設計，第一階段隨機抽取區縣，第二階段抽取鄉、鎮和街道，第三階段抽取個體。第一階段隨機抽取產生了50個區縣（primary sampling units，PSU）；第二階段，在每個區縣內按鄉、鎮、街道分層，在每層內按概率比例規模抽樣（PPS）方法，整群抽取2個鄉、1個鎮和2個街道，共產生了250個鄉鎮街道；在抽取的每個鄉、鎮和街道中，隨機抽取60名年齡45歲及以上的成年人，共15000個體。該調查將在過去的3個月內發生1次或2次及以上跌倒性傷害為反應變量，記為1，未發生任何跌倒性傷害為0，反應變量為二分類變量，可用復雜抽樣logistic回歸分析該數據；但該數據顯然具有地區和個體的層次結構，因此，合適的統計分析方法應是考慮復雜抽樣2水平logistic回歸模型。這里個體為1水平，250個鄉、鎮、街道為2水平，最高水平區縣（PSU）在這里不作為水平考慮，但在分析中可作為層變量進行調整處理。

本文分析了跌倒后造成的傷害與體重指數（BM I）及一些感興趣變量的關系，因變量為在過去三個月內是否發生過跌倒性傷害，自變量有體重指數、性別、年齡、健康狀況、受教育程度、雇傭狀況和婚姻狀況等7個變量。該數據擬合復雜抽樣2水平隨機截距logistic回歸模型，分析了250個鄉、鎮和街道，由于有缺失數據，實際分析的樣本量為12086例。

調查采用了不等概率抽樣，需要根據抽樣設計對樣本進行抽樣加權，加權過程涉及到個體抽樣權重，未包括不應答權重和事后分層加權。三階段抽樣，用w代表權重，腳注1、2和3分別代表縣、鄉鎮街道和個體各階段的抽樣權重，用i表示某一樣本個體，s代表某一樣本個體所在的層。樣本區縣的抽樣權重（wsi1）值為分層簡單隨機抽樣下樣本區縣抽樣概率的倒數，其計算公式如下：

樣本鄉鎮街道的抽樣權重（wsi2）值為與人口數成比例的PPS抽樣下樣本鄉鎮街道抽樣概率的倒數，用下式計算：

樣本個體的抽樣權重（wsi3）值為樣本個體抽樣概率的倒數。可計算如下：

依據上述各階段抽樣權重，最終樣本個體的抽樣權重為下式（5）：

wsi＝wsi1×wsi2×wsi3＝樣本個體i所在區縣分層抽樣比倒數×

樣本個體i所在鄉鎮街道45歲及以上人口百分比構成，如無法得到實際數據，可用該縣或省的構成（Pop45）計算。

計算鄉鎮街道（2水平單位）πj，個體被抽中的概率πij及它們的權重，即概率的倒數，其計算公式為：

ws2＝wsi1×wsi2＝樣本個體i所在區縣分層抽樣比倒數×

Stata提供了分析復雜抽樣2水平隨機截距logistic回歸模型的分析模塊gllamm和說明書，其未加權最大似然估計（MLE），stata分析語句為：gllamm injury bmi age gender health marriage education employed，i（ID_level_2）link（logit）fam ily（binom）nip（12）adapt

穩健標準誤估計（Robust standard errors）：gllamm，robust

其加權偽最大似然估計（PMLE）stata分析語句為：

gllamm injury bm i age gender health marriage education employed，i（ID_level_2）cluster（w t2）link（logit）family（binom）pweight（w t1_sw）nip（12）adapt

調整PSU后的穩健標準誤估計：gllamm，robust cluster（psu）

其中，因變量為二分類injury，自變量有：體重指數BM I（啞變量）、健康狀況health（啞變量）、婚姻狀況marriage（啞變量）、受教育程度education（啞變量）和雇傭狀況employed（啞變量）。i指出2水平變量為ID_level_2；cluster定義2水平變量的權重為w t2；link指出連接函數為logit；family指出擬合二項分布（binom），若擬合Poisson分布，則family（Poisson）；pweight定義1水平個體權重為w t1_sw，其為經縮放后的值。前面提到1水平權重不能直接應用，要進行權重的縮放，下面是應用式（3）和（4）兩種權重縮放方法的stata語句。

方法1：

方法2：

應用2水平隨機截距logistic回歸模型分析了鄉鎮街道（2水平）的個體（1水平）跌倒性傷害與感興趣的因素關系，用區縣變量（PSU）分層調整。分析時沒有將區縣作為一個水平來對待，主要考慮區縣的變異可能不大，但在標準誤估計時，調整了區縣（PSU）的聚集效應。由于缺失數據的存在，最終的分析在50個PSU，250個鄉鎮（2水平）和12086個體（1水平）中進行；最終結果報告了未加權2水平logistic回歸的參數估計-基于模型標準誤估計和穩健標準誤估計，以及應用權重縮放方法的PMLE參數估計，結果見表1。從表1可見，未加權2水平logistic回歸的參數估計，基于模型和穩健標準誤估計結果很接近，對跌倒性傷害有統計學顯著影響的變量：健康狀況中等和差、未被雇傭和未婚，年齡為負相關，即年齡越大，發生跌倒性傷害的危險性越小；復雜抽樣2水平logistic回歸參數PMLE估計，權重縮放方法1和2估計結果非常接近；與未加權的方法比，對跌倒性傷害有統計學顯著影響的變量與未加權的結果基本一致，但未婚失去了統計顯著性，兩種縮放方法估計的P值在0.05～0.10之間。其它分析變量，如體重指數、性別和受教育程度與跌倒性傷害的發生沒有統計學聯系。加權復雜抽樣PMLE估計，大多數參數估計值比未加權估計增加，但其標準誤比未加權估計值增加明顯，導致大部分變量P值增大，估計保守。

表1 2水平隨機截距logistic回歸和復雜抽樣2水平隨機截距logistic回歸擬合

討 論

流行病學抽樣調查一般應用多階段抽樣，第一階段先抽取地區或群（clusters），第二階段抽取群下面的亞群（subclusters）；最后階段抽取基本抽樣單位，如個體。這種抽樣方法產生了多水平數據，基本抽樣單位為1水平或最低水平，最高水平為PSU。對這類資料的分析不但考慮復雜抽樣的權重（包括抽樣權重、不應答權重和事后分層權重），而且考慮數據的多水平結構，即應用復雜抽樣數據多水平模型分析方法。本文應用我國某省行為危險因素調查的數據，應用Rabe-Hesketh描述的考慮權重的全PMLE參數估計方法，通過適應積分法估計任何水平模型的參數和標準誤。Stata軟件gllamm程序分析發現：健康狀況中等和差、未被雇傭和未婚與跌倒性傷害的關系與以前的文獻或研究一致［14］；但一般認為年齡大發生跌倒性傷害的可能性大［15］，但本結果與之相反，這也許與本研究選取的個體年齡大有關。文獻認為男性較女性更易發生跌倒性傷害［16］，本研究未加權的分析結果是男性發生跌倒性傷害的危險性小；加權的2水平隨機截距logistic回歸模型顯示男性與跌倒性傷害是正相關；雖然加權和未加權的結果均沒有統計顯著性，但加權后的結果合理，符合目前對跌倒性傷害發生機制的認識。這也部分說明了復雜抽樣的層次結構數據，分析時納入權重的多水平分析方法的合理性。

復雜抽樣數據的統計分析也有百年的歷史［17］，1990年代以后更有新統計理論和方法發展，復雜抽樣的統計方法已包括：列聯表資料的對數線性模型和相關的方法、廣義線性模型（logistic回歸和Poisson回歸）、生存分析、一般線性混合模型、結構方程模型、隱變量模型，這些方法充分考慮了抽樣權重，對參數標準誤和可信區間的估計更準確。目前SAS、Stata、SPSS和SUDAN統計軟件包都包含了復雜抽樣數據的統計描述、一般線性回歸、logistic回歸（二分類、有序和無序）、Poisson回歸（零膨脹、負二項）、生存分析等，可以滿足大部分復雜抽樣數據的統計分析。

但是目前流行的統計軟件還都沒有包括復雜抽樣多水平模型擬合程序。當抽樣權重納入模型時，一般構造偽似然函數（pseudo-maximum-likelihood）來估計模型參數。在實踐中的一個主要問題是大多數流行病學抽樣調查僅給出了基本抽樣單位或1水平的抽樣權重，而沒有高水平的權重；但在進行復雜抽樣數據多水平模型擬合時，需要高水平的抽樣權重。因此，在今后的流行病學調查設計中和實施時，應注意收集計算各水平權重指標數據，以便在數據分析時做加權處理。

1.Rao JNK.Interplay between sample survey theory and practice；anappraisal.Survey Methodology，2005，31：117-138.

2.呂筠，何平平，李立明.復雜抽樣調查數據實例分析.中華流行病學雜志，2008，29（8）：832.

3.繆凡，童峰.復雜抽樣數據的logistic回歸分析方法及其應用.中國衛生統計，2008，25（6）：577-579.

4.胡躍華，匡翔宇，金承剛，等.復雜抽樣Poisson回歸分析方法及應用.中國衛生統計，2012，29（5）：650-653.

5.楊珉，李曉松主編.醫學和公共衛生研究常用多水平統計模型.北京，北京大學醫學出版社，2007.

6.賈改珍，閆陽，徐天和，等.多水平模型在大學生預防艾滋病健康教育影響因素分析中的應用.中國衛生統計，2013，30（1）：37-39.

7.Gebremariam MK，Andersen LF，Bielland M，et al.Does the school food environment influence the dietary behaviours of Norwegian 11-yearolds.The HEIA study.Scand JPublic Health，2012，40（5）：491-497.

8.Nansel TR，Lipsky LM，Lannotti RJ.Cross-sectional and longitudinal relationships of bodymass index with glycemic control in children andadolescentswith type 1 diabetesmellitus，2013，100（1）：126-132.

9.Steven G.Heeringa，Wagner J，Torres M，et al.Sample designs and samplingmethods for the Collaborative Psychiatric Epidem iology Studies（CPES）.Int.J.Methods Psychiatr.Res.，2004，13（4）：221-240.

10.Sophia Rabe-Hesketh.Multilevelmodeling of complex survey data.J. R.Statist.Soc.，2006，169：805-827.

11.Moshe Feder，Gad Nathan，Danny Pfeffermann.Survey Methodology，2000，26（1）：53-65.

12.Longford NT.Model-basedmethods for analysis of data from 1990 NAEP Trial State Assessment.Research and Development Report NCES 95-696.Washington DC：National Center for Education Statistics.

13.Pfeffermann D，Holmes CJ，Goldstein DJ，et al.Weighting for unequal selection probabilities in multilevelmodels.J.R.Statist.Soc.B，1998，60：23-40.

14.Roe B，Howell F，Riniotis K，et al.Older people and falls：health status，quality of life，lifestyle，care networks，prevention and views on service use follow ing a recent fall.JClin Nurs，2009，18：2261-2272.

15.Hausdorff JM，Rios DA，Edelberg HK.Gait variability and fall risk in community-living older adults：A 1-year prospective study.A rch Phys Med Rehabil，2001，82：1050-1056.

16.Stevens JA，Sogolow ED.Gender differences for non-fatal unintentional fall related injuries among older adults.Injury Prevention，2005，11：115-119.

17.Fisher RA.StatisticalMethods for Research Work.Oliver and Boyd，Edinburgh，1925.

（責任編輯：劉 壯）

App lication of M ultilevel M odeling to Com plex Sam ple Survey Data

Yu Shicheng，Liao Jiaqiang，Yu Mo，et al（Chinese Center for Disease Control and Prevention（102206），Beijing）

ObjectiveTo illustrate the principal and application ofmultilevelmodeling of complex survey data thatwere derived from multistage sampling，unequal sampling probabilities and different features of post-stratification.MethodsWeighted and un-weighted random intercept logistic regressionmodelswere applied to complex survey data of behavioral risk factors in a province to look at the association of fall injuries w ith some factors of interest.ResultsThere were 12086 subjects（level 1）aged 45 years or above nested w ithin 250 villages，towns and sub-districts（level 2）from 50 counties/districts（PSU）.Un-weighted results showed that variables significantly and positively associated w ith the risk of fall injurieswere fair or poor health，unemployed situation，unmarried；age was significantly and negatively associated w ith the risk of fall injuries，or one less likely got injured when getting older.The results from 2-level random intercept logisticmodel demonstrated that the variables associated w ith the risk of fall injuries were sim ilar to those from un-weighted models，but the variable of unmarried m itigated its significance to be insignificant.Body mass index，beingmale，educational levelwere notassociated w ith the risk of fall injuries from the analyses.ConclusionIn contrast to the results from un-weighted methods，statistically significant variables from weightedmethodswere analogous to those from weighted ones；however，estimates using full pseudo-maximum-likelihood estimation（PMLE）weremore conservative as opposed to un-weighted ones.As for gender，weighted result was in consistent w ith the currentunderstanding of themechanism for the developmentof fall injuries，therefore，it soundedmore reasonable to employ multilevel modeling for the complex survey data.

Complex survey data；Multilevelmodel；Multistage sampling；Random intercept logistic regression

*：淮河流域癌癥綜合防治項目（1310800003）

1.中國疾病預防控制中心公共衛生監測與信息服務中心（102206）

2.四川大學華西公共衛生學院衛生統計教研室（610041）

3.北京師范大學社會發展與公共政策學院（100875）

△通信作者：肖革新，E-mail：biocomputer＠126.com