王青青 彭雯潔 崔 靖 張 婷 周小林,3△ 余紅梅,4,5△

【提 要】 目的 本文通過對貝葉斯基準(zhǔn)劑量估計(jì)軟件(bayesian benchmark dose,BBMD)進(jìn)行剖析,為進(jìn)一步理解和掌握該軟件奠定一定的理論基礎(chǔ)。方法 從特性、功能以及運(yùn)行3個(gè)方面對BBMD軟件進(jìn)行了介紹,以連續(xù)性劑量-反應(yīng)模型為例結(jié)合中國輻射防護(hù)研究院的數(shù)據(jù)實(shí)現(xiàn)了基準(zhǔn)劑量估計(jì)的分析。結(jié)果 BBMD以圖形和表格的形式輸出模型參數(shù)估計(jì)、后驗(yàn)概率p值、模型權(quán)重、劑量-反應(yīng)曲線、模型參數(shù)相關(guān)系數(shù)矩陣、基準(zhǔn)劑量反應(yīng)、模型的BMD估計(jì)值等結(jié)果,與BMDS或PROAST相比,BBMD具有基于最新研究的BMD分析附加功能、計(jì)算效率更高、劑量反應(yīng)建模更加可靠的優(yōu)點(diǎn)。結(jié)論 BBMD可以提供概率估計(jì),極大地滿足了當(dāng)前概率風(fēng)險(xiǎn)評估的需求。

過去幾十年里,關(guān)于人類健康風(fēng)險(xiǎn)評估中劑量-反應(yīng)(dose-response)關(guān)系的評定,通常采用無可見有害作用水平(no observed adverse effect level,NOAEL)。但傳統(tǒng)的NOAEL在評估劑量-反應(yīng)關(guān)系時(shí)存在一定缺陷。譬如,結(jié)果過分依賴于實(shí)驗(yàn)劑量的選擇;受樣本量大小的影響;NOAEL基于一個(gè)點(diǎn)的數(shù)據(jù),未考慮到劑量-反應(yīng)關(guān)系的斜率(曲線);當(dāng)一個(gè)實(shí)驗(yàn)未能發(fā)現(xiàn)NOAEL時(shí),安全閾值劑量只能基于最低可見有害作用水平(lowest observable adverse effect level,LOAEL),然后乘以10倍的不確定系數(shù)加以外推;由于不同實(shí)驗(yàn)提供的反應(yīng)水平上的作用劑量未必相同,所以NOAEL在可比性方面存在缺陷[1]?；贜OAEL的缺陷,1984年美國科學(xué)家Crump[2]提出了基準(zhǔn)劑量估計(jì)(benchmark dose,BMD),指引起預(yù)定頻率的不良健康劑量的統(tǒng)計(jì)學(xué)下限值?；鶞?zhǔn)劑量估計(jì)克服了NOAEL的不足,與NOAEL方法相比,BMD方法具有許多重要優(yōu)勢,但它需要更復(fù)雜的回歸算法將各種劑量反應(yīng)模型擬合到輸入數(shù)據(jù)中。因此,有必要利用完善的軟件來促進(jìn)BMD方法的實(shí)施。

目前已有兩大成熟的軟件用于基準(zhǔn)劑量的估計(jì),一個(gè)是由美國環(huán)境保護(hù)署于2000年研發(fā)的基準(zhǔn)劑量軟件BMDS(Benchmark Dose Software,2.6.0.1)。另一個(gè)軟件是PROAST(Programme for Analysing Dose-Response Data),由荷蘭國家公共衛(wèi)生與環(huán)境研究所(RIVM)推行,采用R軟件編程可應(yīng)用于各種系統(tǒng)(Windows、Mac OS、UNIX、Linux),并于2014年初更新至最新版本。BMDS和PROAST都有各自的優(yōu)勢,在一些技術(shù)細(xì)節(jié)上也大多不同。譬如,兩種軟件對劑量-反應(yīng)模型和連續(xù)數(shù)據(jù)分布的默認(rèn)假設(shè)不同,BMDS默認(rèn)變量服從正態(tài)分布、而PROAST則默認(rèn)變量服從對數(shù)正態(tài)分布;BMDS適用于Window系統(tǒng)下,具有良好的操作界面能夠分析多種類型的劑量反應(yīng)關(guān)系,包括最常用的計(jì)數(shù)資料和連續(xù)性資料。PROAST不僅可以分析計(jì)數(shù)資料和連續(xù)性資料,還可以分析等級資料?？偟膩碚f,這兩個(gè)軟件包都適用于劑量反應(yīng)分析和推導(dǎo)BMD及其統(tǒng)計(jì)下限(lower confidence limit of the benchmark dose,BMDL),然而這兩種軟件方法都基于概率論的統(tǒng)計(jì)方法(即極大似然估計(jì))進(jìn)行劑量反應(yīng)模型擬合和參數(shù)估計(jì)。

