基于SAS的調(diào)整人群歸因危險(xiǎn)度評(píng)估法在口岸流行病學(xué)調(diào)查中的應(yīng)用研究

裘炯良 鄭劍寧 施惠祥

(1.寧波出入境檢驗(yàn)檢疫局 浙江寧波 315012;2.北侖出入境檢驗(yàn)檢疫局)

1 前言

口岸流行病學(xué)調(diào)查往往涉及口岸人群中風(fēng)險(xiǎn)因素與傳染性疾病的相關(guān)性研究，如口岸居民區(qū)中房前屋后的瓦罐比率與蚊媒的相關(guān)性，以及其導(dǎo)致輸入性登革熱局部暴發(fā)流行的風(fēng)險(xiǎn)性高低評(píng)估等。對(duì)于單個(gè)風(fēng)險(xiǎn)因素與傳染病相關(guān)性研究可以應(yīng)用人群歸因危險(xiǎn)度(Population Attributable Risk，PAR)法來(lái)進(jìn)行綜合評(píng)估［1-3］。然而，要同時(shí)分析多個(gè)風(fēng)險(xiǎn)因素與傳染病的關(guān)聯(lián)性，這個(gè)指標(biāo)則無(wú)法發(fā)揮作用。為此，本研究引入調(diào)整人群歸因危險(xiǎn)度法(PARc)，希望對(duì)出入境檢驗(yàn)檢疫機(jī)構(gòu)科學(xué)估算由于同時(shí)暴露于多種風(fēng)險(xiǎn)因素所引起的傳染病發(fā)生頻率，以及預(yù)測(cè)防制工作中消除這些因素后該病可被減少或消除的有效程度具有一定的指導(dǎo)意義。

2 原理與方法

2.1 人群歸因危險(xiǎn)度

指人群中因暴露于某因素所致疾病或死亡所占人群總發(fā)病或死亡的百分比。它把病原因子的相對(duì)危險(xiǎn)度和人群中暴露者的比例結(jié)合起來(lái)，對(duì)某一暴露因素造成的危害程度進(jìn)行綜合評(píng)價(jià)，是流行病學(xué)調(diào)查分析中一項(xiàng)很有價(jià)值的指標(biāo)。1953年由Levin ML首次提出，故亦稱Levin氏歸因危險(xiǎn)度(Levin’s Attributable Risk)。

Levin氏法測(cè)定PAR的公式：

其中，D表示疾病(disease)，E表示暴露(exposure)，P表示概率(probability)。

未調(diào)整時(shí)的PAR估計(jì)：

式中：Ip：人群某病發(fā)生率

Iu：非暴露者某病發(fā)生率

RR：相對(duì)危險(xiǎn)度(在病例對(duì)照研究中可以O(shè)R估計(jì))

P(E)：總?cè)巳褐斜┞墩叩谋壤?/p>

Walter已證明ln(1-PAR)呈正態(tài)分布，因此可應(yīng)用正態(tài)分布理論來(lái)計(jì)算PAR的可信限和進(jìn)行顯著性檢驗(yàn)。

2.2 調(diào)整人群歸因危險(xiǎn)度

人群歸因危險(xiǎn)度僅適用于單變量分析，如資料中有混雜因素或是多個(gè)因素同時(shí)存在時(shí)，上述公式就不再適用，而須計(jì)算調(diào)整控制一個(gè)或多個(gè)變量后的PAR指標(biāo)。這種僅考慮某一危險(xiǎn)因素而把其余危險(xiǎn)因素均作為干擾因素進(jìn)行校正后所得的PAR稱為調(diào)整人群歸因危險(xiǎn)度，記為PARc。假設(shè)研究總體病例數(shù)為x，對(duì)照數(shù)為n-x;K個(gè)危險(xiǎn)因素(A為研究因素，C為其余校正因子)的不同暴露水平交叉分組構(gòu)成 J層危險(xiǎn)狀態(tài)(j=0，1，2，….j-1)。一個(gè)個(gè)體處于第j層是病例的概率為Ij=xj/nj。

式中IC表示當(dāng)C因素維持原有水平而A因素降為基準(zhǔn)水平時(shí)，該人群的預(yù)期發(fā)病率。

PARc計(jì)算：若令j*表示第j層A因素降為基準(zhǔn)水平，而C因素保持原有水平時(shí)的暴露狀態(tài)。令，式中表示第j層暴露狀態(tài)相對(duì)于第j*層暴露狀態(tài)的相對(duì)危險(xiǎn)度，亦即校正C因子的效應(yīng)后第j層A因素暴露狀態(tài)的相對(duì)危險(xiǎn)度;表示第j層暴露狀態(tài)C因素的相對(duì)危險(xiǎn)度。由此，A因素的調(diào)整人群歸因危險(xiǎn)度可改寫為：

