https://mp.weixin.qq.com/s/XZV_5iioPDHnMQfEPCIlMg
BKMR ®
#首先清理缓存。
rm(list = ls())
#运行R包，如果没有下载要先下载。
library(bkmr)
library(ggplot2)
#给数据赋值，如果要自己进行研究，数据的地址以及数据的变量需要对应的自行选择替换，data1是导入数据，covar是定义协变量，expose是定义暴露，y为骨关节炎也就是本次要研究的对象。
data1<-read.csv(“C:\Users\Administrator\Desktop\bkmr.csv”)
covar <- data.matrix(data1[, c(“age”, “RIAGENDR”, “RIDRETH1” , “DMDEDUC2” , “DMDMARTL” , “INDFMPIR” , “BMXBMI” , “alq101” , “LBXCOT”, “phy”)])
expos <- data.matrix(data1[, c(“Ba”, “Cd”, “Co”,“Cs”, “Mo”, “Pb”,“Sb”, “Tl”, “Tu”)])
Y <- data1$OA
#对暴露进行标准化处理
scale_expos <- scale(expos)
#设置种子数以便于文章结果复现，可以设置任何数字，相同数字复现结果一致。
set.seed(1000)
#利用fields包中的cover.design函数，基于标准化后的变量scale_expos生成50个节点的设计矩阵。
knots50 <- fields::cover.design(scale_expos, nd = 50)$design
#这段代码是贝叶斯核机回归的本体，使用kmbayes函数对数据进行基于贝叶斯的回归模型拟合，其中Y是响应变量，Z和X是解释变量，iter表示迭代次数，family表示使用二项式分布的模型，est.h表示估计超参数，verbose表示是否输出详细信息，varsel表示是否进行变量选择，knots表示指定的节点设计矩阵。
fitkm <- kmbayes(Y, Z = scale_expos, X = covar, iter = 10000, family = “binomial”, est.h = TRUE, verbose = FALSE, varsel = TRUE,knots = knots50)
#绘图，分别绘制参数β，误差方差以及变量之间相关系数的跟踪图。
TracePlot(fit = fitkm, par = “beta”)
TracePlot(fit = fitkm, par = “sigsq.eps”)
TracePlot(fit = fitkm, par = “r”, comp = 12)
#提取参数后验变量的重要性指数。
ExtractPIPs(fitkm)
#计算单变量预测结果。
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm,q.fixed=0.5)
#绘图展示，包括散点图以及置信区间，不要看这段代码长就害怕，这段主要是用+连接的绘图函数，相当多的内容是在设置图表参数。
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96se,
ymax = est + 1.96se)) +
geom_hline(yintercept = 0, lty = 2, col = “brown”)+
geom_smooth(stat = “identity”) +
facet_wrap(~variable, ncol = 4) +
xlab(“Urinary metals (Ln, ug/g creatinine)”) +
ylab(“Estimated risk in OA”)+
theme(plot.title = element_text(hjust = 0.5,size = 12, family=“serif”),axis.text=element_text(size=12,family=“serif”),axis.title.x=element_text(size=12,family=“serif”),axis.title.y=element_text(size=12,family=“serif”),strip.text=element_text(size=12,color=“black”, family=“serif”))+
theme(legend.title=element_text(size=12,family=“serif”))
#保存图片到指定路径。
ggsave(filename=“C:/Users/Administrator/Desktop/b1-1.tiff”,plot=plot_1,width =5, height = 6)
#加载office包以使用office软件。
library(eoffice)
#将图形插入PPT。
graph2ppt(file=“C:/Users/Administrator/Desktop/b1-1.tiff”)
#计算整体风险的综合评估并显示结果。
risks.overall <- OverallRiskSummaries(fit = fitkm, qs = seq(0.1, 0.9, by = 0.05), q.fixed = 0.5)
risks.overall
#再次绘图。
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96sd,
ymax = est + 1.96sd)) +
coord_cartesian(ylim = c(-0.3,0.2),xlim = c(0.1,0.9))+
geom_hline(yintercept = 0, lty = 2, col = “brown”) +
geom_pointrange()+
xlab(“Metals (Ln, ug/g creatinine)”) +
ylab(“Estimated OA risk”)+
theme(plot.title = element_text(hjust = 0.5,size = 12, family=“serif”),axis.text=element_text(size=12,family=“serif”),axis.title.x=element_text(size=12,family=“serif”),axis.title.y=element_text(size=12,family=“serif”),strip.text=element_text(size=12,color=“black”, family=“serif”))+
theme(legend.title=element_text(size=12,family=“serif”))
结果解读
图片
可见图A是各个金属暴露物的单独效应，红色虚线为参照基准，可以得知Cd和Co是占比较重的金属暴露物，与WQS的柱状图可以对应上。图B是所有金属拟合当成一种暴露的图，横坐标是百分位数，纵坐标是OA发生风险，红色虚线是基准参照，可知随着金属暴露百分位数的上升，OA的发生风险也上升。
