混合物分析：BKMR死机了？20倍速度bkmrhat包试一试

bkmrhat包解析

bkmr包有4个缺点，1)没有并行同时运行多个MC链2）没有跨多个MC链的推断3）参数的有限的后验总结4）有限的诊断。其实，bkmrhat就是用future包的多线程写入原来的包，达到经典bkmr包速度的4-24倍(根据CPU核心线程不同不同)，甚至更快。Bkmrhat包地址 https://github.com/kdevick/bkmr-cma。主要函数如下，相当于bkmr增加了paralled功能。

R实战

案例与BKMR案例相同，dataset1数据集(文末获取附件)。结果解释参考BKMR章节，一样。

#安装bkmr包
pacman::p_load(bkmr,readxl,ggplot2)
setwd("D:/聂个人文件/统计/高级统计/BKMR")
dat <- readxl::read_excel("dataset1.xls", sheet = 1, col_names =  TRUE)
head(dat)
 
#数据集中，有 500 行、一个结果变量Y、七个暴露X1到X7和一个协变量Z。为了拟合 BKMR 模型，我们定义了一个暴露矩阵和一个协变量矩阵
#将关注的结局（Y），需要调整的协变量（covar）以及关注的暴露（expos）分别封装进3个矩阵
covar <- data.matrix(dat[, "Z"])
#covar <- data.matrix(df[,c("cor1","cor2","cor3",..., "corn")])
#expos <- data.matrix(df[,c("Mg","Ca","Sr","Ba","Cu","Zn","Fe")]
expos <- data.matrix(dat[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7")])
Y <- dat$Y
#感兴趣的暴露变量用 表示z,整后的协变量用表示x。为了与文献保持一致，相应地重命名曝光和协变量
colnames(covar) <- "x"
colnames(expos) <- paste0("z", 1:ncol(expos))
#ln暴露Z, 可选；　log(x) 为自然对数=ln，log10(x)为10为底对数
#for (i in 1:7) {expos[,i]<-log(expos[,i])}
#标准化暴露变量Z
scale_expos <- scale(expos)

#------ bkmrhat 多核 加速4-10倍！默认直接varsel = TRUE 模式；二阶段法需要手动PIP筛选
set.seed(123)
#bkmrhat :kmbayes_parallel函数,同等工作iter*nchain；时间差不多啊！不变量选择 varsel = FALSE 会慢很多
#---- 预分析，后续chain=4，iter=5000；当chain=20，iter=1000；
#单核对多核CPU默认支持4线程nchains=4，其实也可以8,20等
system.time(kmfit6 <- suppressMessages(kmbayes_parallel(nchains=4, Y, Z = scale_expos, X = covar, iter = 250, verbose = FALSE, varsel = TRUE)))
#用户  系统  流逝 
#0.14  0.01 20.43
system.time(kmfit7 <- suppressMessages(kmbayes_parallel(nchains=20, Y, Z = scale_expos, X = covar, iter = 50, verbose = FALSE, varsel = TRUE)))
#8线程14.20s，20线程7s

#----------- 2.1 多污染模型 联合效应-----------
#----------- 2.1.1  汇总多污染模型 联合效应
#-------kmbayes_combine 合并
fitkmccomb = kmbayes_combine(kmfit6)
summary(fitkmccomb)
#---OverallRiskSummaries + fitkmccomb
mean.difference <- suppressWarnings(OverallRiskSummaries(fit = fitkmccomb, y = fitkmccomb[["y"]], Z = fitkmccomb[["Z"]], X =fitkmccomb[["X"]], 
                                                         qs = seq(0.25, 0.75, by = 0.05), 
                                                         q.fixed = 0.5, method = "exact"))
mean.difference
with(mean.difference, {
  plot(quantile, est, pch=19, ylim=c(min(est - 1.96*sd), max(est + 1.96*sd)), 
       axes=FALSE, ylab= "Mean difference", xlab = "Joint quantile")
  segments(x0=quantile, x1=quantile, y0 = est - 1.96*sd, y1 = est + 1.96*sd)
  abline(h=0)
  axis(1)
  axis(2)
  box(bty='l')
})

#----------- 2.1.2  4条MC chain分别展示，多污染模型 联合效应，黑白版
#-OverallRiskSummaries_parallel + kmfit6
risks.overall <- OverallRiskSummaries_parallel(kmfit6, y = Y, Z = scale_expos, X =covar,
                                      qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "exact")
risks.overall
#与所有暴露都固定为其中值时相比，当所有暴露都处于第 75 个百分点时，结果增加了 2.37 个单位
#天花板效应表明非线性暴露-响应，
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd)) + 
  geom_hline(yintercept = 0, lty = 2, col = "brown") +
  geom_pointrange()