若已知各層在C因素暴露條件下A因素的相對(duì)危險(xiǎn)度和總體中病例在暴露狀態(tài)各層分布的比例，即可算得PARc%值。上式中ρj即為第j層暴露狀態(tài)下病例數(shù)占總病例數(shù)的比例。

以上為針對(duì)研究總體的PARc計(jì)算推導(dǎo)。對(duì)于流行病學(xué)研究的樣本人群，其病例是總體病例的一個(gè)隨機(jī)抽樣，那么，觀察樣本中病例在各層的分布比例ρj＇可作為總體病例在各層的分布比例ρj的估計(jì)值;而對(duì)于發(fā)病率很低的傳染性疾病，觀察資料的校正相對(duì)危險(xiǎn)度可用多因素Logistic回歸模型得出因素的OR值來(lái)估計(jì)。

3 實(shí)例演算

3.1 資料來(lái)源

本文應(yīng)用的數(shù)據(jù)來(lái)自寧波口岸某病流行病學(xué)調(diào)查資料。采用三階段整群隨機(jī)抽樣方法，選擇對(duì)照是以與病例同年齡(10歲為1組)、同性別、同居住地作為1：1頻數(shù)匹配條件，從自然人群中選擇。抽取的143名健康人，其中133例接受流行病學(xué)調(diào)查，應(yīng)答率為93.0%;調(diào)查方式均為面訪;研究的目的是確定某流行病與相關(guān)因素的病因?qū)W聯(lián)系。

3.2 資料均衡性分析

對(duì)頻數(shù)匹配的因素(年齡、性別與居住地)進(jìn)行了均衡性檢驗(yàn)，結(jié)果表明：病例組平均年齡為26歲(四分位間矩：22-33)，略高于對(duì)照組的24歲(四分位間矩：23-31)，將年齡變量轉(zhuǎn)換為以10歲為1組的等級(jí)變量后，經(jīng)χ2檢驗(yàn)，P＞0.05;對(duì)性別及居住地因素的檢驗(yàn)結(jié)果顯示，病例與對(duì)照組之間的分布均勻，差異未達(dá)統(tǒng)計(jì)學(xué)顯著性(P值分別為0.230和 0.262)。

3.3 調(diào)整人群歸因危險(xiǎn)度計(jì)算

表1為應(yīng)用多元Logistic回歸模型法計(jì)算獲得該流行病4種風(fēng)險(xiǎn)因素及風(fēng)險(xiǎn)度估計(jì)值。

表1 寧波口岸某病風(fēng)險(xiǎn)因素及風(fēng)險(xiǎn)度

表2為該病風(fēng)險(xiǎn)因素的調(diào)整人群歸因危險(xiǎn)度值，分別為：風(fēng)險(xiǎn)因素A所致的流行病危害占11.8%，風(fēng)險(xiǎn)因素B占30.8%，風(fēng)險(xiǎn)因素C占15.9%，風(fēng)險(xiǎn)因素D占38.15%。

表2 寧波口岸某流行病風(fēng)險(xiǎn)因素的PARc(%)

4 SAS程序?qū)崿F(xiàn)

筆者自行編制相應(yīng)的SAS程序代碼，如表3所示。在SAS9.2統(tǒng)計(jì)軟件的“PROGRAM EDITOR”窗口中輸入表3的程序代碼，單擊“RUN”按鈕，在“OUTPUT”窗口即可直接查看“數(shù)學(xué)模型統(tǒng)計(jì)學(xué)檢驗(yàn)”、“數(shù)學(xué)模型統(tǒng)計(jì)學(xué)檢驗(yàn)結(jié)果提示”、“寧波口岸某病風(fēng)險(xiǎn)因素及風(fēng)險(xiǎn)度”及“調(diào)整人群歸因危險(xiǎn)度(PARc%)”等4項(xiàng)結(jié)果。表3第04行所列的是口岸流行病學(xué)調(diào)查中獲得的幾組風(fēng)險(xiǎn)因素原始數(shù)據(jù)(為節(jié)省篇幅，僅列出部分)，在實(shí)際工作中可以根據(jù)現(xiàn)場(chǎng)調(diào)查結(jié)果進(jìn)行數(shù)據(jù)替換。

