R语言贝叶斯核机回归用于估算多污染混合物的健康影响

多污染物暴露的健康影响在环境流行病学至关重要。混合物成分可能与健康存在复杂的非线性和非加性关系有时候单用线性来解释可能会捕捉不到阳性结果或者出现错误。核机回归KMR也称为高斯过程回归是机器学习文献中流行的工具。KMR的核心思想是灵活地建模大量变量与特定结果因变量之间的关系。通用建模框架是BKMR 不直接为每个污染物设定回归系数而是将健康结局回归到一个灵活的函数 h(z)上。这个函数通过核函数文中主要使用高斯核来定义不同个体暴露特征Exposure Profile之间的“相似性”。下面咱们使用R语言来演示一下先导入R包library(bkmr)bkmr包自带有一个SimData函数咱们使用它来生成一个模拟数据假设有4个暴露变量50例样本这里默认的是family “gaussian”如果你的分布是其他分布在这里选择对应的dat-SimData(n50,M4)列表中的x是咱们的协变量z是暴露变量y是结局变量咱们把它提取出来y-dat$y Z-dat$Z X-dat$X下面进行马尔可夫链蒙特卡洛MCMCiter表示迭代次数verbose指示是否在模型拟合期间打印中间诊断信息的标志varsel并表示是否对预测变量进行变量选择fitkm-kmbayes(yy,ZZ,XX,iter10000,verboseFALSE,varselTRUE)解析一下模型summary(fitkm)MCMC 方法需要检查模型的收敛性这种围绕中轴的收敛性是比较好的TracePlot(fitfitkm,parbeta)TracePlot(fitfitkm,parsigsq.eps)计算后验概率PIP 值介于 0 到 1 之间越接近 1 越重要通常 PIP 0.5 视为重要。res-ExtractPIPs(fitkm)res变量重要性可视化ggplot(res,aes(xreorder(variable,PIP),yPIP))geom_bar(statidentity,fillsteelblue)geom_hline(yintercept0.5,linetypedashed,colorred)coord_flip()labs(title变量重要性 (PIP),x暴露物,y后验包含概率)theme_minimal()要绘制单变量暴露-反应关系需要先生成预测值生成模型预测值pred.resp.univar-PredictorResponseUnivar(fitfitkm)绘图library(ggplot2)ggplot(pred.resp.univar,aes(z,est,yminest-1.96*se,ymaxest1.96*se))geom_smooth(statidentity)facet_wrap(~variable)ylab(h(z))这里我们看到z1以及z2值的增加与结果值的增加相关Y这是有前提的对于其他预测变量Z固定在第50百分位就是中位数水平协变量在x保持不变。查看两个暴露物之间是否存在交互作用min.plot.dist 1表示网格距离pred.resp.bivar-PredictorResponseBivar(fitfitkm,min.plot.dist1)ggplot(pred.resp.bivar,aes(z1,z2,fillest))geom_raster()facet_grid(variable2~variable1)scale_fill_gradientn(coloursc(#0000FFFF,#FFFFFFFF,#FF0000FF))xlab(expos1)ylab(expos2)ggtitle(h(expos1, expos2))这样有时难看明白咱们可以变换一下subset-pred.resp.bivar%%filter(variable1z1,variable2z2)ggplot(subset,aes(z1,est,coloras.factor(z2)))geom_line()labs(title交互作用分析: Exp1 vs Exp2,subtitle不同颜色代表 Exp2 的不同水平 (分位数),xExp1,yh(z),colorExp2 (Z-score))theme_minimal()如果不同颜色的线不平行即斜率不同则提示存在交互作用。另一种方法是研究单个预测变量的预测-响应函数对于固定在不同分位数的第二个预测变量以及固定在特定值的剩余预测变量。pred.resp.bivar.levels-PredictorResponseBivarLevels(pred.resp.dfpred.resp.bivar,ZZ,qsc(0.1,0.5,0.9))ggplot(pred.resp.bivar.levels,aes(z1,est))geom_smooth(aes(colquantile),statidentity)facet_grid(variable2~variable1)ggtitle(h(expos1 | quantiles of expos2))xlab(expos1)评估混合暴露总体效应就是评估所有暴露物同时增加对结局的影响累积效应是通过比较预测变量的值来计算预测变量的整体效应h当所有预测变量都处于某个百分位而当它们都处于第50百分位时。通常比较所有暴露物都在某一分位数如 25th 到 75th时的结局差异risks.overall-OverallRiskSummaries(fitfitkm,yy,ZZ,XX,qsseq(0.25,0.75,by0.05),q.fixed0.5,methodexact)risks.overall绘图ggplot(risks.overall,aes(quantile,est,yminest-1.96*sd,ymaxest1.96*sd))geom_pointrange(colordarkgreen)geom_hline(yintercept0,linetypedashed)labs(title混合暴露总体效应,subtitle所有暴露物同时增加时的累积效应,x分位数,y效应估计 (相对于中位数))theme_minimal()如果曲线呈现单调上升或下降且置信区间不包含 0说明混合暴露整体对健康有显著影响。