基準(zhǔn)劑量估計(jì)的非參數(shù)貝葉斯方法研究與應(yīng)用*

顧彩姣 王 濤 王 彤

基準(zhǔn)劑量估計(jì)的非參數(shù)貝葉斯方法研究與應(yīng)用*

顧彩姣1,2王 濤2王 彤1△

目的比較兩種基準(zhǔn)劑量估計(jì)的非參數(shù)貝葉斯方法在不同劑量反應(yīng)數(shù)據(jù)情形下的表現(xiàn)，再與參數(shù)模型作比較。方法介紹基于加權(quán)過(guò)程和隨機(jī)過(guò)程的非參數(shù)貝葉斯方法的基本原理，通過(guò)模擬分析與實(shí)例研究比較各方法結(jié)果的差異。結(jié)果模擬研究顯示，兩種非參數(shù)方法NPB1和NPB2的估計(jì)結(jié)果都與真實(shí)BMD值較為接近，NPB2方法的BMDL覆蓋率更接近真實(shí)水平，因此更為可??；實(shí)例中與參數(shù)模型結(jié)果相比，非參數(shù)方法得到的BMD估計(jì)值都落入或非常接近由參數(shù)模型得到的BMD估計(jì)值的范圍，而非參數(shù)方法估計(jì)的BMDLs值大多比參數(shù)模型算得的BMDLs值要小。結(jié)論兩種非參數(shù)貝葉斯方法都可提供合理的擬合值，尤其是在傳統(tǒng)的參數(shù)模型無(wú)法得到擬合結(jié)果的情況下，也可很好的應(yīng)用；NPB2方法在估計(jì)結(jié)果和軟件運(yùn)算速度方面均略優(yōu)于NPB1。

基準(zhǔn)劑量 BMDS軟件 非參數(shù)方法 狄利克雷分布 布朗運(yùn)動(dòng)

過(guò)去幾十年一直采用未觀察到有害作用劑量(no-observed-adverse-effect-level,NOAEL)來(lái)評(píng)估非致癌物的健康風(fēng)險(xiǎn)[1]。但這種方法的使用依賴于大量定義好的限制條件，比如劑量選擇、劑量間隔、在選定劑量水平時(shí)使用的動(dòng)物數(shù)量等，且該方法也未考慮劑量反應(yīng)曲線的形狀等信息。因此，科學(xué)家們開(kāi)始探索其替代方法?；鶞?zhǔn)劑量法(Benchmark Dose,BMD)是由美國(guó)科學(xué)家Crump于1984年最先提出的[2]，該方法提出后引起了科學(xué)家們的廣泛關(guān)注，并將大量歷史資料進(jìn)行了運(yùn)算驗(yàn)證，他們認(rèn)為BMD彌補(bǔ)了NOAEL的不足之處?；鶞?zhǔn)劑量方法也被用在發(fā)育毒物以及神經(jīng)毒物終點(diǎn)事件研究中，很多文獻(xiàn)認(rèn)為基準(zhǔn)劑量法在毒理學(xué)中的應(yīng)用已經(jīng)多于其在風(fēng)險(xiǎn)評(píng)估中的應(yīng)用了[3-5]。

目前文獻(xiàn)研究中基準(zhǔn)劑量的估計(jì)方法主要有：針對(duì)不同類型反應(yīng)數(shù)據(jù)的參數(shù)模型；模型平均法；半?yún)?shù)方法和非參數(shù)方法等。參數(shù)模型在使用中需進(jìn)行模型選擇，因此存在不確定性[6-7]；模型選擇法是在一定模型空間內(nèi)將不確定性移除，但未克服不確定性[7]；形式較為靈活的半?yún)?shù)方法和非參數(shù)方法近年來(lái)使用較多，半?yún)?shù)方法計(jì)算復(fù)雜[8-9]，因此本文主要研究非參數(shù)方法[10]。

原理與方法

1.基準(zhǔn)劑量

2.基于加權(quán)過(guò)程的非參數(shù)貝葉斯方法[12](NPB1)

該方法的主要思想是將反應(yīng)概率函數(shù)p(d)定義為關(guān)于d的光滑、單調(diào)非減、取值范圍在0到1之間的函數(shù)，即：

(1)

πΔ(Δ)=Dirichlet(α1,…,αm+1)

(2)

本文中超參數(shù)α1,…,αm+1指定為(1,1,…,1)/(m+1)。

