彭雯潔 崔 靖 張 婷 王青青 張 超 袁 苗 周小林△ 余紅梅,4△

【提 要】 目的 比較BIC估計法與MCMC近似法兩種后驗概率法在貝葉斯基準劑量估計中的穩(wěn)健性,并為山西省洪洞縣兒童羥基代謝物可接受劑量的制定提供參考建議。方法 首先介紹基于BIC估計法和MCMC近似法計算后驗權(quán)重的原理,模擬研究選用Integrated Risk Information System數(shù)據(jù)庫中不同劑量-反應數(shù)據(jù)集共30個,分析比較兩種方法的優(yōu)劣,并在實例研究中采用權(quán)重法進行數(shù)據(jù)整合。結(jié)果 模擬研究結(jié)果顯示在所研究的30個數(shù)據(jù)集中BIC估計法在BMR為0.01時有4個數(shù)據(jù)集出現(xiàn)BMDL預測失敗的情況,在BMR為0.001時有1個數(shù)據(jù)集出現(xiàn)BMD預測失敗的情況,以及6個數(shù)據(jù)集出現(xiàn)BMDL預測失敗的情況。MCMC近似法計算的BMD/BMDL在每一種模型都有70%以上的數(shù)據(jù)集高于BIC估計法得到的BMD/BMDL。實例分析表明符合洪洞縣兒童體內(nèi)羥基代謝物劑量-反應關(guān)系的模型有l(wèi)inear(P=0.13,β=14.3%)、logistic(P=0.06,β=9.5%)、Weibull(P=0.14,β=10.6%)、multistage(P=0.15,β=31.1%)、Hill(P=0.21,β=34.6%)。在BMR為0.001的情況下,洪洞縣兒童體內(nèi)八種羥基代謝物(2-OHN、1-OHN、9-OHF、2-OHF、2-OHphe、1-OHphe、1-OHBaP、3-OHBaP)的可接受劑量(μmol/mol)依次為0.577 μmol/mol、1.546 μmol/mol、8.135 μmol/mol、0.359 μmol/mol、0.120 μmol/mol、0.098 μmol/mol、0.044 μmol/mol、0.003μmol/mol。結(jié)論 MCMC近似法在BMD估計中具有較好的穩(wěn)定性和魯棒性。

基準劑量(benchmark dose,BMD)是計算在有害物質(zhì)背景值的基礎(chǔ)上引起基準反應(benchmark response,BMR)不良健康效應的劑量,是危害表征的指標之一[1]。同時BMD是有毒有害化學物質(zhì)在某一特定環(huán)境中可接受的最高界值,若環(huán)境中有毒有害化學物質(zhì)超過這一界值就會造成人群健康壽命年的降低。BMDL(the lower confidence limit of BMD)是在規(guī)定了置信區(qū)間之后BMD的下限,它為環(huán)境工作者在制定界值范圍時起到了規(guī)范作用,是經(jīng)濟衛(wèi)生領(lǐng)域提高效益-收益的指標之一。1984年美國環(huán)境保護署(environmental protection agency,EPA)提出了最優(yōu)模型法,即在提出的所有模型中選擇赤池信息量(akaike information criterion,AIC)最小的模型計算BMD[2]。隨著計算機的發(fā)展和貝葉斯算法的出現(xiàn),1986年Crump提出了貝葉斯基準劑量[3](bayesian benchmark dose,BBMD),BBMD是利用貝葉斯框架提供了一種通過模型參數(shù)的先驗分布來整合先驗信息的方法。這對于提高低質(zhì)量數(shù)據(jù)的劑量-反應建模的可靠性具有很大的潛力。同時BBMD也克服了傳統(tǒng)方法對試驗時間、試驗樣本和試驗劑量的高要求性,從而計算出在特定要求下的BMD[4]。但有時在處理低暴露弱相關(guān)數(shù)據(jù)時,單一模型并不能完全闡明數(shù)據(jù)劑量-反應特征,若要選取最優(yōu)模型,就會不可避免地放棄其他模型提供信息的可能,從而導致估計模型的不確定性[5]。為充分考慮模型的不確定性,2018年Khao將模型平均法應用在貝葉斯基準劑量中,提出了貝葉斯模型平均法(bayesian model averaging,BMA)的思想[5]。BMA是一種在BBMD方法的基礎(chǔ)上以模型的后驗概率為權(quán)重,將備選模型的不確定性考慮在內(nèi)的統(tǒng)計學方法。它能夠綜合考慮不同的備選模型的權(quán)重,使分析更具有科學性和完整性[6]。BMA得以應用的關(guān)鍵在于計算每個模型所占的權(quán)重,目前的研究方法有基于貝葉斯信息準則(Bayesian information criterion,BIC)和馬爾可夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)兩種[7]。本文主要目的在于比較這兩種方法的優(yōu)劣,并計算山西省臨汾市洪洞縣兒童羥基內(nèi)暴露的可接受范圍,為環(huán)境管理部門制定相關(guān)標準提供理論依據(jù)。

原理與方法