#----------- 2.1.3  4条MC chain分别展示，多污染模型 联合效应，彩色版SD区间
risks.overall = OverallRiskSummaries_parallel(kmfit6, y = Y, Z = scale_expos, X =covar ,qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "exact")
nchains = length(unique(risks.overall$chain))
with(risks.overall, {
  plot.new()
  plot.window(ylim=c(min(est - 1.96*sd), max(est + 1.96*sd)), 
              xlim=c(min(quantile), max(quantile)),
              ylab= "Mean difference", xlab = "Joint quantile")
  for(cch in seq_len(nchains)){
    width = diff(quantile)[1]
    jit = runif(1, -width/5, width/5)
    points(jit+quantile[chain==cch], est[chain==cch], pch=19, col=cch) 
    segments(x0=jit+quantile[chain==cch], x1=jit+quantile[chain==cch], y0 = est[chain==cch] - 1.96*sd[chain==cch], y1 = est[chain==cch] + 1.96*sd[chain==cch], col=cch)
  }
  abline(h=0)
  axis(1)
  axis(2)
  box(bty='l')
  legend("bottom", col=1:nchains, pch=19, lty=1, legend=paste("chain", 1:nchains), bty="n")
})

#----------- 2.1.4  4条MC chain分别展示，多污染模型 联合效应，彩色SE置信区间带版
regfuns_par = PredictorResponseUnivar_parallel(kmfit6, y = Y, Z = scale_expos, X =covar ,qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "exact")
nchains = length(unique(risks.overall$chain))
# single variable
with(regfuns_par[regfuns_par$variable=="z1",], {
  plot.new()
  plot.window(ylim=c(min(est - 1.96*se), max(est + 1.96*se)), 
              xlim=c(min(z), max(z)),
              ylab= "Predicted Y", xlab = "Z")
  pc = c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999")
  pc2 = c("#0000001A", "#E69F001A", "#56B4E91A", "#009E731A", "#F0E4421A", "#0072B21A", "#D55E001A", "#CC79A71A", "#9999991A")
  for(cch in seq_len(nchains)){
    ribbonX = c(z[chain==cch], rev(z[chain==cch]))
    ribbonY = c(est[chain==cch] + 1.96*se[chain==cch], rev(est[chain==cch] - 1.96*se[chain==cch]))
    polygon(x=ribbonX, y = ribbonY, col=pc2[cch], border=NA)
    lines(z[chain==cch], est[chain==cch], pch=19, col=pc[cch]) 
  }
  axis(1)
  axis(2)
  box(bty='l')
  legend("bottom", col=1:nchains, pch=19, lty=1, legend=paste("chain", 1:nchains), bty="n")
})

#---------- 2.2.2 4四条chain
risks.singvar <- SingVarRiskSummaries_parallel(
  kmfit6, y = Y, Z = scale_expos,X= covar,
  qs.diff = c(0.25, 0.75), 
  q.fixed = c(0.25, 0.50, 0.75))
#二步法选12457，色彩图；分4条chain
subset(risks.singvar, variable %in% c("z1", "z2", "z4", "z5", "z7"))
ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd, 
                          col = q.fixed)) + 
  geom_pointrange(position = position_dodge(width = 0.75)) + coord_flip()

#-----------2.3   计算PIP;需要 varsel = TRUE，变量选择
multidiag = kmbayes_diagnose(kmfit6, warmup=0, digits_summary=2)
lapply(kmfit6, function(x) t(ExtractPIPs(x)))
#合并不同起始值的4条链BKMR拟合= fitkmccomb
fitkmccomb = kmbayes_combine(kmfit6)
# 相比fitkm多了2对象， fitkmccomb[["chain"]] fitkmccomb[["iters"]]
summary(fitkmccomb)
ExtractPIPs(fitkmccomb)

Tips

目前未见比较好的文章有详细运用BKMRhat包各项功能的案例，有见到的小伙伴还请文章底部留言告知一下。谢谢。

小结

BKMRhat包提供了相对BKMR包高达20倍以上速度提升，可计算总效应趋势(多污染物模式)、相对重要性(PIP)；也可分MC链展示分链多污染物模式、分链单污染物模式。

数据和代码索要过期不再分享