在貝葉斯方法里，寬度參數(shù)h是未知的，可指定一個(gè)先驗(yàn)信息。在將劑量水平標(biāo)化為單位區(qū)間后，假定寬度參數(shù)小于1，并使用Beta分布：

πh(h)=Beta(α,β)。

(3)

將反應(yīng)率記為p(d;Δ,h)，表明概率值對(duì)參數(shù)Δ和h的依賴性。此時(shí)參數(shù)估計(jì)的似然函數(shù)為：

(4)

π(Δ,h|y1,…,ym)∝L(Δ,h|y1,…,

ym)πΔ(Δ)πh(h)

(5)

其中ni為每個(gè)劑量下的實(shí)驗(yàn)動(dòng)物數(shù)，yi為每個(gè)劑量組內(nèi)結(jié)局為陽(yáng)性反應(yīng)的動(dòng)物數(shù)。為采用MARKOV CHAIN MONTE CARLO(MCMC)方法從后驗(yàn)信息中抽樣，標(biāo)準(zhǔn)的Metropolis-Hastings 算法可從近似的后驗(yàn)中得到大量的L樣本。令(Δ(1),h(1)),…,(Δ(L),h(L))為后驗(yàn)信息中抽出的L個(gè)樣本。由這些后驗(yàn)樣本可得到全部的后驗(yàn)曲線p(d;Δ(l),h(l))，l=1,…,L，劑量反應(yīng)關(guān)系的后驗(yàn)估計(jì)為：

(6)

貝葉斯方法的有用之處在于后驗(yàn)樣本可提供所有重要未知量的推斷，未知量

BMD(l)=BMD(p(d,Δ(l)),h(l),γ)，l=1,…,L

(7)

BMD可由后驗(yàn)估計(jì)均值得到：

(8)

3.基于隨機(jī)過(guò)程的非參數(shù)貝葉斯方法[12](NPB2)

該方法假定劑量反應(yīng)關(guān)系是由基于d的隨機(jī)過(guò)程得到的，該隨機(jī)過(guò)程需確保每個(gè)樣本都滿足光滑、單調(diào)遞增，p(d)取值范圍在0到1之間的一般特征。本文使用的方法是由一種特殊的高斯過(guò)程-布朗運(yùn)動(dòng)(brownian motion, BM)得到的。

假定ω(d)是關(guān)于d的標(biāo)準(zhǔn)布朗運(yùn)動(dòng)，則有ω(d1)=0，且當(dāng)di

(9)

整合的BM形式已經(jīng)被用在二分回歸問(wèn)題的非參數(shù)貝葉斯建模過(guò)程中，近年來(lái)研究者們也都探討過(guò)整合指數(shù)BM方法的一些性質(zhì)。整合隨機(jī)過(guò)程的指數(shù)形式為遞增的非負(fù)函數(shù)，也已被用在了光滑單調(diào)回歸方法中。

劑量反應(yīng)關(guān)系的模型為：

p(d)=p0+(1-p0)z(d)

(10)

對(duì)照組的反應(yīng)率應(yīng)滿足0≤p0=p(0)≤1，為了保持范圍限制，如0≤p(d)≤1，我們將BM的形式截取為：

(11)

為實(shí)現(xiàn)該方法，我們使用離散的BM形式，給定增量ω(dj)-ω(dj-1)，j=1,…,m。則整合的指數(shù)BM的形式近似為：

(12)

對(duì)于每個(gè)j，ω(xj)-ω(xj-1)～N(0,λ(xj-xj-1))，λ為方差參數(shù)。同時(shí)，假定基線反應(yīng)p0的先驗(yàn)分布為Beta(a,b)，此時(shí)，似然和先驗(yàn)相互獨(dú)立，為：

(13)

之后BMD和BMDL值的估計(jì)方法與NPB1類似。

模擬研究

1.模擬數(shù)據(jù)設(shè)置