2009年美國國家研究委員會(United States National Research Council,NRC)實(shí)現(xiàn)了一項(xiàng)里程碑式的工作——發(fā)表了《科學(xué)與決策:推進(jìn)風(fēng)險(xiǎn)評估》[3],強(qiáng)調(diào)了在風(fēng)險(xiǎn)評估中以概率的形式量化特定風(fēng)險(xiǎn)劑量以支持監(jiān)管決策的重要性。2014年,世界衛(wèi)生組織(World Health Organization,WHO)發(fā)布了一份關(guān)于使用概率框架來協(xié)調(diào)風(fēng)險(xiǎn)表征方法的指導(dǎo)文件(IPCS 2014)[4]。2015年Fang[5]提出了貝葉斯基準(zhǔn)劑量估計(jì)(bayesian benchmark dose,BBMD),2018年Shao與Shapiro[6]一同開發(fā)了貝葉斯基準(zhǔn)劑量估計(jì)軟件,該軟件實(shí)現(xiàn)了基于貝葉斯推理的基準(zhǔn)劑量估計(jì)。

BBMD軟件是在BMDS軟件的基礎(chǔ)上新開發(fā)的一款網(wǎng)頁版本的貝葉斯基準(zhǔn)劑量分析系統(tǒng)。BBMD軟件建立在貝葉斯框架上,通過模型參數(shù)的先驗(yàn)分布來合并先驗(yàn)信息等特點(diǎn),使得BBMD系統(tǒng)與當(dāng)前流行的BMD估計(jì)軟件(如BMDS/PROAST)有著根本的不同。BBMD除了估計(jì)連續(xù)數(shù)據(jù)、計(jì)數(shù)資料的傳統(tǒng)BMD外,同時(shí)還能進(jìn)行批量的BMD分析。它屬于開源軟件,可以免費(fèi)使用,網(wǎng)址:https://benchmarkdose.org,或者h(yuǎn)ttps://benchmarkdose.com。

BBMD軟件介紹

1.特性方面

BBMD包括兩個(gè)模塊,負(fù)責(zé)計(jì)算和數(shù)據(jù)存儲的后端模塊和通過web瀏覽器與用戶交流的前端模塊。BBMD的內(nèi)部網(wǎng)站組件是由Python網(wǎng)絡(luò)框架(django)自動生成的,設(shè)計(jì)簡潔且實(shí)用。限制存儲在數(shù)據(jù)庫中的個(gè)人身份識別信息(personally identifiable information,PII)的數(shù)量,用戶的電子郵件是BBMD網(wǎng)站中唯一存儲的個(gè)人身份識別信息。

(1)基本原理:通過馬爾可夫鏈蒙特卡羅迭代擬合生成模型,并產(chǎn)生相應(yīng)的后驗(yàn)樣本,用戶給定基準(zhǔn)劑量反應(yīng)后可估計(jì)BMD。

(2)分布類型:BBMD軟件根據(jù)數(shù)據(jù)類型分別提供了8種模型。計(jì)數(shù)資料和連續(xù)資料一共16種模型分布。

計(jì)數(shù)資料的8種劑量-反應(yīng)模型分布:

①Q(mào)uantal-linear model:f(dose)=a+(1-a)×(1-e-b×dose);0≤a≤1,b≥0

③Probit model:f(dose)=Φ(a+b×dose);50≤a≤-50,b≥0

④Weibull model:f(dose)=a+(1-a)×(1-e-b×doseg);0≤a≤1,b≥0,g≥restriction

⑤Multistage-2ndmodel:f(dose)=a+(1-a)×(1-e-b×dose-c×dose2);0≤a≤1,b≥0,c≥0

0≤a≤1,b～Uniform(restriction,15),-5≤c≤15

⑦LogProbit model:f(dose)=a+(1-a)×Φ(c+b×log(dose));

0≤a≤1,b～Uniform(restriction,15),-5≤c≤15

0≤a≤1,b～Uniform(restriction,15),-5≤c≤15,0≤g≤1

連續(xù)變量的8種劑量-反應(yīng)模型:

①Linear model:f(dose)=a+b×dose;a≥0,b～Uniform(blower,bupper)②Power model:f(dose)=a+b×doseg;a≥0,b～Uniform(blower,bupper),g≥restriction

⑤Exponential 2 model:f(dose)=a×eb×dose;a≥0,bI～Unifrom(0,50);bD～Uniform(-50,0)

