BKMR模型实现混合物分析

四叶草^0^

已于 2023-06-12 14:39:17 修改

阅读量1.7w

点赞数 33

文章标签： r语言数学建模

于 2023-06-12 14:37:32 首次发布

原文链接：https://www.bilibili.com/read/cv22070728

版权

文章介绍了如何运用BKMR模型，一种贝叶斯核函数回归方法，来研究环境流行病学中的混合物分析，特别是针对重金属暴露与不孕症关联的研究。通过R语言的bkmr包，作者展示了数据预处理、模型拟合和结果解释的过程，揭示了不同重金属对不孕症的潜在影响及可能的交互作用。

摘要由CSDN通过智能技术生成

混合物分析时近几年环境流行病学最热门的方法，很多污染物是同源的，且存在一定的交互作用。在传统的统计分析方法当中，因变量之间存在高度的共线性关系、复杂的非线性关系以及暴露变量的交互作用导致回归模型的拟合效果不佳甚至结果错误。

BKMR模型介绍

BKMR又称贝叶斯核函数回归，是哈佛大学J.F. Bobb等2015年在Biostatistics上发表的方法，该团队后续在2018年发布BKMR包（R）可被调用处理应用问题。其数学知识包括了回归分析层、核函数层以及变量选择三个部分；BKMR无需设置参数表达形式，允许非线性效应和交互作用的存在，根据放进模型的暴露变量生成核函数，再利用贝叶斯抽样和分析方法，生成混合物组分和暴露变量的关系曲线（dose-response curves），公式如下：

$Y_{i}=h(z_{i1},...,z_{iM})+x_{i}\beta +\epsilon _{i}$

其中， $Y$ 为个体 $i$ 的事件， $Z$ 为暴露因素， $X$ 为协变量； $z_{iM}$ :第 $M$ 个暴露因素， $h()$ :暴露响应函数， $\beta$ ：对协变量的效应， $\epsilon$ ：残差项。

BKMR模型R语言实战

为研究金属暴露与不孕症之间的关联，通过NHANES数据库挖掘了五年的研究人群资料数据如下

协变量：RIDAGEYR,INDFMPIR......PAQ600年龄、性别、是否喝酒等

暴露因素变量（重金属）：URXUCD,URXUCO......URXUSB,URXUSN镉、钴、锡等

结局变量：Infertility是否患有不孕症

数据准备

#加载相关包以及数据
rm(list=ls())
setwd("F:/NHANS")
pacman::p_load(bkmr,readxl,ggplot2)
dat = readxl::read_excel('data_bmkr.xlsx',sheet=1,col_names = TRUE)
#将分类变量转换为因子型
dat$DMDEDUC2 = factor(dat$DMDEDUC2)
dat$RHQ078 = factor(dat$RHQ078)
dat$RHQ031 = factor(dat$RHQ031)
dat$PAQ600 = factor(dat$PAQ600)
划分协变量、暴露变量以及结局变量
covar = data.matrix(dat[,c('RIDAGEYR','INDFMPIR','BMXBMI','DMDEDUC2','RHQ078','RHQ031','RHQ010','PAQ600')])
expos = data.matrix(dat[,c('URXUCD','URXUCO','URXUBA','URXUCS',"URXUTL","URXUPB","URXUMO","URXUMN","URXUSB","URXUSN")])
Y = dat$Infertility
#标准化暴露因素
scale_expos = scale(expos)

kmbayes拟合BKMR模型

set.seed(20230404)
fitkm = kmbayes(Y,Z=scale_expos,X=covar,iter=1000,verbose = FALSE,varsel = TRUE,family='binomial',est.h = TRUE)
TracePlot(fit = fitkm,par = 'beta')
TracePlot(fit = fitkm, par = "sigsq.eps")
TracePlot(fit = fitkm, par = "r", comp = 1)

模型训练迭代次数为1000次（需运行十几分钟）；varsel为True，则说明模型进行了变量筛选

多污染模型

即重金属与不孕症之间的总体关联度：将所选择的重金属暴露因素的总体浓度视为因变量，探究其对不孕症的发生之间的影响。

risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
risks.overall
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()

结果解释：整体的影响趋势为先减弱后增强，且整体为抑制性影响。

单污染模型

即研究单个金属变量的浓度对结局不孕症的影响：计算单个微量元素在第25 - 75百分位数的不孕症潜在连续结局变化，并将其他元素固定在第25、50和75百分位数。估计值为零意味着无效、EST被定义为单个金属元素和潜在的连续结果之间的关联。当其他成员被设定为第75、50和25百分位数时，改金属元素对不孕症产生的影响。

risks.sigvar = SingVarRiskSummaries(
  fit = fitkm, qs.diff = c(0.25,0.75),
  q.fixed = c(0.25,0.5,0.75))
subset(risks.sigvar,variable %in% c('URXUPB','URXUTL','URXUSB','URXUCR','URXUCD','URXUCO'))
ggplot(risks.sigvar,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()

结果解释：HG、CS、CD金属与不孕症的发生是负相关关系，但结果不显著

单变量暴露-反应关系

即当其他金属元素固定在第50百分位数时，估计单元素暴露与流产之间的单变量暴露-反应函数（估计值和95%CI）

pred.resp.univar = PredictorResponseUnivar(fit = fitkm)
ggplot(pred.resp.univar,aes(z,est,ymin=est-1.96*se,ymax=est+1.96*se))+
  geom_smooth(stat = 'identity')+
  facet_wrap(~variable,ncol=4)+
  xlab('expos')+
  ylab('h(expos)')

结果解释：当其他金属元素固定在第50百分位数时，Co对不孕症的促进作用逐渐增大，随后减弱至0附近

双变量-交互作用

每对金属混合物与不孕症的二元关联，由贝叶斯核机器回归估计。它显示了当行中的金属混合物固定在其第25、50和75百分位数，其余金属固定在其第50百分位数时，列中的金属与不孕症之间的关系。

expos.pairs = subset(data.frame(expand.grid(expos1=c(1,2,4,5,6),expos2=c(1,2,4,5,6))),expos1<expos2)
expos.pairs
pred.resp.bivar = PredictorResponseBivar(fit=fitkm,min.plot.dist = 0.5,z.pairs = expos.pairs)
ggplot(pred.resp.bivar,aes(z1,z2,fill=est))+
  geom_raster()+
  facet_grid(variable2~variable1)+
  scale_fill_gradientn(colours = c('#0000FFFF','#FFFFFFFF','#FF0000FF'))+
  xlab('expos1')+
  ylab('expos2')+
  ggtitle('h(expos1,expos2)')
#选0.25、0.5、0.75 4分位数的双变量CR
pred.resp.bivar.levels = PredictorResponseBivarLevels(
  pred.resp.bivar,scale_expos,qs=c(0.25,0.5,0.75))
ggplot(pred.resp.bivar.levels,aes(z1,est))+
  geom_smooth(aes(col=quantile),stat='identity')+
  facet_grid(variable2~variable1)+
  ggtitle('h(expos1|quantiles of expos2)')+
  xlab('expos1')