江蘇省計劃生育科學技術研究所(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。
使用SAS9.3軟件隨機分配各節育器對應的4種有關聯的不良事件,按照構建的模型和參數設置節育器與不良事件組合發生的頻數。模擬產生100個對應數據集,信號檢測結果取100個數據集結果的均值。
提取信號檢測計算時所需要的四格表數據,使用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)
將各種方法的檢測結果與模擬設置的真實信號(RR不為1的組合)進行比較,計算靈敏度、特異度、陽性預測值、陰性預測值、約登指數等指標,并比較不同判定閾值A、不同RR對各方法檢測結果的影響。
100個模擬數據集共有30792825條記錄,宮內節育器與不良事件組合96312個,平均每個數據集有30余萬條記錄,涉及900余個組合。模擬數據中RR的分布如表1所示,各組合記錄數的分布如表2所示。

表1 模擬數據庫中宮內節育器與不良事件組合RR的分布

表2 模擬數據庫中宮內節育器-不良事件組合記錄數的分布
考慮到在模擬數據庫參數設置中,將RR分為四級,可針對不同RR分析各檢測方法的靈敏度,結果如表3所示。隨RR增大,各方法的靈敏度均明顯增加。在同一RR下,PRR和ROR法的靈敏度非常接近,而BCPNN和GPS法的靈敏度稍低;基于FDR控制的四種方法相較傳統方法靈敏度均有所提高,但BCPNN和GPS法的提高幅度大于PRR和ROR法。

表3 各檢測方法在不同RR下的靈敏度
模擬數據庫中宮內節育器與不良事件組合記錄數的最小值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法為首選。
模擬研究存在其自身局限性,真實世界中不良反應/事件報告的影響因素眾多,不良事件背景發生率、報告概率等都未可知,很難找到完全理想的模型來構建。本次模擬的數據與真實的宮內節育器不良事件報告數據相比,節育器與不良事件的組合數稍多,各組合涉及報告數的分布較為接近,結果可供相關研究人員參考。后續可更詳盡地考慮其他影響因素,設置合理的模型和參數,開展更為充分的研究。