假定多階段模型為真實(shí)的劑量反應(yīng)關(guān)系，模型包含四個(gè)參數(shù)q=(q0,q1,q2,q3)，劑量反應(yīng)函數(shù)為p(d)=1-exp(-q0-q1d-q2d2-q3d3),該多階段模型用MS(q0,q1,q2,q3)來(lái)表示。假定每個(gè)劑量點(diǎn)di的動(dòng)物數(shù)均為ni=50，劑量反應(yīng)數(shù)據(jù)服從二項(xiàng)分布B(ni,p(di))。在使用貝葉斯方法時(shí)，設(shè)定MCMC抽樣的N=10000，B=5000；NPB1方法中每個(gè)樣本的寬度參數(shù)h均使用beta(0.7,1)，而狄利克雷參數(shù)α1,α2,…,αm+1為(1,1,…,1)/(m+1)；NPB2中，假定p0的先驗(yàn)為beta(0.1,1)，方差參數(shù)λ取0.0001。在函數(shù)p的表達(dá)式已知，額外風(fēng)險(xiǎn)設(shè)置為10%的情況下，我們可以計(jì)算得到每種數(shù)據(jù)下的真實(shí)的BMD值。

參照其他文獻(xiàn)研究[13]，本模擬設(shè)定兩種不同的參數(shù)配置，分別為：q=(0,1,1,3)和q=(0.05,0,3,0)；劑量組數(shù)m=4和m=6兩種情況；劑量設(shè)置為等間距和近似對(duì)數(shù)間距兩種情況。以上參數(shù)設(shè)置構(gòu)成8種不同情形的劑量反應(yīng)數(shù)據(jù)，在每種情況下，都產(chǎn)生R=500的獨(dú)立蒙特卡洛樣本，每個(gè)樣本都用NPB1，NPB2兩種方法估計(jì)，以BMD估計(jì)值與真實(shí)值的差距以及BMDL的覆蓋率為評(píng)價(jià)指標(biāo)，比較這兩種方法在不同劑量反應(yīng)數(shù)據(jù)情形下的表現(xiàn)。BMDL的覆蓋率是由R個(gè)估計(jì)值中，包含了真實(shí)BMD值的100(1-α)%可信區(qū)間所占的比例得到。我們分別取α=0.10和0.05，得到對(duì)應(yīng)的90%和95%可信區(qū)間覆蓋率。

表1 模擬研究結(jié)果

2.模擬結(jié)果

由表1中的結(jié)果可以看到， NPB1和NPB2兩種方法得到的BMD后驗(yàn)估計(jì)均值都與真實(shí)BMD值較為接近。當(dāng)多階段模型參數(shù)設(shè)置為MS(0,1,1,3)、劑量組數(shù)為6，間距為對(duì)數(shù)時(shí)，NPB1方法的覆蓋率明顯低于正常水平，表明此時(shí)的BMD估計(jì)值有較大偏倚。而NPB2方法的覆蓋率更接近真實(shí)水平，且運(yùn)算速度較快，模擬一次的時(shí)間是NPB1方法的1/10左右，就本次研究來(lái)講，NPB2方法較NPB1方法更為可取。

實(shí)例分析

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

從Nitcheva等人2007年分析過(guò)的數(shù)據(jù)集中，選取9組癌癥數(shù)據(jù)進(jìn)行分析[14]。

表2 選自Nitcheva et al.(2007)的9組癌癥數(shù)據(jù)實(shí)例

2.實(shí)例分析結(jié)果

對(duì)于每個(gè)數(shù)據(jù)集，都采用BMDS軟件中包含的9種二分反應(yīng)模型[15]，以及NPB1和NPB2分別來(lái)估計(jì)劑量反應(yīng)關(guān)系，BMD和BMDL值見(jiàn)表3。

表3 參數(shù)方法與非參數(shù)方法估計(jì)結(jié)果的比較

表3中是兩種非參數(shù)方法和BMDS軟件模型包中參數(shù)模型的估計(jì)結(jié)果。BMDS軟件的結(jié)果包括兩部分：由AIC值最小原則選出的最優(yōu)參數(shù)模型的估計(jì)結(jié)果；滿足模型擬合優(yōu)度檢驗(yàn)界值的所有參數(shù)模型估計(jì)出的BMDs和BMDLs的范圍。9組實(shí)例由兩種非參數(shù)方法得到的估計(jì)結(jié)果比較接近，而NPB2方法的估計(jì)值略低于NPB1方法得到的估計(jì)值。而兩種非參數(shù)方法的BMD估計(jì)值大多都在實(shí)驗(yàn)劑量范圍之內(nèi)，而已有證據(jù)都表明在這個(gè)范圍內(nèi)可以觀察到10%的額外風(fēng)險(xiǎn)發(fā)生。