⑥Exponential 3 model:f(dose)=a+eb×doseg;

a≥0,bI～Unifrom(0,50);bD～Uniform(-50,0),g≥restriction

⑦Exponential 4 model:f(dose)=a×[c-(c-1)×e-b×dose];a≥0,b≥0,c≥0

⑧Exponential 5 model:f(dose)=a×[c-(c-1)×e-(b×dose)g];a≥0,b≥0,c≥0,g≥restriction

其中,參數(shù)“a”表示背景劑量水平下的效應(yīng);參數(shù)“b”表示效力,bI表示數(shù)據(jù)存在上升趨勢、bD數(shù)據(jù)存在下降趨勢;數(shù)據(jù)存在上升趨勢時(shí)blower=0、bupper=[(效應(yīng)指標(biāo)最大最小值之差)/(效應(yīng)指標(biāo)最大最小值所對應(yīng)的劑量水平之差)]×5;數(shù)據(jù)存在下降趨勢時(shí)blower=[(效應(yīng)指標(biāo)最大最小值之差)/(效應(yīng)指標(biāo)最大最小值所對應(yīng)的劑量水平之差)]×5,bupper=0;參數(shù)“c”表示漸進(jìn)項(xiàng);參數(shù)“d”表示劑量;參數(shù)“g”表示與劑量有關(guān)的功率參數(shù),其中restriction為用戶定義值可取0、0.5、1。

2.功能方面

BBMD軟件可以對連續(xù)型變量、計(jì)數(shù)資料進(jìn)行貝葉斯基準(zhǔn)劑量估計(jì)、概率參考劑量估計(jì)以及處理批量的BMD分析。BBMD軟件輸出結(jié)果以文本和圖形兩種形式呈現(xiàn),輸出內(nèi)容包含相關(guān)的統(tǒng)計(jì)學(xué)描述、原始劑量響應(yīng)數(shù)據(jù)和90%置信區(qū)間的擬合曲線圖等。

3.運(yùn)行/操作方面

(1)BBMD軟件基本分析步驟

①輸入一個(gè)劑量-反應(yīng)數(shù)據(jù)集、設(shè)定馬爾可夫鏈蒙特卡羅算法和選擇劑量-反應(yīng)模型;

②完成第一步的操作后系統(tǒng)會進(jìn)行模型擬合并生成模型參數(shù)的估計(jì)值、擬合分布圖、并進(jìn)行交叉模型的比較;

③BMR設(shè)定、先驗(yàn)?zāi)Ｐ蜋?quán)重的設(shè)定;

④輸出結(jié)果。操作步驟如圖1。

圖1 BBMD軟件基本分析步驟流程圖

(2)BBMD軟件選項(xiàng)卡

注冊并登錄個(gè)人賬戶后,點(diǎn)擊右上角的“New analysis”選項(xiàng),進(jìn)入分析編輯界面,會話界面有9個(gè)選項(xiàng)卡。

①Dataset:Dataset子窗戶包括:dataset type和dataset兩個(gè)小部分。Dataset type包含連續(xù)型變量、計(jì)數(shù)資料兩大類,每一大類又分別細(xì)分為個(gè)體數(shù)據(jù)、匯總數(shù)據(jù);dataset(數(shù)據(jù)錄入模塊),包括劑量(dose)、研究對象(n)、發(fā)病率(incidence),單擊save dataset,可以保存當(dāng)前數(shù)據(jù)集。數(shù)據(jù)集保存成功后,數(shù)據(jù)將以圖形和表格的形式顯示在子對話框的右邊。

②MCMC setting:包含四個(gè)子模塊,迭代次數(shù)(markov chain iterations/per chain)、馬爾可夫鏈數(shù)(number of markov chains)、預(yù)熱百分比(warmup percent,%)、種子數(shù)(random seed)。

Markov chain iterations/per chain指MCMC鏈的長度,即每個(gè)MCMC鏈中的后驗(yàn)樣本的數(shù)量,默認(rèn)值為30,000,允許范圍:10000到50000之間的任意整數(shù)。Number of markov chains,要采樣的馬爾可夫鏈數(shù),默認(rèn)值為1,允許的范圍為1～3條。warmup percent,將放棄每個(gè)馬爾可夫鏈中的樣本百分比,默認(rèn)值為50%,允許的范圍為10%～90%。random seed,MCMC算法中使用的隨機(jī)種子數(shù),隨機(jī)生成允許范圍:0～99,999。

根據(jù)迭代次數(shù)、鏈數(shù)、預(yù)熱百分比可以算出最終的后驗(yàn)樣本數(shù):

