東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計教研室(210009) 田 野 閔 捷 劉 沛 王燦楠

在進行污染物的長期膳食暴露評估時，合理地估計每人每天長期攝入某種物質(zhì)平均量(usual intake，簡稱日常攝入量)往往是風(fēng)險評估過程中的重點。日常攝入量的計算涉及該物質(zhì)在各種食物中的濃度和人群消費各種食物的量，前者可用各種食物污染物檢測的平均濃度估計，不考慮其變異;后者需要采集人群長期消費各種食物的信息，實際工作中是難以實現(xiàn)的，如何利用短期2～7天的24小時膳食回顧調(diào)查數(shù)據(jù)來估計人群長期暴露量就顯得尤為重要。許多專家及機構(gòu)提出了估計及改進的方法〔1-3〕，目前大家普遍接受的是貝塔二項-正態(tài)分布模型(Beta binomial-normal，BBN)，該方法是2005年由De Boer和Van der Voe等人提出的〔4〕。其基本思想是考慮調(diào)查人群中，個體間變異可近似反映人群的長期飲食情況，個體內(nèi)變異則被認(rèn)為僅反映調(diào)查個體在調(diào)查的n天內(nèi)短期膳食結(jié)構(gòu)的波動，采用合適的統(tǒng)計模型去除個體內(nèi)變異，保留個體間變異的暴露量，近似地作為人群的長期暴露量。本文利用BBN模型對江蘇省各年齡、性別組人群氰戊菊酯的長期膳食安全風(fēng)險水平進行評估，以期為農(nóng)藥膳食長期暴露風(fēng)險評估提供科學(xué)依據(jù)。

BBN模型原理

BBN模型包括攝入頻率和陽性暴露量兩部分。其中，貝塔二項分布被用來構(gòu)建攝入頻率模型，對數(shù)正態(tài)分布被用來構(gòu)建陽性暴露量模型。

1.攝入頻率模型

假設(shè)n和npos分別代表每個調(diào)查對象的總調(diào)查天數(shù)和攝入食物中含有該污染物的天數(shù)(即陽性攝入天數(shù))。由于每人每天攝入的食物中可能包含污染物，也可能不包含污染物，因此陽性攝入天數(shù)npos服從陽性攝入概率為p，重復(fù)數(shù)為n的二項分布;其中，陽性攝入概率p的取值在0到1之間，服從貝塔分布;那么，陽性攝入天數(shù)就服從貝塔二項分布:

貝塔分布的均數(shù)和方差分別為

令π=α/(α+β)，φ=1/(α+β+1)，則該分布的均數(shù)和方差可表示為:

通過最大似然估計法擬合貝塔二項分布模型。可以得到π和φ的估計值π和φ，從而反推得到α和β的估查對象在某一調(diào)查日食用某種事物的陽性攝入頻率分布即可用貝塔分布

2.陽性暴露量模型

長期膳食暴露評估主要反應(yīng)人群中某污染物日常攝入情況，假設(shè)陽性暴露量即調(diào)查對象每人每天膳食污染物暴露量不為零，經(jīng)過對數(shù)變換后近似服從正態(tài)分布，變換后的陽性暴露量，通過最大似然估計法擬合線性隨機效應(yīng)模型(linear random effects model):

其中，transf(yposij)表示經(jīng)對數(shù)變換后的陽性暴露量，ci和uij表示人間效應(yīng)和人內(nèi)效應(yīng)，其二者分布分別滿足正態(tài)分布 N(0，σ2B)和N(0，σ2w)。通過最大似然估計法可以得到三個參數(shù)的估計值m，分別代表個體平均值、人間效應(yīng)估計值及人內(nèi)效應(yīng)估計值。

3.計算長期暴露量分布

長期暴露量是通過對陽性暴露量分布和攝入頻率分布進行Monte Carlo抽樣，然后進行整合獲得。其具體過程為:

首先通過模擬獲得對數(shù)轉(zhuǎn)換后的陽性暴露量m+zi，其中zi從正態(tài)分布N(0)中隨機抽取10萬次，然后反對數(shù)變換獲得陽性暴露量:

實例分析