與不同的參數(shù)模型結(jié)果相比，兩種非參數(shù)方法得到的BMD估計(jì)值都落入或非常接近由參數(shù)模型得到的BMD估計(jì)值的范圍；而在可信區(qū)間方面，非參數(shù)方法估計(jì)的BMDLs值大多比參數(shù)模型算得的BMDLs值要小，這很可能是因?yàn)榉菂?shù)方法考慮的功能空間更大。對(duì)于某些劑量反應(yīng)數(shù)據(jù)類型，在常用的參數(shù)模型無(wú)法給出合理擬合值的情況下，非參數(shù)方法仍可較好的使用，這也充分說(shuō)明了非參數(shù)方法在其使用中的靈活性和應(yīng)用范圍的廣泛性。

討 論

在擬合毒理實(shí)驗(yàn)動(dòng)物數(shù)據(jù)的劑量反應(yīng)曲線時(shí)，使用參數(shù)模型通常要根據(jù)AIC最小原則進(jìn)行模型選擇。Webster等人2012年的一篇文章對(duì)BMD估計(jì)中模型不確定性及其影響做了詳細(xì)研究。模擬研究表明[16]，當(dāng)每組劑量實(shí)驗(yàn)動(dòng)物數(shù)比較少(低于50)時(shí)，通過(guò)AIC原則選出正確模型的概率均明顯很低，且模型包含的參數(shù)越多，每組動(dòng)物數(shù)越少，被選為正確模型的概率越低；當(dāng)樣本量增加時(shí)，該概率值會(huì)升高，但樣本大小為1000時(shí)，選出正確模型的概率最高也只能達(dá)到77%。另外，通過(guò)最優(yōu)模型計(jì)算得到的BMDL值與指定的真實(shí)模型的BMD值的比較發(fā)現(xiàn)，BMDL值低于BMD值的概率值遠(yuǎn)低于理想水平95%。一旦模型指定錯(cuò)誤，它估計(jì)的BMDL值便會(huì)大于選定的BMR值，通常會(huì)超出6倍以上。在實(shí)際的風(fēng)險(xiǎn)評(píng)估過(guò)程中，這種超出率被認(rèn)為是無(wú)法接受的。

非參數(shù)方法在進(jìn)行劑量反應(yīng)曲線擬合過(guò)程中的使用非常靈活，在某些數(shù)據(jù)情形下參數(shù)模型無(wú)法使用，非參數(shù)方法依然可以得到很好的估計(jì)值，并且解決了模型不確定性的問(wèn)題，作為目前可用的參數(shù)模型的補(bǔ)充與替代方法，有很重要的意義和使用價(jià)值。常用的非參數(shù)方法有多種，如分位數(shù)回歸，保序回歸，非參數(shù)貝葉斯方法等。分位數(shù)回歸在實(shí)驗(yàn)劑量點(diǎn)較少時(shí)，通常估計(jì)偏倚較大；保序回歸BMDL的覆蓋率也比貝葉斯方法的低；而非參數(shù)貝葉斯方法可通過(guò)指定先驗(yàn)來(lái)合并一些背景信息，尤其是在劑量點(diǎn)較少時(shí)，該特點(diǎn)讓非參數(shù)貝葉斯方法應(yīng)用更廣。

在實(shí)際應(yīng)用中，當(dāng)樣本量很大(N=1000)時(shí)相應(yīng)的模型選擇的正確率較高，可根據(jù)參數(shù)模型估計(jì)相應(yīng)的BMD和BMDL值，但BMDS軟件中包含的參數(shù)模型較多，需要用各個(gè)模型分別擬合數(shù)據(jù)，后根據(jù)模型擬合優(yōu)度檢驗(yàn)結(jié)果及AIC原則等模型選擇方法選出最優(yōu)模型，整個(gè)計(jì)算過(guò)程比較繁瑣，還要承擔(dān)模型錯(cuò)誤指定的風(fēng)險(xiǎn)，因此，推薦使用本課題提到的非參數(shù)貝葉斯方法。雖然兩非參數(shù)貝葉斯方法原理部分較難理解，但其相應(yīng)的R程序均已編寫完成，使得應(yīng)用非常方便、省時(shí)，結(jié)果卻更加保守可靠。而NPB2方法與NPB1方法相比，在擬合結(jié)果精度和軟件運(yùn)算速度方面優(yōu)勢(shì)更明顯，更推薦NPB2方法。