30000×1×(1-50%)=15000

最后單擊save run setting,保存MCMC設(shè)置。

③Model settings:包含create a model和add standard models兩大模塊,其中create a model包含model type和model name兩個(gè)小模塊。model type中有8個(gè)模型類別可供選擇,對于具有power參數(shù)的模型,需要對power選擇一個(gè)限制值。當(dāng)前系統(tǒng)中有5個(gè)選項(xiàng):0、0.25、0.5、0.75和1。

④Execute model fit:模型運(yùn)行的3個(gè)先決條件:數(shù)據(jù)創(chuàng)建完畢、MCMC設(shè)備已保存、至少選擇了一個(gè)模型。

⑤Model fit result:模型運(yùn)行后,輸出模型的相應(yīng)結(jié)果。

⑥Risk at dose:根據(jù)第一步dataset中設(shè)定的dose的劑量組,分別估計(jì)每個(gè)劑量組的風(fēng)險(xiǎn)。

⑦BMD estimates:設(shè)定BMR的值(通常為5%、10%),分別估計(jì)BMD、BMDL。

⑧RfD estimates:運(yùn)用BMD estimates估計(jì)得到BMDL后,根據(jù)BMDL估計(jì)參考劑量(reference dose,Rfd)值。

⑨Report:通過Email以word、excel格式將結(jié)果傳送至用戶。

應(yīng)用實(shí)例

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

目前國內(nèi)關(guān)于運(yùn)用BBMD軟件估計(jì)BMD的相關(guān)報(bào)道較少,本文利用中國輻射防護(hù)研究院山西省洪洞縣數(shù)據(jù),研究2-羥基芴與紅細(xì)胞分布寬度(red cell volume distribution width,RDW)之間的劑量-反應(yīng)關(guān)系。

在該數(shù)據(jù)中,選取部分成年居民的暴露信息,將2羥基芴作為暴露變量(dose),紅細(xì)胞分布寬度標(biāo)準(zhǔn)差(RDW-SD)作為反應(yīng)變量(response)。選用BBMD軟件估計(jì)2-羥基芴與紅細(xì)胞分布寬度之間的劑量-反應(yīng)關(guān)系。詳見表1。

表1 暴露資料分布表

(1)數(shù)據(jù)輸入

單擊dataset選項(xiàng)卡,選擇continuous(individual)數(shù)據(jù)類型,依次輸入dose和response的數(shù)值。

(2)MCMC設(shè)定

默認(rèn)系統(tǒng)的設(shè)定:迭代次數(shù)30000次,預(yù)熱百分比50%,馬爾可夫鏈1條,隨機(jī)種子數(shù)21624。

(3)模型的選擇

單擊add standard models后,系統(tǒng)默認(rèn)選擇連續(xù)變量的8個(gè)標(biāo)準(zhǔn)模型。

(4)基準(zhǔn)反應(yīng)的設(shè)定

在 BMD estimation method里選擇混合方法估計(jì)(Hybrid method(tails)),adversity measure框里選擇對照組百分比(control group percentile),adversity value框里輸入0.56,Benchmark response value設(shè)置為0.1。最終BBMD軟件根據(jù)模型權(quán)重對8個(gè)模型進(jìn)行權(quán)重均衡后計(jì)算2-羥基芴的健康指導(dǎo)值[7]。

2.結(jié)果呈現(xiàn)

BBMD系統(tǒng)以圖形和表格的形式輸出相應(yīng)的結(jié)果。本文演示了連續(xù)變量的劑量-反應(yīng)關(guān)系,軟件輸出如下結(jié)果:模型參數(shù)估計(jì)、后驗(yàn)概率p值、模型權(quán)重和劑量-反應(yīng)曲線(圖2)、模型參數(shù)的分布、后驗(yàn)樣本描繪、模型參數(shù)相關(guān)系數(shù)矩陣、基準(zhǔn)劑量反應(yīng)(圖3);BMD估計(jì)的概率分布圖(圖5)和模型的BMD估計(jì)值、額外風(fēng)險(xiǎn)和附加風(fēng)險(xiǎn)值(圖6)。

圖2 模型擬合結(jié)果

圖3 模型參數(shù)分布圖

圖4 模型估計(jì)結(jié)果

圖5 模型概率分布函數(shù)圖

圖6 風(fēng)險(xiǎn)概率分布函數(shù)圖

