FAERS不良事件信号挖掘实践：从数据清洗到ROR/PRR计算

简介：这是一份面向药物警戒研究人员、医学数据分析师及临床药理学者的研究型文档，系统展示基于FAERS数据库开展培塞利珠单抗不良事件信号挖掘的完整流程。内容涵盖研究背景、数据清洗、标签标准化、人口学筛选、比率比与95%置信区间计算、信号判定阈值等核心方法，并分系统讨论神经系统、皮肤、消化系统等不良事件信号及亚组分析结果。资源为单个docx文件，整包约93KB，目录层级清晰，便于按模块参考。目前已有62人学习使用，适合需开展真实世界药物安全性研究、撰写相关论文或了解信号检测技术路径的读者借鉴。全文还涉及严重不良事件特征与常见/罕见信号对比，可直接复用于研究设计。

1. FAERS数据库中培塞利珠单抗不良事件信号挖掘：先搞清这是给谁用、要输出什么

FAERS数据库中培塞利珠单抗不良事件信号挖掘，核心不是把统计公式跑一遍，而是从几十万条自发报告里找出真正值得医学团队跟进的风险信号。常见场景是：某药企负责药物警戒的分析师拿到一个季度更新包，业务方指定某个生物制剂，要求出一份“培塞利珠单抗不良事件信号清单”，里面不能只有 P 值和倍数，还要能说清这个信号是真实暴露关系，还是背景医药环境里的噪声。

这件事适合三类人：需要给注册或安全性报告提供证据的警戒分析师、帮临床团队做上市后安全性监测的研究人员，以及想用公开报告数据练手信号挖掘的医药信息从业者。它解决的是一类非常实际的问题——FAERS 报告数量大、字段乱、重复多，直接算出来的信号往往虚高或漏报，必须先做数据治理，再用比例失衡法做初筛，最后用医学逻辑兜底。整条链路里，数据清洗比统计难，结果解释比跑脚本难。

2. 拿到 FAERS 原始季度包：从压缩包到一份可分析的事件表

2.1 FAERS 原始文件的五张核心表与关键字段

FAERS 的数据发布形式是季度压缩包，解压后通常是一批按模块命名的文件，最常见的是 DEMO、DRUG、REAC、INDI、THER、OUTC、RPSR。做不良事件信号挖掘，多数情况只用到 DEMO、DRUG、REAC 三张表，少数场景会看 INDI 适应症和 OUTC 结局，用来做分层分析。

在实际处理中，我先用命令建一个干净的工作目录，把季度包统一解压，然后用 shell 确认文件是否能正常读取，避免把 Excel 导出的“伪 CSV”混进来。这一步很基础，但值得用脚本固化下来，因为季度包文件名和内部文件后缀并不是每个季度都完全一致。

mkdir -p faers_q3 && cd faers_q3 unzip -o faers_ascii_*.zip # Windows 下如果文件是 .txt 或 .dat，先批量改名成 .csv for f in *.txt *.dat; do mv -- "$f" "${f%.*}.csv" 2>/dev/null; done ls -1 DEMO*.csv DRUG*.csv REAC*.csv

解压后不要急着建模。先用wc -l看一下行数，用head看字段名，确认原始列名是大小写混用还是全大写，这会影响后续 R 代码的列名处理。常见的 DEMO 字段包括 primaryid、caseid、caseversion、event_dt、age、age_cod、sex、gndr_cod、reporter_country；DRUG 表里和挖掘直接相关的是 primaryid、drugname、role_cod；REAC 表里最关键的是 primaryid 和 pt，pt 是 MedDRA 首选术语，也就是不良事件的标准表达。

2.2 为什么先做案件去重和角色过滤，再做药物名匹配

FAERS 不是一份“干净的研究型数据集”，它是企业提交报告、患者报告、文献报告汇集的运营库。同一个不良事件可能由多个来源重复上报，也可能在后续季度包里用新的 caseversion 更新。如果不做去重，同一个反应会被重复计数，信号强度会被明显抬高。常规做法是用 caseid 和 caseversion 做版本判断，保留每个 caseid 的最后一个版本。

药物角色过滤同样关键。DRUG 表里的 role_cod 记录了该药物在报告中的角色，常见的有 PS 主怀疑药、SS 次怀疑药、C 伴随药、I 相互作用药。挖掘培塞利珠单抗信号时，如果把伴随药也计入暴露组，那其他药物的不良反应会被错误归因到这个生物制剂上，导致假信号成片出现。我一般会在初筛阶段只保留 PS 和 SS，必要时单独跑一版只看 PS 的敏感性分析。

药物名匹配是另一个容易翻车的地方。FAERS 的 drugname 列同时存在通用名和商品名，培塞利珠单抗的英文写法有 certolizumab、certolizumab pegol，还可能加上制造商名称或剂型后缀。处理逻辑是维护一个目标药物词典，用正则做不区分大小写的匹配，而不是==严格相等。

