重磅：暴露组学混合物分析abkmr包,10万数据1小时搞定,取代bkmr包
背景：
混合物分析越来越热门了，尤其是现在组学当道,新方法发高分是硬核需求。bmkr就是切中要害的新方法。
bkmr可以解决暴露组学或营养、环境、生化指标多组学数据复杂关系的分析，兼具可解释性，获取总体效应及单独效应和重要性，是一个不可多得的创新性分析思路，让你在旧数据上开新花。 实战医学统计3年在国内首次系统介绍了BKMR， 混合物分析新思路-BKMR套路解析 , BKMR死机了？20倍速度bkmrhat包试一试, 混合物分析前沿方法-5维度与工具,好几篇都阅读量不知不觉间破万了，听闻许许多多的硕士博士SCI都借鉴了这几篇文章。后续还介绍了WQS、Qgc、cvek方法，有缘者见了。感兴趣的可以点红色链接去回看。
贝叶斯核机器回归（BKMR）是由哈佛大学的Brent Coull教授团队2015年提出，2018年编写了bkmr包。该团队发表了众多mixture analysis的方法学及论著。 不过bmkr包在使用的时候有个最大的毛病，样本量超过1000，就开始卡，如果去到1万电脑不好的就需要运行10多个小时，样本量为10万的数据集是完全不可行的。 超大样本量如何用bmkr，一直是个大问题？今天热辣出炉的文章就是介绍一个国人开发的新R包，aBKMR，解决卡死问题，效率是bkmr包的70万倍！也就是说潜力无限，可以适配暴露组学。究竟如何使用，follow me。
一、算法
其实我一直都想做这样的事情，py在ML、轨迹领域运行效率更高，而R又可以调用py，为何不修改底层代码让他快速？群里小伙伴分享后拿到文章一看，发现作者思路居然也是这样！那索性就正好好好深入来学习，把心得也分享一下。本文下载地址，https://www.cell.com/the-innovation/fulltext/S2666-6758(25)00451-5 ，bkmr原方法地址https://academic.oup.com/biostatistics/article/16/3/493/269719?login=false。
1BKMR算法
我们先来看最熟悉的bkmr，知己知彼百战百胜。具体方法学文章，背景链接处 混合物分析新思路-BKMR套路解析 。Bkmr是解决传统回归模型在处理多重共线性、非线性剂量-反应关系以及复杂交互作用时的局限性。
2 aBKMR算法
拆解源码，核心算法与bkmr一模一样，但是加速核心在于矩阵运行用py取代R，且全采样变为择优节点。用近似快速换精度耗时。本质没有严格的统计学方法原创创新，毕竟15年的原文方法学惊为天人，很难超越。
3 算法差异与为什么这么快
abmkr可以运行超过70万倍速度，核心在于bkmr旧版R包的单线程问题，py是多线程高效率。下面我从2个py和R的替代视角来解读，为什么能快这么多。下面表格是我总结的。
3.1 矩阵运行从R→ py
原始 bkmr: 所有的矩阵运算（求逆 solve、Cholesky分解 chol、矩阵乘法) 都在 R 内部完成。当样本量较大时（例如N>500），R 的单线程矩阵运算速度会成为瓶颈，也就是大家平时说的卡死了。abkmr: 引入了 reticulate 包，不会在Rstudio中用py的同学看我推文，Rstudio完美运行python。
通过 load_py_function.R 加载了 numpy 和 scipy。它将最耗时的步骤（协方差矩阵的逆）封装在 Python 函数中。例如 a_sigsq.eps.update 中，原始代码计算 crossprod(mu, Vinv) %*% mu，而 abkmr 调用 py$compute_muVinv。目的: 利用 Python NumPy 库底层的优化（通常链接到 BLAS/LAPACK 库）来加速高维矩阵运算。这个过程py发挥了比R运算高效碾压优势。
3.2 节点选择从全秩→固定优选
原始 bkmr: 当使用高斯预测过程（Gaussian Predictive Process, GPP）来处理大数据时，通常需要选择“节点” (Knots)。原始包通常建议用户随机抽取子集或使用均匀网格。abkmr: 新增了 sam_py_r 函数。该算法不再随机抽取，而是计算样本之间的距离，通过最小化几何空间填充标准（geometric space-filling criterion）来选择最具代表性的节点。随机抽样“卡” 是因为它的空间利用率低，迫使你使用成百上千个节点，导致矩阵计算量爆炸。
举个例子，BKMR 模型的核心计算量在于对核矩阵（Kernel Matrix）求逆或分解。如果节点数量为k，这个计算的时间复杂度约为 O(K^3)。也就是说，节点数量翻一倍，计算时间会变成原来的 8 倍。随机选得不好必须增加节点数导致矩阵运算量爆炸，MCMC 跑不动，恶性循环，就出现“卡”住了的问题。
abkmr的加权节点采样思路非常棒！它像铺地砖一样，计算样本之间的距离，特意选择那些能最大化彼此距离或均匀撑开整个数据空间的点。这能在使用较少节点（降低计算维度）的情况下，尽可能保留原始数据的分布信息，从而提高近似计算的精度。
3.3 算法小结
总结就是abkmr是 bkmr 的 "高性能计算实现版"。
他完整复用了 bkmr 的统计模型公式、先验分布设定和 MCMC 采样流程，加速核心是将耗时的线性代数运算从 R 搬到了 Python (NumPy)。改进了在大数据近似计算时选择“节点”的算法（加权采样），以精度换速度。
尽管高斯近似过程 (GPP) 是 Banerjee 就在 2008 年提出的经典方法；加权节点采样是基于经典的“空间填充设计 (Space-filling design)”算法；但是，abkmr 的贡献在于将这些复杂的算法从 R 移植到了 Python 中进行加速，并针对 BKMR 的数据特性（如重复观测值）设计了具体的“加权”逻辑。这点还是十分值得点赞，造福科研人员。
二、bkmr案例分析
官网案例R代码图片和代码如下。
# BKMR / aBKMR 示例流程（模拟数据 -> 拟合 -> PIP -> 总体效应 -> 单变量/交互）
# 说明：
# - 使用 aBKMR 提供的近似方法（knot + NN）加速 BKMR 拟合
# - 展示 PIP（变量选择概率）、总体混合暴露效应、单变量暴露-反应、双变量交互等
# install.packages("devtools")
devtools::install_github("Guo-yi-y/A-BKMR")
# 1) 加载包
library(aBKMR)
library(bkmr)
library(ggplot2)
library(dplyr)
library(glue)
library(ggtext)
# 2) 模拟数据
set.seed(2026) # 随机种子：保证每次运行结果可复现
# sim_data() 生成：
# - y：结局（连续型，gaussian）
# - Z：暴露矩阵（n x M）
# - X：协变量矩阵/数据框
# 参数说明（取决于 aBKMR 的 sim_data 实现）：
# n=500：样本量
# M=5：暴露个数（z1~z5）
# ind=1:3：可能表示对前3个暴露施加特定关联/非线性结构
# Zgen="realistic"：暴露分布更“真实”
# family="gaussian"：连续型结局
dat <- sim_data(n = 500, M = 5, ind = 1:3, Zgen = "realistic", family = "gaussian")
#dat$y：长度为 n 的向量（结局）
#dat$Z：n × M 的矩阵（暴露）
#dat$X：n × p 的矩阵/数据框（协变量
df <- data.frame(
y = dat$y,
dat$X,
dat$Z
)
str(df)
# 3) 选择 knots（用于近似 BKMR，提高速度）
# sam_py_r() 说明（你的注释已写，这里补充中文要点）：
# R: 全部暴露数据矩阵（dat$Z）
# nd: 在暴露空间中抽取多少个候选点/样本用于选 knot（经验上 50~200 常见；越大越慢）
# num_nn: 近邻数，用于寻找“代表性” knots（建议 ≤ 100；越大越慢）
# w: 是否使用权重（当 R 中存在重复暴露观测时可设 TRUE）
#
# 输出 k_id：被选作 knots 的行索引
k_id <- sam_py_r(R = dat$Z, nd = 100, num_nn = 100, w = FALSE)
# 4) 拟合 aBKMR 模型
# a_kmbayes() 关键参数说明：
# y: 结局
# Z: 暴露矩阵
# X: 协变量（可为空；若含分类变量建议提前 dummy 化）
# knots: 指定 knots（加速近似）
# est.h = TRUE: 一定要 TRUE，很多后续的暴露-反应/总体效应函数需要 h 的估计
# varsel = TRUE: 开启变量选择（得到 PIP）
# iter: MCMC 迭代次数（推荐 10000；真实研究通常要更长并做收敛诊断；我测试1000）
# family="gaussian": 连续型结局
#
# 注意：
# - 若 X 含分类变量：可用 model.matrix 转为哑变量（见下方可选代码）
km <- a_kmbayes(
y = dat$y,
Z = dat$Z,
X = dat$X,
knots = dat$Z[k_id, ],
est.h = TRUE,
varsel = TRUE,
iter = 1000,
family = "gaussian"
)
# 若 X 包含分类变量，可选做法（示例）：
# X <- model.matrix(~ . , data = X) %>% .[, -1]
# 5) 提取并绘制 PIP（后验包含概率）
# ExtractPIPs(km) 返回每个暴露的 PIP
# PIP 越大，说明该暴露更可能与结局相关（在 varsel 框架下）
PIP_new <- ExtractPIPs(km) %>%
arrange(desc(PIP))
# 绘制 PIP 柱状图
ggplot() +
geom_col(
data = PIP_new,
aes(x = reorder(variable, PIP, decreasing = TRUE), y = PIP),
fill = "#48C5B9"
) +
ylab("Estimated PIPs") +
xlab("Exposures") +
theme_bw()
# 6) 总体混合暴露效应（Overall risk summary）
# 思路：
# - 构造 newz：给每个暴露设定一串“共同变化”的分位点（如 25% -> 75%）
# - a_OverallRiskSummaries_vary() 计算当混合暴露整体从低到高移动时，结局变化
# - com_je() 用于给出 Joint effect（联合效应）的点估计和区间
# 6.1 构造分位点序列（25% 到 75%，步长 0.01）
# 说明：你也可以改成更宽（0.1）以减少点数、加快绘图
quants <- seq(0.25, 0.75, 0.01)
# 6.2 为每个暴露生成对应分位点值，并拼成 newz（每列一个暴露）
# newz 的每一行表示“某个整体暴露水平”，暴露共同从低 -> 高
newz <- lapply(1:ncol(dat$Z), function(i) {
qf <- data.frame(X = quantile(dat$Z[, i], quants))
names(qf) <- names(data.frame(dat$Z))[i]
qf
}) %>%
bind_cols()
# 6.3 计算总体风险汇总
# data.comps = km$data.comps：aBKMR 对数据结构/标准化信息的缓存，后续函数常需要
overall_res <- a_OverallRiskSummaries_vary(km, newz = newz, data.comps = km$data.comps)
# 给绘图数据加回 quants（对应横轴）
overall_res$risk_overall$quants <- quants
# 6.4 计算联合效应（JE）
# 注意：你这里把 quants*100，用意是“每 1% 增加”的尺度表达
# 实际含义取决于 com_je() 的定义（通常是 X 的尺度）
je_res <- com_je(
overall_res$preds_e$postmean,
overall_res$preds_e$postvar,
X = quants * 100
)
# 6.5 绘制总体效应曲线（点+置信区间）
ggplot() +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
geom_pointrange(
data = overall_res$risk_overall,
aes(
x = quants,
y = est_p,
ymin = est_p - 1.96 * sd_p,
ymax = est_p + 1.96 * sd_p
),
color = "#48C5B9",
shape = 16
) +
ylab("Relative risk of outcome") +
xlab("Quantiles of exposures") +
annotate(
"text",
x = 0.5,
y = 0.4,
size = 4,
label = glue(
"Joint effect = {sprintf('%.3f', je_res$est)} ",
"({sprintf('%.3f', je_res$lower_quant)}, {sprintf('%.3f', je_res$upper_quant)})"
)
) +
theme_bw()
# 7) “指定暴露变化方向”的总体效应（例：z1 上升、z2 下降，其它固定）
# 说明：
# - 有时你想模拟“暴露谱的某种替换/对冲”：
# z1 从 25%->75%，同时 z2 从 75%->25%，其他暴露固定在中位数
newz <- data.frame(
z1 = quantile(dat$Z[, 1], seq(0.25, 0.75, 0.01)),
z2 = quantile(dat$Z[, 2], seq(0.75, 0.25, -0.01)),
z3 = median(dat$Z[, 3]),
z4 = median(dat$Z[, 4]),
z5 = median(dat$Z[, 5])
)
overall_res <- a_OverallRiskSummaries_vary(km, newz = newz, data.comps = km$data.comps)
overall_res$risk_overall$quants <- seq(0.25, 0.75, 0.01)
je_res <- com_je(
overall_res$preds_e$postmean,
overall_res$preds_e$postvar,
X = seq(0.25, 0.75, 0.01) * 100
)
ggplot() +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
geom_pointrange(
data = overall_res$risk_overall,
aes(
x = quants,
y = est_p,
ymin = est_p - 1.96 * sd_p,
ymax = est_p + 1.96 * sd_p
),
color = "#48C5B9",
shape = 16
) +
annotate(
"text",
x = 0.5,
y = 0.4,
size = 4,
label = glue(
"Joint effect = {sprintf('%.3f', je_res$est)} ",
"({sprintf('%.3f', je_res$lower_quant)}, {sprintf('%.3f', je_res$upper_quant)})"
)
) +
theme_bw()
# 8) 单变量暴露-反应曲线（Univariate Predictor-Response）
# 说明：
# - a_PredictorResponseUnivar()：在其他暴露固定（通常在中位数/均值）时，
# 让某一个暴露在分位点上变化，计算结局变化（近似法 method="approx"）
# - 结果里常包含每个暴露的曲线，以及一个“该暴露的总体效应估计”摘要（est）
singlevar_res <- a_PredictorResponseUnivar(
km,
quants = seq(0, 1, 0.1),
method = "approx",
data.comps = km$data.comps
)
# 整理成可 facet 的长数据，并把每个暴露的“效应估计+区间”拼到 facet 标题上
singlevar_data <- lapply(seq_along(singlevar_res), function(i) {
singlevar_res[[i]]$risk_overall %>%
mutate(
var_est = glue(
"{vars} [",
"{sprintf('%.3f', singlevar_res[[i]]$est$est)} ",
"({sprintf('%.3f', singlevar_res[[i]]$est$lower_quant)}, ",
"{sprintf('%.3f', singlevar_res[[i]]$est$upper_quant)})",
" ]"
)
)
}) %>%
bind_rows()
# 绘制单变量曲线（这里用 geom_smooth + stat="identity" 画已计算好的曲线）
ggplot() +
geom_smooth(
data = singlevar_data,
aes(x = z, y = est, ymin = est - 1.96 * se, ymax = est + 1.96 * se),
color = "#48C5B9",
stat = "identity"
) +
facet_wrap(~var_est) +
ylab("Relative risk of outcome") +
xlab("Exposures")
# 9) 双变量交互（Bivariate Predictor-Response）
# 说明：
# - a_PredictorResponseBivar() 先得到两两暴露的响应面信息
# - PredictorResponseBivarLevels() 把“另一个暴露”固定在特定分位点（如 25%、75%）
# - com_2interaction() 计算交互检验（P_interaction），你这里针对 z1-z2
inter <- a_PredictorResponseBivar(
fit = km,
data.comps = km$data.comps,
min.plot.dist = 1
)
# 设定 qs = c(0.25, 0.75)：让 expos2 分别取 25% 和 75% 水平，画 expos1 的曲线
inter2 <- PredictorResponseBivarLevels(inter, Z = dat$Z, qs = c(0.25, 0.75))
# 计算 z1 与 z2 的交互检验结果（返回 P 值等）
inter2_res <- com_2interaction("z1", "z2", inter2)
# 生成 facet 内标注文本（使用 ggtext::geom_richtext 支持下标/斜体）
annotations <- data.frame(
variable1 = c("z1", "z2"),
variable2 = c("z2", "z1"),
x = c(0, 0),
y = c(-1, -1),
label = c(
glue("<i>P</i><sub>interaction</sub> = {sprintf('%.3f', inter2_res$P[1])}"),
glue("<i>P</i><sub>interaction</sub> = {sprintf('%.3f', inter2_res$P[2])}")
)
)
# 只画 z1-z2 相关面板
ggplot(
inter2 %>% filter(variable1 %in% c("z1", "z2") & variable2 %in% c("z1", "z2")),
aes(z1, est)
) +
geom_smooth(aes(col = quantile), stat = "identity") +
facet_grid(variable2 ~ variable1) +
ggtitle("h(expos1 | quantiles of expos2)") +
ylab("Relative risk of outcome") +
xlab("expos1") +
geom_richtext(
data = annotations,
aes(x = x, y = y, label = label),
inherit.aes = FALSE,
fill = NA,
label.color = NA,
color = "red",
size = 5,
hjust = 0,
vjust = 0
)
# 10) 单变量效应在“其它暴露固定水平”下的差异（分层/趋势检验）
# 说明：
# - a_SingVarRiskSummaries()：
# 对某个暴露，在其它暴露固定为不同分位点（q.fixed）时，
# 计算该暴露从 qs.diff[1] -> qs.diff[2] 的效应（类似“低->高”的差值效应）
# - trend()：你这里对每个 q.fixed 的效应做趋势检验，得到 P_trend
multi_inter_res <- a_SingVarRiskSummaries(
km,
y = dat$y,
Z = dat$Z,
data.comps = km$data.comps,
qs.diff = c(0.25, 0.75), # 暴露差：25% -> 75%
q.fixed = c(0.25, 0.50, 0.75),# 其它暴露固定在这三档
method = "approx"
)
# 只取 z1、z2 的结果（示例）
multi_2inter <- multi_inter_res %>% filter(variable %in% c("z1", "z2"))
# 对 z1（前3行）与 z2（后3行）分别做趋势检验
# 这里假设 multi_2inter 的行顺序恰好是 z1 的3档在前，z2 的3档在后
multi_inter1 <- sprintf("%.3f", trend(multi_2inter$est_p[1:3], multi_2inter$sd_p[1:3])$p_trend)
multi_inter2 <- sprintf("%.3f", trend(multi_2inter$est_p[4:6], multi_2inter$sd_p[4:6])$p_trend)
ggplot(
multi_2inter,
aes(
x = variable,
y = est_p,
ymin = est_p - 1.96 * sd_p,
ymax = est_p + 1.96 * sd_p,
col = q.fixed
)
) +
geom_pointrange(position = position_dodge(width = 0.75)) +
coord_flip() +
annotate(
"text",
x = 2.4,
y = 0.8,
size = 6,
label = glue("~italic(P)[trend] == '{multi_inter1}' "),
parse = TRUE
) +
annotate(
"text",
x = 1.4,
y = 0.7,
size = 6,
label = glue("~italic(P)[trend] == '{multi_inter2}'"),
parse = TRUE
)
三、发散思维
abkmr是为了解决bkmr“跑不动”的问题，它牺牲了完美的精度（通过优化的节点选择将损失降到最低），换取了处理大规模数据的能力。统计算法几乎都是bkmr的底层算法， 并不算统计方法创新，要知道bkmr横空出世的时候，是从0到1的创造。这里必须要强调申明一下。下载地址，https://github.com/Guo-yi-y/A-BKMR。
abkmr核心的、计算量大的矩阵运算仅仅是从R过度到py，通过 reticulate 接口外包给了 Python (NumPy) 执行，并增加了新的节点采样策略。我测试下来的确速度飞快，但是时间有限，后续再深度使用再看，目前速度和美观性都有大幅度提升。abkmr额外完善了图片的p值，包括线性趋势检验判断单变量暴露-反应关系或总体联合效应是否呈现显著的线性上升或下降趋势、双变量交互作用的统计检验比较变化率是否显著不同。
abkmr能更好广泛适用于环境混合暴露评估、营养摄入研究、生活方式分析等多个领域，并具备拓展至重复测量资料、中介分析及分布滞后非线性模型（DLNM）等复杂场景的潜力。但是问题也存在，以精度换速度，在小样本<500的时候可能存在偏倚问题，>1万人的时候abkmr的确是最优解。毕竟样本量在那，牺牲写精度可以理解。
py取代R， 方法我之前推文介绍了。Rstudio完美运行python。这步的确是py的碾压优势，对于py小白十分受用。由于模型底层算法代码一样，相信这个极速R包可以未来完全取代笨重的bkmr包的运算。发R包新思路告诉你们了啊，底层代码R转py ，下一篇也许就是你的。
是否大家照着这个思路改改其他R包，也可以投innovation？问题来了，这样的改造是否严格算统计学原创？还是算R包原创？我非常困惑 。
文末可以留言一起自由讨论下，请不要对任何人人身攻击，因为这个现象越来越普遍，我借着这个文章觉得值得好好讨论下，毕竟30分的innovation是国产之光。 求各位小伙伴们指点迷津。
四、小结
abkmr是大样本量混合物分析的神器，运行急速，完全替代bkmr，>1万人时强烈推荐??????????。
<500人时更推荐用bkmr，待定？