顧彩姣 王 濤 王 彤

基準劑量估計的非參數(shù)貝葉斯方法研究與應用*

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

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

基準劑量 BMDS軟件 非參數(shù)方法 狄利克雷分布 布朗運動

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

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

原理與方法

1.基準劑量

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

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

(1)

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

(2)

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

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

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

(3)

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

(4)

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

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

(5)

其中ni為每個劑量下的實驗動物數(shù)，yi為每個劑量組內結局為陽性反應的動物數(shù)。為采用MARKOV CHAIN MONTE CARLO(MCMC)方法從后驗信息中抽樣，標準的Metropolis-Hastings 算法可從近似的后驗中得到大量的L樣本。令(Δ(1),h(1)),…,(Δ(L),h(L))為后驗信息中抽出的L個樣本。由這些后驗樣本可得到全部的后驗曲線p(d;Δ(l),h(l))，l=1,…,L，劑量反應關系的后驗估計為：

(6)

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

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

(7)

BMD可由后驗估計均值得到：

(8)

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

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

假定ω(d)是關于d的標準布朗運動，則有ω(d1)=0，且當di

(9)

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

劑量反應關系的模型為：

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

(10)

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

(11)

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

(12)

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

(13)

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

模擬研究

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

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

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

表1 模擬研究結果

2.模擬結果

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

實例分析

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

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

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

2.實例分析結果

對于每個數(shù)據(jù)集，都采用BMDS軟件中包含的9種二分反應模型[15]，以及NPB1和NPB2分別來估計劑量反應關系，BMD和BMDL值見表3。

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

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

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

討 論

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

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

在實際應用中，當樣本量很大(N=1000)時相應的模型選擇的正確率較高，可根據(jù)參數(shù)模型估計相應的BMD和BMDL值，但BMDS軟件中包含的參數(shù)模型較多，需要用各個模型分別擬合數(shù)據(jù)，后根據(jù)模型擬合優(yōu)度檢驗結果及AIC原則等模型選擇方法選出最優(yōu)模型，整個計算過程比較繁瑣，還要承擔模型錯誤指定的風險，因此，推薦使用本課題提到的非參數(shù)貝葉斯方法。雖然兩非參數(shù)貝葉斯方法原理部分較難理解，但其相應的R程序均已編寫完成，使得應用非常方便、省時，結果卻更加保守可靠。而NPB2方法與NPB1方法相比，在擬合結果精度和軟件運算速度方面優(yōu)勢更明顯，更推薦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]張俊杰,張海燕.基準劑量評估系統(tǒng)的研究與實現(xiàn).北京林業(yè)大學,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]顧劉金,陳瓊姜,吳立仁,等.基準劑量統(tǒng)計軟件BMDS簡介.中國衛(wèi)生統(tǒng)計，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

國家自然科學基金(81473073)

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

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

郭海強)