1.資料情況主要涉及的數(shù)據(jù)有膳食調(diào)查數(shù)據(jù)、污染物常規(guī)監(jiān)測數(shù)據(jù)、研究人群人口學(xué)基本情況數(shù)據(jù)。膳食調(diào)查及人口學(xué)數(shù)據(jù):來源于2002年全國居民營養(yǎng)與健康狀況調(diào)查中的江蘇省數(shù)據(jù);污染物監(jiān)測數(shù)據(jù):來源于2000～2006年全國14個省/地區(qū)食品污染物監(jiān)測網(wǎng)數(shù)據(jù)及2005～2006年海關(guān)出口農(nóng)產(chǎn)品監(jiān)測數(shù)據(jù)中的江蘇省數(shù)據(jù)。本研究資料中氰戊菊酯檢出率為3.10%，按照WHO關(guān)于污染物未檢出值處理的原則，本文中未檢出值均采用0替換〔5〕。目前，F(xiàn)AO/WHO食品添加劑聯(lián)合專家委員會(JECFA)制定的食品中每日允許攝入氰戊菊酯量(ADI)為2μg/kg bw/day。將模型結(jié)果與之進行比較，評價人群的氰戊菊酯攝入水平。

2.結(jié)果與分析

應(yīng)用SAS9.1統(tǒng)計軟件進行數(shù)據(jù)整理和膳食暴露量的計算，結(jié)果如下。

(1)BBN模型參數(shù)

利用SAS的Proc nlmixed過程的最大似然估計法迭代獲得貝塔分布參數(shù)和的估計值。在對陽性暴露量進行方差分解時，首先將陽性暴露量進行對數(shù)變換，使其近似正態(tài)分布，然后利用SAS的Proc mixed過程將殘差分解為個體間效應(yīng)和個體內(nèi)效應(yīng)，通過最大似然估計法獲得個體平均值m、人間效應(yīng)估計值及人內(nèi)效應(yīng)估計值。計算的氰戊菊酯膳食暴露BBN模型參數(shù)見下表。

表1 氰戊菊酯膳食暴露BBN評估模型參數(shù)

(2)氰戊菊酯長期膳食暴露評估結(jié)果

如表2可見，氰戊菊酯的長期暴露量在平均值及各百分位數(shù)均低于其每日允許攝入量ADI(2 μg/kg bw/day)，說明氰戊菊酯殘留在長期膳食暴露評估中是相對安全的。

表2 不同人群氰戊菊酯長期膳食暴露量(單位μg/kg bw/day)

(3)氰戊菊酯膳食暴露量的不確定性分析

為了解人群氰戊菊酯暴露變異的各百分位數(shù)的可信區(qū)間，采用Bootstrap方法有放回的從原樣本中進行抽樣，并對抽取的每一個Bootstrap樣本進行Monte Carlo模擬以計算暴露量百分位數(shù)，經(jīng)過n次Bootstrap抽樣計算百分位數(shù)分布及可信區(qū)間，代表百分位數(shù)的不確定性。本文取200次Bootstrap抽樣進行不確定性分析。

由表3可見，各暴露水平的暴露量均位于95%可信區(qū)間內(nèi)，說明暴露量百分位數(shù)能夠達(dá)到穩(wěn)定。由于本次研究采用的是全比例抽樣，其不確定性主要存在于攝入頻率p～Beta(α，β)和人間效應(yīng)zi～N(0)的抽樣中，同時Monte Carlo模擬較好的量化了變異性，故在可信區(qū)間內(nèi)變異較小。

討 論

1.長期膳食暴露評估模型

評價污染物長期膳食暴露有諸多不同的數(shù)學(xué)模型，其共同思想是考慮在大規(guī)模人群短期橫斷面膳食調(diào)查情況下，數(shù)據(jù)資料中變異主要包括個體間變異和個體內(nèi)變異。不同個體間的變異可以近似反映人群的長期膳食狀況，而個體內(nèi)變異則被看做僅反映調(diào)查個體在橫斷面調(diào)查中短期的膳食結(jié)構(gòu)波動，因此，通過采用合適的數(shù)學(xué)模型來去除個體內(nèi)變異，而保留個體間變異，從而達(dá)到用污染物短期的膳食暴露量近似代表長期的膳食暴露量的目的。本研究在計算長期膳食暴露評估時采用歐盟MCRA(Monte Carlo Risk Assessment)推薦的BBN模型。將長期的污染物膳食陽性攝入頻率用短期陽性攝入頻率的貝塔二項分布模型來表示;同時，用正態(tài)分布構(gòu)建長期陽性暴露量模型。從分布中抽樣，進行Monte Carlo模擬獲得污染物的長期暴露量。利用BBN模型進行長期膳食暴露評估的優(yōu)點是將個體作為研究對象，通過統(tǒng)計方法進行模擬，使結(jié)果更符合實際并獲得人群的暴露量分布，其對膳食調(diào)查數(shù)據(jù)及污染物殘留數(shù)據(jù)的信息利用程度遠(yuǎn)大于點評估模型。國外一些研究進行了模型的比較也證明了這一結(jié)論〔6〕。當(dāng)然，此種暴露評估方法也有其局限性，如模擬次數(shù)足夠多才能使最終結(jié)果穩(wěn)定及對數(shù)據(jù)質(zhì)量要求較高等。另外該模型適用于污染物檢出率不高的暴露評估，如果檢出率高，例如鉛、鎘等的污染物，則BBN中可不考慮Beta分布。

