基于错误发现率控制的不良反应信号检测方法比较*

江苏省计划生育科学技术研究所(210036)

施雯慧 许豪勤△ 巴 磊

【提 要】 目的 通过蒙特卡罗模拟比较基于错误发现率(FDR)控制的信号检测方法与传统方法的检测效能,验证基于错误发现率控制的方法是否能解决传统方法假阳性信号比例高的问题。方法 基于不良反应监测数据的特点,建立自发呈报系统模型,合理设置参数,产生模拟数据,使用R软件进行数据分析,进行传统方法和基于错误发现率控制的方法的比较,并考察判定阈值对信号检测结果的影响。结果 传统频数法(PRR和ROR法)的灵敏度高于贝叶斯法(BCPNN和GPS法),特异度稍低;基于FDR控制的方法较传统方法灵敏度提高,特异度降低,且FDR控制对贝叶斯法的影响较频数法更大。结论 若检测要求较高的灵敏度和较大的ROC曲线下面积,可优先选择基于FDR控制的GPS法,若要求较高的阳性预测值,则传统BCPNN法为首选。

药品不良反应/医疗器械不良事件自发呈报系统目前是国际上药品/医疗器械上市后监测的最主要手段[1-2],自20世纪60年代以来,随着监测体系和相关法规的不断完善,各国的自发呈报系统收集到的不良反应报告数大幅增长,以我国为例,截至2016年国家药品不良反应监测中心已累计收到药品不良反应/事件报告近1075万份[3]、可疑医疗器械不良事件报告近168万份[4]。如何从海量报告数据中及时、准确地发现不良反应信号,是当前药物警戒研究领域的热点问题。目前常用的信号检测方法主要是基于不相称测定原理,可分为：(1)频数法,如比例报告比值比法(proportional reporting ratio,PRR)、报告比值比法(reporting odds ratio,ROR);(2)贝叶斯法,如贝叶斯置信传播神经网络(Bayesian confidence propagation neural network,BCPNN)、经验贝叶斯伽马泊松缩减(empirical Bayes gamma Poisson shrinker,GPS)。这些方法应用已经较为成熟,但存在多重假设检验问题[5]。信号检测过程一般是根据目标药品/医疗器械和目标不良反应/事件,计算各方法下的不相称测定指标(如PRR、ROR、IC等),再依据判定规则判断该目标组合是否为可疑信号。假设数据库中有i种药品/医疗器械和j种不良反应/事件,则相应的检验进行了i×j次,如将每次比较的错误率控制在0.05的水准,则会增大总Ⅰ型错误率,产生较多的假阳性信号。

传统的多重比较方法(如Bonferroni法、sidak法等)控制的是总Ⅰ型错误率[6],如将其定为0.05,那么每次检验的水准就极低,则会导致只能发现少数强关联组合,检验效能降低。FDR(false discovery rate)称为错误发现率或阳性结果错误率,表示阳性检验结果中判断错误的比例,由Benjaminni和Hochberg首先提出[7]。相比传统假设检验的检验水准取值固定,FDR能灵活取值,作为假设检验错误率的控制指标;此外,相比总Ⅰ型错误率主要用于控制I类错误,FDR的意义更为明确,可用于评价筛选出来的差异变量,因而常用于高微阵列数据分析的多重比较[8]。

FDR在药品不良反应信号检测领域的应用始于法国的Ahmed[9-11],但在医疗器械不良事件信号检测方面的效果尚未可知。考虑到信号检测方法没有金标准,本文采用蒙特卡罗模拟的方法来比较四种常用信号检测方法(PRR、ROR、BCPNN、GPS)及基于FDR控制后的四种方法在宫内节育器不良事件数据中的检验效能。

模型构建

Emmanuel Roux等的研究表明[12],自发报告系统的药物暴露频数服从Poisson分布,报告数nij在一定时间Δt内服从δij的Poisson分布:

nij～Poi(δij)=Poi(RRijIjEiPij)

其中RRij为药物i与不良事件j组合的相对危险度,Ij为不良事件j的背景发生率,Ei为药物i的使用人数,Pij为药物i与不良事件j组合的报告概率。