1.貝葉斯基準劑量估計模型

基準劑量估計首先要確定所研究的劑量-反應數(shù)據(jù)類型和基準反應,對于不同的數(shù)據(jù)類型,EPA介紹了不同種類的模型,本次著重于二分類數(shù)據(jù)的研究,主要涉及的模型包括以下八種:

linear模型:f(d)=a+(1-a)×[1-exp(-b×d)]

Probit模型:f(d)=Φ(a+b×d)

Weibull模型:f(d)=a+(1-a)×[1-exp(-b×dg)]

multistage模型:f(d)=a+(1-a)×[1-exp(-b×d-c×d2)]

Log-Probit模型:f(d)=a+(1-a)×Φ[b+g×log(d)]

其中:d為劑量(dose);a,b為單個模型的斜率和截距參數(shù):0≤a≤1,b≥0;g為背景參數(shù)。

對定性數(shù)據(jù)基準反應的確定,一般選取額外風險或附加風險反應類型。由于額外風險在應用中不容易受到背景值的影響,表現(xiàn)更加穩(wěn)定,所以在毒理學研究中一般選取BMR=0.01的額外風險計算,但應用于人群數(shù)據(jù)時BMR則取值更低,一般為0.001。

BMRadd=f(BMD)-f(0)

其中:f為二分類劑量-反應模型;f(BMD)為當前值;f(0)為背景值。

2.貝葉斯模型平均法

BMA就是使用BIC和MCMC計算后驗分布,得到模型后驗預測P值和后驗權(quán)重。利用模型后驗預測P值篩選出可能的多個模型,以后驗權(quán)重對單一模型加權(quán)整合,最后得出所求的BMD。

BICK=-2LK+pKlog(N)

3.權(quán)重法

多環(huán)芳烴羥基代謝物屬于低暴露弱相關(guān)數(shù)據(jù),雖能夠造成人體血液系統(tǒng)的損害,但目前沒有單一代謝物與血液指標一一對應的關(guān)系,所以本次采用權(quán)重法進行數(shù)據(jù)整合[8]

統(tǒng)計分析均采用R 3.6.3,基準劑量估計采用BMD和BBMD軟件進行計算。

4.模型的驗證與評價

研究分為模擬研究和實例研究兩個部分,模擬研究分別計算出兩種方法的BMD以及BMDL值,設置BMR為0.01和0.001。本次研究借鑒Shao[5]等人的評價指標觀察兩種方法中BMD/BMDL,來評價基準劑量估計方法的可靠性,比值越低表明模型的參考區(qū)間越精確,說明其方法的精確性和可靠性越高。同時觀察兩種方法中數(shù)據(jù)集預測失敗的例數(shù),失敗例數(shù)越低則魯棒性越高。

模擬研究

1.數(shù)據(jù)來源

選用Integrated Risk Information System數(shù)據(jù)庫(https://www.epa.gov/iris)中不同劑量-反應數(shù)據(jù)集共30個,分別采用基于BIC估計法與MCMC近似法計算BMD并做單一模型的比較。

2.兩種模型評價

使用BIC估計法與MCMC近似法分別計算30個數(shù)據(jù)集中的BMD和BMDL。由表1可知,30個數(shù)據(jù)集中BIC在BMR為0.01時出現(xiàn)4個數(shù)據(jù)集BMDL預測失敗的情況,分別為Weibull、L-logistic和Hill模型。在BMR為0.001時出現(xiàn)1個數(shù)據(jù)集BMD預測失敗的情況,以及6個數(shù)據(jù)集BMDL預測失敗的情況。MCMC近似法在兩種BMR的情況下均未出現(xiàn)預測失敗的情況。剔除預測失敗的數(shù)據(jù)集之后,計算每個模型在兩種方法中得到的BMD/BMDL,表2顯示每一種模型都有70%以上的BIC估計法得到的BMD/BMDL高于MCMC近似法得到的BMD/BMDL(原始結(jié)果見https://github.com/motobndr/penny)。

表1 MCMC與BIC法對30個數(shù)據(jù)集計算BMD與BMDL的成功個數(shù)

表2 MCMC法對預測成功的數(shù)據(jù)集計算BMD/BMDL的合格率

實例分析

1.數(shù)據(jù)來源

選取2016-2018年山西省臨汾市洪洞縣焦化廠周圍居住范圍在2km以內(nèi)的5～12歲的全部兒童。最終收集到408名合格兒童的血常規(guī)和尿中羥基代謝物(2-OHN、1-OHN、9-OHF、2-OHF、2-OHphe、1-OHphe、1-OHBaP、3-OHBaP)。

2.分析結(jié)果

對所測得的八種羥基代謝物濃度進行正態(tài)性檢驗,得出八種羥基代謝產(chǎn)物濃度均不符合正態(tài)分布(P<0.001),以中位數(shù)所占百分比代表每種代謝物濃度對血液系統(tǒng)影響,其權(quán)重依次為:0.053、0.142、0.747、0.033、0.011、0.009、0.004、0.0003(表3)。按權(quán)重計算408名兒童的綜合羥基代謝物濃度,得出泛化的劑量-反應關(guān)系并按總?cè)藬?shù)等分為10組(表4)。