表3 不同人群氰戊菊酯長期膳食暴露量不確定性分析(單位μg/kg bw/day)

2.關(guān)于陽性暴露量的正態(tài)性轉(zhuǎn)換

BBN模型中，需要利用正態(tài)分布構(gòu)建陽性攝入量模型。人群的陽性攝入量的原始數(shù)據(jù)往往不滿足正態(tài)性，常常需要用變量變換的方法將偏態(tài)資料正態(tài)化。在BBN模型中常用的變量變換有兩種:對數(shù)變換和冪變換〔7〕。經(jīng)過變換后的變量，一般需對其分布進行擬合優(yōu)度檢驗，常用的檢驗有Anderson-Darling檢驗、Kolmogorov-Smirnov檢驗、Chisquare檢驗等。本研究中，陽性暴露量的原始數(shù)據(jù)表現(xiàn)為正偏態(tài)的基本特征，本研究采用對數(shù)變換對陽性攝入量的原始數(shù)據(jù)進行轉(zhuǎn)換。

3.Bootstrap抽樣次數(shù)的確定

關(guān)于Bootstrap抽樣次數(shù)的確定，一般根據(jù)研究目的的不同取不同的值，如果用來計算均值、方差、標(biāo)準(zhǔn)誤等統(tǒng)計量，200次抽樣即可，但如果計算極端百分位數(shù)如99.9th的可信區(qū)間，Bootstrap的模擬次數(shù)應(yīng)當(dāng)增加到至少500次〔8〕。本研究Monte Carlo模擬的10萬個人天數(shù)作為一個Bootstrap樣本含量，根據(jù)需要估計的百分位數(shù)情況，將Bootstrap樣本數(shù)定為200次、500次、1000次分別進行模擬，比較抽樣結(jié)果發(fā)現(xiàn)可信區(qū)間寬度沒有明顯差異，最終確定本研究中Bootstrap模擬次數(shù)為200次比較合適。需要說明的是，對不確定度更多的應(yīng)從專業(yè)上加以評價，如果風(fēng)險管理者認(rèn)為數(shù)據(jù)的不確定度過大，可以拒絕概率估計的結(jié)論，或者采取一定措施(如提高數(shù)據(jù)質(zhì)量)來減少不確定度〔9〕。

1.Nusser SM，Carriquiry AL，Dodd KW，et al.A semi-parametric transformation approach to estimating usual intake daily intake distributions.J Am Stat Assoc，1996，91:1440-1449.

2.Slob W.Modeling Long-term exposure of the whole population to chemicals in food.RISK ANAL，1993，13(5):525-530.

3.Slob W.Probabilistic dietary exposure assessment taking into account variability in both amount and frequency of consumption.Food and Chemical Toxicology，2006，44(7):933-951.

4.De Boer WJ，van der Voet H.MCRA，Release 5，a web-based program for Monte Carlo Risk Assessment.Biometris and RIKILT，Wageningen University and Research Centre，2006.

5.王緒卿，吳永寧.食品污染監(jiān)測低水平數(shù)據(jù)處理問題.中華預(yù)防醫(yī)學(xué)雜志，2002:837-852.

6.Boon PE，Bonthuis M.Comparison of different exposure assessment methods to estimate the long-term dietary exposure to dioxins and ochratoxin A.Food and Chemical Toxicology，2011，49(9):1979-1988.

7.De Boer WJ，van der Voet H.Comparison of two models for the estimation of usual intake addressing zero consumption and non-normality.Food Additives and Contaminants，2009，26(11):1433-1449.

8.Cullen AC，F(xiàn)rey HC.Probabilistic Techniques in Exposure Assessment:A Handbook for Dealing with Variability and Uncertainty in Models and Inputs.Society for Risk Analysis，1999.

9.Thompson KM.Variability and uncertainty meet risk management and risk communication.Risk Anal，2002，22:647-654.