2.3 用 R 把清洗流程落成可复用脚本

为了可复现，我会把 DEMO、DRUG、REAC 三张表全部读成字符型，避免 R 自动把数字列转成 numeric，导致 primaryid 的小差异被四舍五入。下面是清洗流程的骨架，适合考完数据结构后直接改造成你自己的版本。

library(data.table) dem <- fread("DEMO.csv", colClasses = "character") drug <- fread("DRUG.csv", colClasses = "character") reac <- fread("REAC.csv", colClasses = "character") # 统一小写列名，后续代码才不用猜大小写 setnames(dem, tolower(names(dem))) setnames(drug, tolower(names(drug))) setnames(reac, tolower(names(reac))) # 案件去重：同一个 caseid 只保留 caseversion 最大的记录 dem <- dem[order(caseid, caseversion)] dem <- dem[!duplicated(caseid, fromLast = TRUE)] # 过滤药物角色：只保留主怀疑药和次怀疑药 drug <- drug[role_cod %in% c("PS", "SS")] # 药物名词典，注意顺序和后缀处理 target_dict <- c("certolizumab", "certolizumab pegol") drug[, is_target := grepl(paste(target_dict, collapse = "|"), drugname, ignore.case = TRUE)] # 找出所有包含目标药物的报告号 target_case_ids <- unique(drug[is_target == TRUE, primaryid]) # 保存清洗后的中间结果，供信号计算使用 saveRDS(dem, "dem_clean.rds") saveRDS(drug, "drug_clean.rds") saveRDS(reac, "reac_clean.rds")

这段脚本里需要注意三个点：案件去重用的fromLast = TRUE，可以保留下重复 caseid 的最后一行，也就是版本号更大的那条；药物名匹配用了grepl而没用==，是为了容忍 “certolizumab pegol solution” 这种带剂型后缀的写法；读取时设置colClasses = "character"，是为了保住 primaryid 的原始精度，避免后续 merge 时因类型转换丢数据。

清洗之后，下一步要面对的是“用什么算法算信号”的问题。这一步选错方法，后面的所有产出都会失去说服力。

3. 信号算法选型：ROR、PRR 与贝叶斯方法各自的适用边界

3.1 四格表是信号挖掘的公共底盘

无论是比例失衡法还是贝叶斯方法，最终都要落到同一个四格表上：目标药物暴露组中该事件报告数 a，目标药物暴露组中其他事件报告数 b，背景组中该事件报告数 c，背景组中其他事件报告数 d。

从 FAERS 里计算 a、b、c、d 时，一个容易被忽略的细节是“报告数”到底是按报告去重还是按反应行去重。我的规则是：同一份报告中的同一 PT 只计一次，不同 PT 各计一次。这样能避免某一报告里列了多个相同 PT 导致的重复计数。另一个细节是暴露组的定义，这里的目标药物暴露组是所有含培塞利珠单抗且 role_cod 为 PS 或 SS 的报告，背景组是去掉这些报告之后的全库报告。

3.2 ROR 和 PRR 的阈值与 95% 置信区间

ROR 和 PRR 是最常用的频率学派方法，非常适合作为 FAERS 信号挖掘的初筛器。ROR 的点估计为 ad/(bc)，表示目标药物目标事件相对于背景组的报告失衡程度；PRR 为 (a/(a+b)) / (c/(c+d))，表达目标药组中目标事件占比与背景组占比的比值。

我在项目里常用的判定条件是：ROR 的 95% 置信区间下限大于 1，且目标事件报告数 a 至少为 3。PRR 的常规阈值是 PRR 大于等于 2，卡方值大于等于 4，且 a 大于等于 3。用小样本时需要对四格表做 0.5 校正，也就是每个格值加 0.5 再带入公式，避免分母为 0 导致无穷大信号。

这两个方法简单透明，适合批量跑全库。但它们的局限也明显：当某事件在背景组只有寥寥几例时，ROR 会特别大，置信区间宽得没有意义。这时候就需要后面的贝叶斯方法做二次确认。

3.3 贝叶斯方法：BCPNN 和 EBGM 在稀疏数据里更稳

BCPNN 和 EBGM 都属于贝叶斯收缩方法，逻辑是当观察数很少时，把信号估计向全局背景收缩，避免一两例数据就催生一个强信号。BCPNN 计算的是信息成分 IC，以 2 为底取对数，当 IC 的 95% 置信区间下限大于 0 时提示信号；EBGM 则以几何均值做点估计，常用判定是 EBGM 大于 2 且 EB05 大于 1。

