貝葉斯基準劑量統計軟件*

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

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

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

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

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

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

BBMD軟件介紹

1.特性方面

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

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

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

計數資料的8種劑量-反應模型分布:

①Quantal-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

連續變量的8種劑量-反應模型:

①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

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

2.功能方面

BBMD軟件可以對連續型變量、計數資料進行貝葉斯基準劑量估計、概率參考劑量估計以及處理批量的BMD分析。BBMD軟件輸出結果以文本和圖形兩種形式呈現,輸出內容包含相關的統計學描述、原始劑量響應數據和90%置信區間的擬合曲線圖等。

3.運行/操作方面

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

①輸入一個劑量-反應數據集、設定馬爾可夫鏈蒙特卡羅算法和選擇劑量-反應模型;

②完成第一步的操作后系統會進行模型擬合并生成模型參數的估計值、擬合分布圖、并進行交叉模型的比較;

③BMR設定、先驗模型權重的設定;

④輸出結果。操作步驟如圖1。

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

(2)BBMD軟件選項卡

注冊并登錄個人賬戶后,點擊右上角的“New analysis”選項,進入分析編輯界面,會話界面有9個選項卡。

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

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

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

根據迭代次數、鏈數、預熱百分比可以算出最終的后驗樣本數:

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

最后單擊save run setting,保存MCMC設置。

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

④Execute model fit:模型運行的3個先決條件:數據創建完畢、MCMC設備已保存、至少選擇了一個模型。

⑤Model fit result:模型運行后,輸出模型的相應結果。

⑥Risk at dose:根據第一步dataset中設定的dose的劑量組,分別估計每個劑量組的風險。

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

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

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

應用實例

1.數據來源

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

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

表1 暴露資料分布表

(1)數據輸入

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

(2)MCMC設定

默認系統的設定:迭代次數30000次,預熱百分比50%,馬爾可夫鏈1條,隨機種子數21624。

(3)模型的選擇

單擊add standard models后,系統默認選擇連續變量的8個標準模型。

(4)基準反應的設定

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

2.結果呈現

BBMD系統以圖形和表格的形式輸出相應的結果。本文演示了連續變量的劑量-反應關系,軟件輸出如下結果:模型參數估計、后驗概率p值、模型權重和劑量-反應曲線(圖2)、模型參數的分布、后驗樣本描繪、模型參數相關系數矩陣、基準劑量反應(圖3);BMD估計的概率分布圖(圖5)和模型的BMD估計值、額外風險和附加風險值(圖6)。

圖2 模型擬合結果

圖3 模型參數分布圖

圖4 模型估計結果

圖5 模型概率分布函數圖

圖6 風險概率分布函數圖

2-羥基芴暴露發生紅細胞分布寬度偏低的8個模型擬合的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個模型進行權重均衡,最終得到的BMDL的平均值為1.80 μmol/mol,BMD的平均值為6.91 μmol/mol。詳見圖4。

討 論

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

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

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