2-羥基芴暴露發(fā)生紅細(xì)胞分布寬度偏低的8個(gè)模型擬合的BMDL分別為0.31 μmol/mol、4.40 μmol/mol、0.18 μmol/mol、4.32 μmol/mol、2.34 μmol/mol、4.72 μmol/mol、0.25 μmol/mol、0.28 μmol/mol;BMD分別為0.45 μmol/mol、7.01 μmol/mol、0.38 μmol/mol、6.85 μmol/mol、6.51 μmol/mol、7.08 μmol/mol、0.40 μmol/mol、0.42 μmol/mol。通過對8個(gè)模型進(jìn)行權(quán)重均衡,最終得到的BMDL的平均值為1.80 μmol/mol,BMD的平均值為6.91 μmol/mol。詳見圖4。

討 論

本文通過特性、功能以及運(yùn)行3個(gè)方面對BBMD軟件進(jìn)行了介紹,利用中國輻射防護(hù)研究院收集的數(shù)據(jù),利用BBMD軟件實(shí)現(xiàn)了基準(zhǔn)劑量估計(jì)的分析,使讀者對BBMD有一個(gè)初步的認(rèn)識,對理解和掌握該軟件的基本操作及其性能有重要的意義。

迄今為止,已有不少軟件可以進(jìn)行基準(zhǔn)劑量估計(jì)。譬如,世界公認(rèn)的兩大基準(zhǔn)劑量估計(jì)軟件BMDS和PROAST,還有@Risk[8]等其他專業(yè)軟件。BBMD屬于一種新興的基準(zhǔn)劑量估計(jì)軟件,它提供了一種通過先驗(yàn)分布整合先驗(yàn)信息的方法建模并進(jìn)行后續(xù)的基準(zhǔn)劑量估計(jì),這對于提高數(shù)據(jù)的質(zhì)量,建模的可靠性和魯棒性具有一定幫助,同時(shí)也可能是風(fēng)險(xiǎn)評估人員唯一可用的方法。與其他基準(zhǔn)劑量估計(jì)軟件(如BMDS/PROAST)相比,BBMD軟件的亮點(diǎn)在于:第一,具有基于最新研究的BMD分析附加功能。第二,將模型擬合和BMD分析的步驟解耦,使得系統(tǒng)計(jì)算效率更高[5-6]。第三,允許更大的分析靈活性。最重要的是,BBMD采用貝葉斯估計(jì),結(jié)合先驗(yàn)信息,使劑量反應(yīng)建模更加可靠,并且可以為劑量反應(yīng)評估中重要的參數(shù)提供分布估計(jì),這極大地促進(jìn)了概率風(fēng)險(xiǎn)評估的發(fā)展[9-10]。因此,關(guān)于劑量-反應(yīng)關(guān)系評定方面,BBMD是一種有效的替代工具,并且BBMD軟件代表了劑量-反應(yīng)評估領(lǐng)域的最新科學(xué)方法和技術(shù)[5,7]。

不同于其他基準(zhǔn)劑量估計(jì)軟件采用極大似然法估計(jì)模型參數(shù),BBMD基于貝葉斯理論框架,通過模型參數(shù)的先驗(yàn)分布合并先驗(yàn)信息進(jìn)行模型參數(shù)估計(jì)[6,11-12]。BBMD不局限于單一模型結(jié)果,而是以各模型的后驗(yàn)概率為權(quán)重,對單一模型的預(yù)測值進(jìn)行加權(quán)平均得到一個(gè)綜合的模型估計(jì)值。BBMD估計(jì)結(jié)果合理之處在于,將不同模型的不確定性考慮在內(nèi),待估計(jì)量是基于備選模型后驗(yàn)概率的加權(quán)平均,能給予擬合良好的模型更大的權(quán)重,擬合不好的模型較小的權(quán)重。此外,BBMD可以通過設(shè)置不同的先驗(yàn)概率分布,將模型選擇的外部信息結(jié)合到數(shù)據(jù)分析中,既考慮了模型的先驗(yàn)信息,又提供數(shù)據(jù)的樣本信息,分析更具科學(xué)性,在一定程度上可消除劑量-反應(yīng)關(guān)系模型的不確定性。再者,包括美國國家研究委員會、WHO等機(jī)構(gòu)均提倡運(yùn)用貝葉斯框架協(xié)助概率風(fēng)險(xiǎn)評估[3-4]。除此之外,BBMD軟件還在不斷更新與升級,它的開發(fā)是為了支持人類健康風(fēng)險(xiǎn)評估中劑量-反應(yīng)的估計(jì),促進(jìn)BMD方法的實(shí)現(xiàn)。相信在未來,會有越來越多的科研工作者將貝葉斯基準(zhǔn)劑量估計(jì)應(yīng)用于危險(xiǎn)度評價(jià),BBMD軟件也將成為基準(zhǔn)劑量估計(jì)的主流軟件。