[1] Davis JA,Gift JS,Zhao QJ.Introduction to benchmark dose methods and U.S.EPA′s benchmark dose software(BMDS)version 2.1.1..Toxicology and Applied Pharmacology,2011,254:181-191.

[2] Crump KS.A new method for determining allowable daily intakes.Fundamental and Applied Toxicology,1984,4:854-871.

[3] European Chemicals Bureau(ECB).Technical guidance document on risk assessment of chemical substances following european regulations and directives,parts i-iv.http://europa.eu.int,2016-03-20.

[4] U.S.General Accounting Office.Selected federal agencies′ procedures,assumptions,and policies.Washington,2001-08:64-151.

[5] OECD.Current approaches in the statistical analysis of ecotoxicity data:A guidance to application.Environment Directorate,Organisation For Economic Co-Operation and Development.www.univet.hu,2006.

[6] Faustman EM,Bartell SM.Review of noncancer risk assessment:Applications of benchmark dose methods.Hum.Ecol.Risk Assess,1997,3(5):893-920.

[7] Kang SH,Kodell RL,Chen JJ.Incorporating model uncertainties along with data uncertainties in microbial risk assessment.Regul.Toxicol.Pharmacol,2000,32(1):68-72.

[8] Wheeler MW,Bailer AJ.Monotonic Bayesian Semiparametric Benchmark.Risk Analysis,2012,32(7):1207-1218.

[9] Bosch RJ,Wypij D,Ryan LM.A semiparametric approach to risk assessment for quantitative outcomes.Risk Analysis,1996,16(5):657-665.

[10]Bhattacharya R,Lin L.Nonparametric benchmark analysis in risk assessment_ A Comparative study by simulation and data analysis.Statistics and Probability Letters,2010,80:1947-1953.

[11]張俊杰,張海燕.基準(zhǔn)劑量評(píng)估系統(tǒng)的研究與實(shí)現(xiàn).北京林業(yè)大學(xué),2013:1-45.

[12]Guha G,Roy A,Kopylev L,et al.White P.Nonparametric Bayesian Methods for Benchmark Dose Estimation.Risk Analysis,2013,33(9):1608-1619.

[13]Kopylev L,Fox J.Parameters of a Dose-Response Model Areon the Boundary:What Happens with BMDL.Risk Analysis,2009,29(1):18-25.

[14]Nitcheva DK,Piegorsch WW,West RW.On use of the multistage dose-response model for assessing laboratory animal.Regulatory.Toxicology and Pharmacology,2007,48:135-147.

[15]顧劉金,陳瓊姜,吳立仁,等.基準(zhǔn)劑量統(tǒng)計(jì)軟件BMDS簡(jiǎn)介.中國(guó)衛(wèi)生統(tǒng)計(jì)，2005,22(4):257-258.

TheStudyandApplicationofNonparametricBayesianMethodsforBenchmarkDoseEstimation

Gu Caijiao,Wang Tao,Wang Tong

(DepartmentofHealthStatistics,ShanxiMedicalUniversity(030001),Shanxi)

ObjectiveComparing the performance of the two nonparametric Bayesian methods for benchmark dose estimating under different dose response data,then comparing them with traditional parametric methods.MethodsIntroduce the basic principle of the nonparametric Bayesian method based on weighted process and stochastic process separately,then compared the estimations through simulate study and instance analysis.ResultsThe simulate study shows that the posterior estimates were reasonably close to the target true BMD value for the two nonparametric methods,and NBP2 is more desirable compared to NPB1.The nine examples indicate that the BMD estimates from the nonparametric approaches generally fall into or very near the interval of those obtained from BMDS and nonparametric approaches tend to produce lower BMDLs than the parametric modeling approaches.ConclusionThe posterior estimates were reasonably close to the target true BMD value for the two nonparametric methods,especially when standard parametric models fail to fit to the data adequately.The NPB2 method is slightly bet-ter than the NPB1 method in the aspect of estimation result and the software operation speed.

Benchmark dose;BMDS software;Nonparametric method;Dirichlet distribution;Brownian motion

國(guó)家自然科學(xué)基金(81473073)

1.山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室(030001) 2.海淀區(qū)東升鎮(zhèn)社區(qū)衛(wèi)生服務(wù)中心預(yù)防保健科

△通信作者：王彤，E-mail:wtstat1@sina.com

郭海強(qiáng))