实际项目里，我的做法是把 ROR/PRR 当成全库批量初筛，把 EBGM 或 BCPNN 当成重点信号的复核算子。FAERS 数据规模大但事件分布极稀疏，很多 PT 在目标药物组只出现 1 到 2 次，频率学方法会给这些事件很高的排名，医学审阅时会浪费大量精力。所以选型本身并不矛盾，而是在不同阶段用不同工具。

下面表是我平时在结果表里同时保留的指标，便于在读结果时直接横向比较。

4. 用 R 算 ROR 和 PRR 并生成信号清单：一份可直接改的完整脚本

4.1 构造目标药物与全库背景的四格表

四格表的构造难点在于“同一份报告是否属于目标药物暴露组”这个标志需要从 DRUG 表合并到 REAC 表。REAC 表行数是全库事件行数，先把每个报告编号是否含培塞利珠单抗这个布尔值打上去，再按 PT 聚合。

这一步不要用merge把整个 DRUG 表摊到 REAC 的每一行，否则数据会膨胀到几百万行，内存直接撑爆。正确做法是只提取目标药报告编号集合，然后根据集合判断 REAC 中每个 primaryid 是否命中。

reac <- readRDS("reac_clean.rds") drug <- readRDS("drug_clean.rds") # 目标药物报告编号集合，DRUG 表里已经做好了 is_target target_case_ids <- unique(drug[is_target == TRUE, primaryid]) # 给 REAC 每一行打标：是否属于目标药暴露组 reac[, has_target := primaryid %in% target_case_ids] # 去重：同一报告同一 PT 只保留一条 reac <- unique(reac[, .(primaryid, pt, has_target)]) # 计算每个 PT 在目标组和非目标组的报告数 cnt <- reac[, .(with_target = sum(has_target), no_target = sum(!has_target)), by = pt] # 计算两组的总报告数 n_target <- uniqueN(reac[has_target == TRUE, primaryid]) n_bg <- uniqueN(reac[has_target == FALSE, primaryid]) # 拼出四格表的 a b c d sig <- cnt[, .( a = with_target, c = no_target, b = n_target - with_target, d = n_bg - no_target ), by = pt]

这里先汇总with_target和no_target，再用总报告数反推 b 和 d，比逐行判断更高效。需要注意uniqueN来自 data.table，它对 primaryid 去重后得到的是真正的报告数，不是事件行数。

4.2 ROR 和 PRR 的计算函数与阈值筛选

拿到 a、b、c、d 后，剩下的就是标准的四格表计算。我用 0.5 校正统一处理零单元格，然后计算 ROR 点估计、95% 置信区间、PRR 和卡方值，一次输出结果表。

# 加 0.5 校正，避免分母为 0 sig[, `:=`(a = a + 0.5, b = b + 0.5, c = c + 0.5, d = d + 0.5)] # ROR 与 95% 置信区间 sig[, ror := (a * d) / (b * c)] sig[, ror_se := sqrt(1 / a + 1 / b + 1 / c + 1 / d)] sig[, ror_low := exp(log(ror) - 1.96 * ror_se)] # PRR 与卡方值 sig[, prr := (a / (a + b)) / (c / (c + d))] N_total <- n_target + n_bg sig[, chi_sq := ((abs(a * d - b * c) - N_total / 2) ^ 2 * N_total) / ((a + b) * (c + d) * (a + c) * (b + d))] # 初筛信号标记 sig[, signal_ror := fifelse(a >= 3 & ror_low > 1, 1, 0)] sig[, signal_prr := fifelse(prr >= 2 & chi_sq >= 4 & a >= 3, 1, 0)] # 输出强调信号：ROR 和 PRR 同时为阳性 sig[signal_ror == 1 & signal_prr == 1, ][order(-a, -ror_low), .(pt, a, b, c, d, ror, ror_low, prr, chi_sq)]

这段代码的核心是每个指标都要有明确的阈值来源，不要随手把 ROR 和 PRR 混着用。代码里fifelse是 data.table 的向量化条件判断，比ifelse快很多，处理几十万行的 PT 表不会卡顿。

4.3 结果表怎么组织：从统计列到医学审阅列

统计结果表只是半成品，真正交付给医学团队时，我会再加上事件分类、适应症、报告年份、严重结局等分层信息，让信号清单可以直接进入人工审阅流程。常见组织方式是给每个 PT 输出一组固定列：PT 名称、目标药组报告数、背景组报告数、ROR、PRR、信号标签、首末报告年份、最常报告的适应症。

这里有一个容易踩的坑：不要把 a 很小的信号往顶层推。即使 ROR 再高，a = 1 时置信区间极宽，信号的可信度很低。我在输出前会强制过滤 a ≥ 3，并且在结果表上保留 a 列，方便医学同事一眼看到证据厚度。只有 a 足够大且比例失衡明显的 PT，才值得进入后续敏感性分析。