国内已有的几项药品不良反应模拟数据研究中均沿用了上述模型,假定报告影响因素包括用药人数、药品上市时间、不良事件背景发生率、不良事件严重程度、报告概率[13-14]。而国家计划生育药具不良反应监测中心既往几年的监测工作情况表明,宫内节育器不良事件报告并不完全符合上述特点。首先,严重的不良事件报告概率未必高于一般不良事件,因为大部分监测点设立于县、乡两级的妇幼保健计划生育服务机构,能收集到的严重伤害事件不多(由于医疗水平的限制,大多数严重伤害事件会被转诊至三级医疗机构);其次,由于宫内节育器产品的特殊性(使用时间长,研发周期长,在市场上流通种类远低于药品),产品上市时间与报告概率的关系也不是很大;因此,参考相关文献和专家意见,假定宫内节育器不良事件自发呈报系统的报告数服从Poisson分布,影响因素包括使用人数、不良事件背景发生率和报告概率。

参数设置

本次模拟假定有40种宫内节育器和30种不良事件,共计有1200种节育器与不良事件组合,假定每种宫内节育器有4种与其有关联的不良事件(相对危险度RR从低到高设置4档,分别为2、3、5和10),则共计有40×4=160个组合为真实信号,剩余组合为虚假信号,相对危险度RR为1。

40种宫内节育器中,按使用人数分为三类:5种为常用,使用人数估计为500万;15种为一般,使用人数估计为80万;20种为较少使用,使用人数估计为10万。30种不良事件中,按背景发生率分为三类,其中10种为1/100,10种为1/1000,10种为1/20000;报告概率设定为4个水平,分别为0.005、0.025、0.05和0.1。

步 骤

1.生成模拟数据库

使用SAS9.3软件随机分配各节育器对应的4种有关联的不良事件,按照构建的模型和参数设置节育器与不良事件组合发生的频数。模拟产生100个对应数据集,信号检测结果取100个数据集结果的均值。

2.数据提取及信号检测

提取信号检测计算时所需要的四格表数据,使用R软件及PhViD包进行信号检测。PhViD包由Ismail Ahmed开发,使用其内置的主要函数PRR()、ROR()、BCPNN()、GPS(),可根据给定数据以PRR、ROR、BCPNN、GPS四种方法快速计算出信号。这些函数中比较重要的参数包括MIN.n11、DECISION、DECISION.THRES等,其中MIN.n11指定目标药物与目标不良事件组合的最小频数,DECISION指定生成信号的判定规则(基于传统规则还是基于FDR控制),DECISION.THRES指定相应规则的阈值。本软件包中使用的FDR估计方法是基于LBE的估计方法[15]。

根据PhViD包的使用说明,基于传统判定规则生成信号时,对应的函数参数设置如下:

(1)PRR(DATABASE,RR0 = 1,DECISION = 3,DECISION.THRES = 0,RANKSTAT = 2)

(2)ROR(DATABASE,RR0 = 1,DECISION = 3,DECISION.THRES = 0,RANKSTAT = 2)

(3)BCPNN(DATABASE,RR0 = 1,DECISION = 3,DECISION.THRES = 0,RANKSTAT = 2,MC = TRUE,NB.MC = 10000)

(4)GPS(DATABASE,RR0 = 1,DECISION = 3,DECISION.THRES =1,RANKSTAT = 2,TRONC = FALSE,TRONC.THRES = 1,PRIOR.INIT = c(alpha1 = 0.2,beta1 = 0.06,alpha2 = 1.4,beta2 = 1.8,w = 0.1),PRIOR.PARAM = NULL)

基于FDR控制生成信号时,对应的函数参数设置如下:

(1)PRR(DATABASE,RR0 = 1,DECISION.THRES=0.05)

(2)ROR(DATABASE,RR0 = 1,DECISION.THRES=0.05)

(3)BCPNN(DATABASE,RR0 = 1,DECISION.THRES=0.05)

(4)GPS(DATABASE,RR0 = 1,DECISION.THRES=0.05)

3.结果比较

将各种方法的检测结果与模拟设置的真实信号(RR不为1的组合)进行比较,计算灵敏度、特异度、阳性预测值、阴性预测值、约登指数等指标,并比较不同判定阈值A、不同RR对各方法检测结果的影响。

结 果

1.模拟数据描述

100个模拟数据集共有30792825条记录,宫内节育器与不良事件组合96312个,平均每个数据集有30余万条记录,涉及900余个组合。模拟数据中RR的分布如表1所示,各组合记录数的分布如表2所示。

表1 模拟数据库中宫内节育器与不良事件组合RR的分布

