BIC與MCMC基準(zhǔn)劑量估計(jì)模型的比較研究*

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

【提 要】 目的 比較BIC估計(jì)法與MCMC近似法兩種后驗(yàn)概率法在貝葉斯基準(zhǔn)劑量估計(jì)中的穩(wěn)健性,并為山西省洪洞縣兒童羥基代謝物可接受劑量的制定提供參考建議。方法 首先介紹基于BIC估計(jì)法和MCMC近似法計(jì)算后驗(yàn)權(quán)重的原理,模擬研究選用Integrated Risk Information System數(shù)據(jù)庫(kù)中不同劑量-反應(yīng)數(shù)據(jù)集共30個(gè),分析比較兩種方法的優(yōu)劣,并在實(shí)例研究中采用權(quán)重法進(jìn)行數(shù)據(jù)整合。結(jié)果 模擬研究結(jié)果顯示在所研究的30個(gè)數(shù)據(jù)集中BIC估計(jì)法在BMR為0.01時(shí)有4個(gè)數(shù)據(jù)集出現(xiàn)BMDL預(yù)測(cè)失敗的情況,在BMR為0.001時(shí)有1個(gè)數(shù)據(jù)集出現(xiàn)BMD預(yù)測(cè)失敗的情況,以及6個(gè)數(shù)據(jù)集出現(xiàn)BMDL預(yù)測(cè)失敗的情況。MCMC近似法計(jì)算的BMD/BMDL在每一種模型都有70%以上的數(shù)據(jù)集高于BIC估計(jì)法得到的BMD/BMDL。實(shí)例分析表明符合洪洞縣兒童體內(nèi)羥基代謝物劑量-反應(yīng)關(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估計(jì)中具有較好的穩(wěn)定性和魯棒性。

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

原理與方法

1.貝葉斯基準(zhǔn)劑量估計(jì)模型

基準(zhǔn)劑量估計(jì)首先要確定所研究的劑量-反應(yīng)數(shù)據(jù)類型和基準(zhǔn)反應(yīng),對(duì)于不同的數(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為單個(gè)模型的斜率和截距參數(shù):0≤a≤1,b≥0;g為背景參數(shù)。

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

BMRadd=f(BMD)-f(0)

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

2.貝葉斯模型平均法

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

BICK=-2LK+pKlog(N)

3.權(quán)重法

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

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

4.模型的驗(yàn)證與評(píng)價(jià)

研究分為模擬研究和實(shí)例研究?jī)蓚€(gè)部分,模擬研究分別計(jì)算出兩種方法的BMD以及BMDL值,設(shè)置BMR為0.01和0.001。本次研究借鑒Shao[5]等人的評(píng)價(jià)指標(biāo)觀察兩種方法中BMD/BMDL,來(lái)評(píng)價(jià)基準(zhǔn)劑量估計(jì)方法的可靠性,比值越低表明模型的參考區(qū)間越精確,說(shuō)明其方法的精確性和可靠性越高。同時(shí)觀察兩種方法中數(shù)據(jù)集預(yù)測(cè)失敗的例數(shù),失敗例數(shù)越低則魯棒性越高。

模擬研究

1.數(shù)據(jù)來(lái)源

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

2.兩種模型評(píng)價(jià)

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

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

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

實(shí)例分析

1.數(shù)據(jù)來(lái)源

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

2.分析結(jié)果

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

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

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

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

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

表5 不同模型的后驗(yàn)預(yù)測(cè)值和對(duì)應(yīng)權(quán)重

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

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

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

討 論

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

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

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