SAS程序首先應(yīng)用多元Logistic回歸模型篩選有統(tǒng)計(jì)學(xué)意義的風(fēng)險(xiǎn)因素共4種，以似然比檢驗(yàn)方法計(jì)算整個(gè)回歸模型的顯著性水平，P值=0.0108;SAS進(jìn)行程序自動(dòng)判定，給出提示“數(shù)學(xué)模型檢驗(yàn)顯示有統(tǒng)計(jì)學(xué)意義!”;然后，SAS程序自動(dòng)運(yùn)算，給出該4種風(fēng)險(xiǎn)因素在校正其他因素前提下的權(quán)重系數(shù)估計(jì)值、標(biāo)準(zhǔn)誤、Wald卡方值、顯著性檢驗(yàn)P值(表1)以及OR值和其95%CI;最后，SAS程序自動(dòng)對(duì)每個(gè)風(fēng)險(xiǎn)因素進(jìn)行卡方分析，給出每個(gè)風(fēng)險(xiǎn)因素不同水平的病例和對(duì)照數(shù)量，同時(shí)計(jì)算相應(yīng)的OR值和95%CI，并以公式5運(yùn)算得出每個(gè)風(fēng)險(xiǎn)因素的調(diào)整人群歸因危險(xiǎn)度(表2)。為節(jié)省篇幅，表3僅給出計(jì)算其中一個(gè)風(fēng)險(xiǎn)因素的PARc%所需的SAS程序源代碼，其他風(fēng)險(xiǎn)因素的PARc%計(jì)算程序同理可得。

表3 實(shí)現(xiàn)調(diào)整人群歸因危險(xiǎn)度計(jì)算的SAS程序代碼

53 proc freq data=par;54-56 table factor4*case/chisq;run;data list_5;57 set list_1 list_2 list_3 list_4;58 keep factor1 factor2 factor3 factor4 case Frequency sum_1;59 if case=1;/sle=0.5 sls=0.5 noint; 60 if factor1^=.or factor2^=.or factor3^=.or factor4^=.;12 proc print data=test NOOBS; 61-62 sum_1=103;proc sort;13 title＇數(shù)學(xué)模型統(tǒng)計(jì)學(xué)檢驗(yàn)＇; 63-65 by factor1;data list_6;set list_1 list_2 list_3 list_4;14-15 ods listing;data t; 66 keep factor1 frequency control;16 set test; 67 if case=0;17 if test= ＇Likelihood Ratio＇; 68 if factor1^= .or factor2^= .or factor3^= .or factor4^= .;18-19 title＇數(shù)學(xué)模型統(tǒng)計(jì)學(xué)檢驗(yàn)結(jié)果提示＇;file print; 69-70 data list_7;set list_6;20 if Prob_ChiSq ＜0.05 then put//+30 ＇結(jié)果提示：數(shù)學(xué)模型檢驗(yàn)顯示有統(tǒng)計(jì)學(xué)意義!＇;71-73 rename frequency=control;proc sort;by factor1;21 else put//+30＇結(jié)果提示：數(shù)學(xué)模型檢驗(yàn)顯示無(wú)統(tǒng)計(jì)學(xué)意義!＇;74-75 data list_all;merge list_5 list_7;22-23 data OR_EST;merge EST OR; 76-80 by factor1;data calcu_1;set list_all;if factor1^=.;class=put(factor1，$1.);24 if effect= ＇factor1＇or effect= ＇factor2＇or effect= ＇factor3＇or effect= ＇factor4 ＇;81-84 data calcu_11;set or_est_1;rename ClassVal0=class;if variable= ＇factor1＇;25 DROP EFFECT DF; 85 data calcu_111;26 PROC PRINT NOOBS; 86 merge calcu_1 calcu_11;27 title＇寧波口岸某病風(fēng)險(xiǎn)因素及風(fēng)險(xiǎn)度＇; 87 by class;28 RUN;ods listing close; 88 if class= ＇1＇then oddsRatioEst=1;29 ods output OddsRatios=OR_1; 89 parc=Frequency/sum_1/oddsRatioEst;30 ods output ParameterEstimates=est_1; 90 ods output Summary=summy;31 ods output GlobalTests=test_1; 91 proc means sum data=calcu_111;32 proc logistic descending data=par; 92-94 var parc;data a;set summy;33 class factor1(ref=first)factor2(ref=first)factor3(ref=first)factor4(ref=first)/param=ref;95-96 parc_1=abs(1-parc_Sum)*100;data calcu_1111;34 model case=sex address age factor1 factor2 factor3 factor4/sle=0.5 sls=0.5 noint;97 merge calcu_111 a;35 data OR_EST_1; 98 keep factor1 frequency control oddsRatioEst parc_1;36 merge EST_1 OR_1; 99 data calcu_11111;37 if effect^= ＇sex＇and effect^= ＇address＇and effect^= ＇age＇; 100 set calcu_1111;38 DROP EFFECT DF; 101 rename frequency=case;39 ods listing close; 102-103 ods listing;PROC PRINT noobs;40 ods output CrossTabFreqs=list_1; 104 itle＇寧波口岸某病風(fēng)險(xiǎn)因素A的調(diào)整人群歸因危險(xiǎn)度(PARc%)＇;41 proc freq data=par; 105-106 RUN;quit;