# 生成给业务方看的核心清单 final_signal <- sig[signal_ror == 1 & signal_prr == 1 & a >= 3, ][order(-ror_low)] write.csv(final_signal, "signal_certolizumab.csv", row.names = FALSE)

最终清单里的 a、b、c、d 列千万别删，它们是医学审阅时判断信号是否真实的重要依据。删掉背景计数只看比例，是最常见的交付失误。

5. FAERS 信号挖掘避坑清单：5 条来自实际数据踩坑的记录

5.1 案件版本重复导致信号虚高

现象：某个 PT 的 a 值比上一季度翻了一倍，但医学上没有任何新的临床证据，信号排名突然冲进前十。原因：季度包里同一 caseid 出现了新版本，旧版本的事件在新版本里又保留了一部分，两个版本同时计入统计，产生了重复计数。解决：在清洗阶段用caseid + caseversion去重，保留每个 caseid 最新版本，并验证去重前后的总报告数差异。

5.2 药物名正则匹配漏掉拼写变体

现象：目标药物组报告数明显少于业务方的预期，很多已知使用培塞利珠单抗的报告没有进入暴露组。原因：drugname 里写了商品名或带剂型的通用名，匹配用的字符串只覆盖了“certolizumab pegol”一种写法。解决：维护一份目标药物别名表，加入商品名、去掉空格和连字符后再做正则匹配，同时用 DRUG 表中的 ingredient 字段做二次校验。

5.3 伴随用药被当成暴露药物

现象：信号清单里出现“感染”风险升高，但拆分数据后发现很多报告是患者同时在用另一款免疫抑制剂，不良事件实际可能是联合用药带来的。原因：DRUG 表里 role_cod 包含伴随用药，而我没有过滤，导致所有含目标药的报告都进了暴露组。解决：初筛只保留 PS 和 SS，后续做敏感性分析时再单独看 PS-only 的结果，两条线对比。

5.4 “未知”事件名称污染 PT 统计

现象：结果表里出现大量Unknown、Not specified之类的 PT，且 a 值不小，信号排名靠前。原因：REAC 表里存在未编码或整理失败的自由文本，被原样当成了 PT 处理。解决：统计前过滤掉不满足 MedDRA 编码规则的 pt，重点检查是否包含空格、标点、纯数字字符串。

5.5 报告日期缺失导致时间分层分析失效

现象：按年份拆分信号时，某一年的事件数暴跌，趋势断档严重。原因：DEMO 表的 event_dt 字段存在大量非标准日期，比如全 0 值、文本日期、只有年份没有月日，直接读入 R 后被解析成 NA。解决：把日期列统一读成字符，用正则解析，提取年月后另存为year列，无法解析的归入“缺失年份”，单独标记而不是直接丢弃。

6. 信号结果怎么验证：敏感性分析、分层审阅和一条固定收尾习惯

信号挖掘不是算出阳性就结束，真正的交付是让结果经得起多次改动参数的考验。我通常在一版信号清单生成后，立刻跑三组敏感性分析：第一组把药物匹配词典从通用名改成通用名加商品名，看 top 信号是否依旧在同一批 PT 里；第二组把角色从 PS+SS 改成仅 PS，看哪些信号的依赖来自次怀疑药；第三组把 a 的过滤阈值从 3 改成 5，看信号强弱排序是否剧烈抖动。

分层审阅时要看的维度包括报告年份、适应症和结局结局。比如某个 PT 在最近两个季度集中出现，而前几年没有，那就比长期零散出现更像真实信号。再比如某个事件集中在“类风湿关节炎”适应症里，而适应症本身是疾病活动相关症状，就要考虑是疾病本身引起的背景事件，而不是药物直接导致。

最终我会把统计信号切回原始报告逐条读一遍。这个习惯非常重要，因为 FAERS 的数值永远只告诉你“有多少报告撞车”，读原文才能判断“这些报告到底是不是同一临床情况”。比如同一个 PT 名称下面的报告，有的可能是同一患者被重复报告，有的可能是给药剂量和合并用药完全不同，只靠统计列根本看不见。每次读完原始报告后，我都会把明显重复或信息矛盾的报告挑出来，打上“不建议作为信号证据”的标签，再更新信号清单。

这个方向值得投入。FAERS 公共数据本身就是一条低成本的上市后警戒学习路径，把它清洗成可复用的分析表，再把 ROR、PRR、EBGM 等信号方法固化成脚本，后续只需要换药物词典和角色过滤规则，就能复用到其他生物制剂项目上。我的经验是：先把清洗流程做成黑匣子里的固定步骤，再让统计指标可复现，最后用医学审阅兜底。希望这篇记录能帮你少踩几个坑。

本文还有配套的精品资源，点击获取