復雜抽樣數(shù)據(jù)多水平模型分析方法及其應用

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

復雜抽樣數(shù)據(jù)多水平模型分析方法及其應用

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

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

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

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

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

模型基本原理

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

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

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

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

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

應用實例及結果解釋

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

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

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

樣本鄉(xiāng)鎮(zhèn)街道的抽樣權重（wsi2）值為與人口數(shù)成比例的PPS抽樣下樣本鄉(xiāng)鎮(zhèn)街道抽樣概率的倒數(shù)，用下式計算：

樣本個體的抽樣權重（wsi3）值為樣本個體抽樣概率的倒數(shù)?？捎嬎闳缦拢?/p>

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

wsi＝wsi1×wsi2×wsi3＝樣本個體i所在區(qū)縣分層抽樣比倒數(shù)×

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

計算鄉(xiāng)鎮(zhèn)街道（2水平單位）πj，個體被抽中的概率πij及它們的權重，即概率的倒數(shù)，其計算公式為：

ws2＝wsi1×wsi2＝樣本個體i所在區(qū)縣分層抽樣比倒數(shù)×

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

穩(wěn)健標準誤估計（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后的穩(wěn)健標準誤估計：gllamm，robust cluster（psu）

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

方法1：

方法2：

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

表1 2水平隨機截距l(xiāng)ogistic回歸和復雜抽樣2水平隨機截距l(xiāng)ogistic回歸擬合

討 論

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

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

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

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

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

3.繆凡，童峰.復雜抽樣數(shù)據(jù)的logistic回歸分析方法及其應用.中國衛(wèi)生統(tǒng)計，2008，25（6）：577-579.

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

5.楊珉，李曉松主編.醫(yī)學和公共衛(wèi)生研究常用多水平統(tǒng)計模型.北京，北京大學醫(yī)學出版社，2007.

6.賈改珍，閆陽，徐天和，等.多水平模型在大學生預防艾滋病健康教育影響因素分析中的應用.中國衛(wèi)生統(tǒng)計，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.中國疾病預防控制中心公共衛(wèi)生監(jiān)測與信息服務中心（102206）

2.四川大學華西公共衛(wèi)生學院衛(wèi)生統(tǒng)計教研室（610041）

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

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