表3 洪洞縣兒童尿中不同羥基代謝物在綜合劑量中的權(quán)重

表4 洪洞縣兒童在不同劑量組下血常規(guī)異常人數(shù)

選取EPA所提供的所有二分類模型,設置MCMC中迭代次數(shù)為30000,選取一條馬爾可夫鏈條,以50%的樣本量作為預熱。計算模型符合的后驗預測值(W),隨機種子數(shù)設置為82063。以0.05

圖1 五種入選模型的劑量-反應關(guān)系

表5 不同模型的后驗預測值和對應權(quán)重

為了得到更加安全的BMD,設置BMR=0.001,計算單一模型和平均法模型的BMD與BMDL(表6)。按權(quán)重分別計算各個羥基代謝物的內(nèi)暴露可接受劑量(表7)。

表6 單一模型和BMA計算綜合內(nèi)暴露可接受劑量(μmol/mol)

表7 八種羥基代謝物內(nèi)暴露可接受劑量(μmol/mol)

討 論

BMD估計目前有參數(shù)法、非參數(shù)法、半?yún)?shù)法和模型平均法四種類型。其中參數(shù)法分為頻率論參數(shù)法和貝葉斯參數(shù)法,頻率論參數(shù)法最為簡單而且效率高,但在多個模型的比較中若AIC都比較接近時就會出現(xiàn)選擇最佳模型而忽略其他模型可能性問題,即模型的不確定性的問題[9]。非參數(shù)法和半?yún)?shù)法雖然不拘泥于模型的選擇而且可以更精準地擬合劑量-反應關(guān)系[10],但非參數(shù)建模相關(guān)的狄利克雷先驗法以及半?yún)?shù)建模相關(guān)的概率核和三次B樣條,其計算方式復雜,檢驗效率低,而且有較高的過擬合風險,所以非參數(shù)法和半?yún)?shù)法一直作為參數(shù)法的補充方法,最關(guān)鍵的是目前沒有一種成型非參數(shù)或者半?yún)?shù)算法可以廣泛地為環(huán)境工作者所用。BMA是貝葉斯參數(shù)法的衍生方法,其采用后驗概率計算參數(shù)模型的權(quán)重,既解決了單個參數(shù)模型不確定性的問題,同時也保留了效率高的優(yōu)點。更重要的一點是BMA在理解和計算上要優(yōu)于非參數(shù)和半?yún)?shù)模型,對于環(huán)境專業(yè)的工作者具有很強的實用性。

利用BIC估計法對30個數(shù)據(jù)集進行計算時,在BMR為0.01時出現(xiàn)4個數(shù)據(jù)集的BMDL估計失敗的情況。在BMR為0.001時出現(xiàn)1個數(shù)據(jù)集的BMD估計失敗的情況以及6個數(shù)據(jù)集的BMDL估計失敗的情況,即隨著BMR的降低,BIC估計法會出現(xiàn)更多無法估計的情況,同時也會降低BIC估計法在低暴露弱相關(guān)數(shù)據(jù)中的應用范圍。目前EPA提出[8]BMR的最小取值為0.001,本次研究表明在不同劑量-反應關(guān)系下MCMC近似法并沒有出現(xiàn)BMD或BMDL值估計失敗的情況,這與Shao的研究一致[11]。BMD/BMDL比值經(jīng)常用于基準劑量模型估計效果的評價,通常BMD/BMDL越低,表示模型的穩(wěn)定性和可靠性越好[5]。模擬結(jié)果顯示,MCMC得到的BMD/BMDL在總體上要低于BIC得到的BMD/BMDL,所以MCMC在實際應用中具有較高的參考價值。

洪洞縣隸屬于山西省臨汾市,地處山西省南部,臨汾盆地北端。其特點是工業(yè)產(chǎn)業(yè)聚集,特別是焦化業(yè)。2018年央視新聞報道,洪洞縣有以三維集團為首的焦化廠常年違規(guī)排放污水、廢氣造成嚴重的環(huán)境污染,也對工廠周圍的村民造成了一定程度的危害[12]。相關(guān)研究表明兒童可能對這些暴露更為敏感。在空氣污染濃度相同的情況下,兒童較弱的新陳代謝能力會導致其體內(nèi)羥基代謝物水平更高[13]。多環(huán)芳烴不僅對兒童智力發(fā)育、行為和免疫功能有影響,而且可能與兒童情緒煩躁、易怒、記憶力降低以及注意力缺乏有關(guān)[14-15]。因此,了解污染區(qū)兒童羥基的內(nèi)暴露風險十分重要。本次研究采用基于五種函數(shù)建立MCMC貝葉斯平均法模型估算羥基代謝物內(nèi)暴露對應的BMD,得出在BMR為0.001時八種羥基代謝物,以期為當?shù)丨h(huán)境部門制定相關(guān)安全劑量提供參考。