5 討論

人群歸因危險(xiǎn)度的計(jì)算相對(duì)比較簡(jiǎn)單，但在自然現(xiàn)象中，往往有混雜因素或是多個(gè)因素同時(shí)存在，則單因素人群歸因危險(xiǎn)度計(jì)算公式就不再適用，而須計(jì)算調(diào)整控制一個(gè)或多個(gè)變量后的PAR指標(biāo)［4］。國(guó)外有學(xué)者曾提出用Mantel-Haenszel法、非條件法、加權(quán)合并法及基于回歸模型的調(diào)整估計(jì)法等方法來(lái)調(diào)整多個(gè)干擾因子的效應(yīng)。但筆者認(rèn)為，以基于回歸模型的調(diào)整估計(jì)法為優(yōu)，主要考慮以下幾點(diǎn)因素：①基于回歸模型的調(diào)整估計(jì)法允許控制調(diào)整較多的混雜和效應(yīng)修正因子;②混雜項(xiàng)、交互項(xiàng)均較易于直接加入模型方程中;③參數(shù)估計(jì)和假設(shè)檢驗(yàn)靈活方便。

本研究正是采用Bruzzi提出的基于回歸模型的調(diào)整估計(jì)法，其原理不是直接擬合PAR模型，而是把基于回歸模型的調(diào)整的RR估計(jì)加入到PAR估計(jì)中去。Bruzzi認(rèn)為，如果病例對(duì)照資料的病例是隨機(jī)樣本，用Logistic回歸模型估計(jì)諸危險(xiǎn)因素與傳染性疾病的相對(duì)危險(xiǎn)度后，即可經(jīng)病例的暴露分布求出PAR。本研究的病例為寧波口岸某病全年所有的病例人群，故與Bruzzi提出的確定病例人群的條件并不矛盾，因此亦可直接應(yīng)用。Bruzzi法的優(yōu)點(diǎn)是僅需知道病例的暴露分布而不需人群暴露狀況，即可估計(jì)危險(xiǎn)因素不同水平的PAR值。

由于整個(gè)運(yùn)算過(guò)程比較繁雜，計(jì)算工作量較大，且靠手工計(jì)算容易出錯(cuò)。為此，本研究應(yīng)用國(guó)際權(quán)威的統(tǒng)計(jì)分析標(biāo)準(zhǔn)軟件——SAS分析平臺(tái)進(jìn)行編程［5］，還在程序中設(shè)定了一些自動(dòng)判斷過(guò)程，實(shí)現(xiàn)了多因素影響的口岸流行病風(fēng)險(xiǎn)因素調(diào)整人群歸因危險(xiǎn)度的全自動(dòng)計(jì)算機(jī)運(yùn)算，這在國(guó)內(nèi)相關(guān)研究中尚屬首次。所得結(jié)果經(jīng)過(guò)了手工計(jì)算及SPSS軟件半自動(dòng)運(yùn)算的雙重驗(yàn)證，證明準(zhǔn)確可靠，為口岸流行病發(fā)生與控制效果的評(píng)價(jià)提供了一種簡(jiǎn)便快捷且客觀有效的評(píng)估工具。

［1］Rebecca J Guy，Handan Wand，David P Wilson，et al.Using population attributable risk to choose HIV prevention strategies in men who have sex with men ［J］.BMC Public Health，2011，11：247.

［2］Giovanni de Simone，Jorge R.Kizer，Marcello Chinali，et al.Normalization for body size and population-attributable risk of left ventricular hypertrophy［J］.Am J Hypertens，2005，18：191-196.

［3］Wendy Y Chen，Bernard Rosner，Susan E Hankinson，et al.Moderate Alcohol Consumption During Adult Life，Drinking Patterns，and Breast Cancer Risk［J］.JAMA，2011，306(17)：1884-1890.

［4］沈月平，武俊青，高爾生.隱函數(shù)Delta法調(diào)整人群歸因危險(xiǎn)度可信區(qū)間的估計(jì)及其應(yīng)用［J］.中國(guó)衛(wèi)生統(tǒng)計(jì)，2004，21(3)：142-145.

［5］高惠璇，譯.SAS系統(tǒng)SAS/STAT軟件、BASE SAS軟件使用手冊(cè)［M］.北京：中國(guó)統(tǒng)計(jì)出版社，1997，1-30.