表2 模拟数据库中宫内节育器-不良事件组合记录数的分布

2.RR对各方法灵敏度的影响

考虑到在模拟数据库参数设置中,将RR分为四级,可针对不同RR分析各检测方法的灵敏度,结果如表3所示。随RR增大,各方法的灵敏度均明显增加。在同一RR下,PRR和ROR法的灵敏度非常接近,而BCPNN和GPS法的灵敏度稍低;基于FDR控制的四种方法相较传统方法灵敏度均有所提高,但BCPNN和GPS法的提高幅度大于PRR和ROR法。

表3 各检测方法在不同RR下的灵敏度

3.判定阈值A对各方法检测能力的影响

模拟数据库中宫内节育器与不良事件组合记录数的最小值1,最大值50671,四分位数分别为4、25、108,选取A分别为1、2、3、4、10、25时,比较各方法的检测能力。

对模拟数据库中的100个模拟数据集分别计算,求得不同A值下各方法检测出信号数的均值。总体上,各方法生成的信号数随A值增加而减少(BCPNN、GPS法中A≤3时例外);基于FDR控制的4种方法比各自对应的传统方法生成的信号数均有所增加,且BCPNN、GPS法增加的幅度更大;相同A值下,PRR、ROR法生成的信号数近似,基于FDR控制的PRR、ROR法生成信号数也近似,但基于FDR控制的BCPNN和GPS法之间信号数的差距较传统BCPNN和GPS法间的差距缩小。

进一步比较判定阈值A对各方法检测能力的影响,结果表明,随A增大,各方法的灵敏度降低、特异度升高、曲线下面积减少;相同A值下,传统频数法的灵敏度高于传统贝叶斯法,基于FDR控制的频数法灵敏度低于FDR控制的贝叶斯法;四种传统方法经FDR控制后,灵敏度均升高,特异度均降低,且FDR控制对贝叶斯方法的影响较频数法更大,见表4。

讨 论

从模拟结果可以看出,频数法(PRR和ROR法)之间的一致性较高,其灵敏度高于贝叶斯法(BCPNN和GPS法),特异度较之稍低,这与既往研究结论类似。而基于FDR控制的四种方法,检测出的信号数较传统方法增加,灵敏度提高,特异度降低,且基于FDR控制的贝叶斯法灵敏度高于基于FDR控制的频数法,说明FDR控制对贝叶斯法的影响较频数法更大。这一点与采用FDR控制后检测信号数应适当减少的预期不太一致。国内郭晓晶在相关预实验中也发现[16],将基于LBE的FDR估计方法与BCPNN法联用,产生的信号数量多于单独使用BCPNN法,具体原因有待进一步研究。

表4 各检测方法在不同A值下的检测能力

而在不同A值下比较各方法检测的灵敏度、特异度、曲线下面积等,可发现A=1时灵敏度高,曲线下面积较大,这一结果与既往的药品不良反应数据模拟研究并不相同。李婵娟等[13]的研究表明,选择A≥3的组合进行信号检测,灵敏度要高于考虑任何A值的检测结果,并由此得出两份以下的报告对检测意义不大的结论。分析原因,可能是节育器不良事件数据和药品不良反应数据的特点不同造成的。药品不良反应报告数据库中,往往是报告例数小的组合占比较大,如我国2010-2011年药品不良反应报告数据库中报告数小于3例的组合占64.35%[14],广东省2002-2007年的药品不良反应报告数据库中报告数小于3例的组合占76.24%[17],而2011-2015年宫内节育器不良事件数据中报告数小于3例的组合仅占25.03%。

综上所述,信号检测方法应用于宫内节育器不良事件数据中时,没有必要限定于A≥3的组合。若检测目的是要求较高的灵敏度和ROC曲线下面积,可优先选择基于FDR控制的GPS法,若重点关注检测出的信号的可靠性,即要求较高的阳性预测值,则传统BCPNN法为首选。

模拟研究存在其自身局限性,真实世界中不良反应/事件报告的影响因素众多,不良事件背景发生率、报告概率等都未可知,很难找到完全理想的模型来构建。本次模拟的数据与真实的宫内节育器不良事件报告数据相比,节育器与不良事件的组合数稍多,各组合涉及报告数的分布较为接近,结果可供相关研究人员参考。后续可更详尽地考虑其他影响因素,设置合理的模型和参数,开展更为充分的研究。