BIC與MCMC基準劑量估計模型的比較研究*

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

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

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

原理與方法

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

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

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為單個模型的斜率和截距參數:0≤a≤1,b≥0;g為背景參數。

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

BMRadd=f(BMD)-f(0)

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

2.貝葉斯模型平均法

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

BICK=-2LK+pKlog(N)

3.權重法

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

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

4.模型的驗證與評價

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

模擬研究

1.數據來源

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

2.兩種模型評價

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

表1 MCMC與BIC法對30個數據集計算BMD與BMDL的成功個數

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

實例分析

1.數據來源

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

2.分析結果

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

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

表4 洪洞縣兒童在不同劑量組下血常規異常人數

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

圖1 五種入選模型的劑量-反應關系

表5 不同模型的后驗預測值和對應權重

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

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

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

討 論

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

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

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