NHANES挖掘培训班来啦，就在4.20-21！
郑老师团队2024年NHANES公共数据库挖掘培训班，由浅入深，零基础可学，欢迎报名！
近年来，使用美国营养健康（NHANES）数据的文章中，有一类统计学方法异军突起，我称之为回归三板斧，即在统计学设计上同时建立广义线性回归，加权位数和回归以及贝叶斯核机回归三种模型，对比结果比较优劣，再进行综合的分析讨论，得出较为严谨详实的结果。
回归三板斧 | 三种回归模型探讨化学物暴露与肥胖的关联（抓紧上车！）开启了回归三板斧的介绍，这期将完成填坑。
广义线性回归及加权分位数回归可点击下方链接查看：
回归三板斧 | NHANES文章惯用第一招：广义线性回归的前世今生
回归三板斧 | NHANES文章惯用第三招：贝叶斯核机回归（附全套R语言代码）
本次结合上周文章介绍回归三板斧最后一招：贝叶斯核机回归。贝叶斯核机回归在本文中以敏感性分析的形式出现，在推文最后小感悟部分将逐行解析文章关于贝叶斯核机回归的R语言代码。
2022年6月，一篇题为：Biological aging mediates the associations between urinary metals and osteoarthritis among U.S. adults的研究论文发表于《BMC Med》，本文为中国学者写作，文章属于中科院分区医学一区，2023年IF=9.3。
这项研究利用美国营养健康（NHANES）的数据，通过多种方法，研究了金属化学物暴露与骨关节炎之间的关系，以及生物衰老在其中的介导作用。结果表明，金属暴露会增加骨关节炎风险，这可能是由生物衰老介导的。
摘要与主要结果
一、摘要
背景：骨关节炎（OA）是一个全球公共卫生问题，主要发生在老年人中。尽管OA的病因尚不清楚，但环境因素被越来越认为是不可忽视的风险因素。本研究旨在评估尿液中金属元素与骨关节炎风险的关联以及生物衰老的中介效应。
方法：基于美国国家健康与营养调查（NHANES），对12584名美国成年人进行了九种尿液金属浓度检测，包括钡（Ba）、镉（Cd）、钴（Co）、铯（Cs）、钼（Mo）、铅（Pb）、锑（Sb）、铊（Tl）和铀（Tu）。分别采用多变量 logistic 回归和加权分位数和（WQS）回归来探讨单一金属和混合金属与骨关节炎风险之间的关联。此外，从不同角度测量了生物老化，包括细胞衰老（端粒长度）和整体衰老（表型年龄和生物年龄）。进行中介分析以研究生物老化对金属与骨关节炎风险关联的中介效应。
结果：在单一暴露模型中，镉（Cd）、钴（Co）和铯（Cs）与骨关节炎风险呈正相关，其比值比（OR）范围为1.48至1.64（所有P值<0.05）。混合暴露分析显示了一致的关联（OR 1.23，95%CI 1.10至1.37），并强调了镉、钴和铯对结果的贡献。此外，镉、钴、铯、铅（Pb）和铊（Tl）与生物老化标志物呈正相关，而所有生物老化标志物与骨关节炎风险都有显著关联。进一步的中介分析显示，单一金属（主要是镉和铯）和混合金属与骨关节炎风险之间的关联是通过上述生物老化标志物进行中介的，中介比例在16.89%至69.39%之间（所有P值<0.05）。此外，这些关联还通过端粒长度-生物年龄路径和端粒长度-表型年龄路径进行了串行中介（中介比例为4.17%-11.67%），表明金属加速细胞衰老导致整体衰老，最终加剧骨关节炎的进展。
结论：这些研究结果表明，暴露于金属元素会增加骨关节炎的风险，而这种关联可能部分地通过生物老化进行中介作用。
二、研究结果
1. 研究人群的基线资料与金属分布
在12584名成年人中，有1356人被诊断为骨关节炎。表1列出了具有或没有骨关节炎的研究参与者的人口特征。总体而言，年龄、性别、种族、婚姻状况、体力活动、饮酒状况、家庭收入与贫困比率、体重指数、血清烟碱、生物年龄、表型年龄和端粒长度在骨关节炎和非骨关节炎的参与者之间存在显著性差异。
金属浓度的分布情况详见附加文件1：表S1。金属的检测率大于75%。金属的Ln转化后的Pearson系数显示，铯和铊之间存在中等相关性（r = 0.58），钡和钴之间存在相关性（r = 0.41），镉和铅之间存在相关性（r = 0.40），而其他相关性相对较差。
2.金属浓度与骨关节炎风险的关联
图2展示了经过肌酐调整的金属浓度与骨关节炎风险之间的关联，采用了加权逻辑回归模型。相较于第一分位数，镉（OR 1.64，95%CI 1.20至2.23）、钴（OR 1.59，95%CI 1.20至2.10）和铯（OR 1.48，95%CI 1.13至1.93）的最高曝光分位数增加了骨关节炎风险（所有趋势P值<0.05）。这些关联也在Ln转化后的金属浓度与骨关节炎风险之间得到证实（所有P值<0.05）。
混合金属的WQS指数与骨关节炎风险呈正相关（OR 1.23，95%CI 1.10至1.37）（图2）。此外，在WQS模型中，最高权重的金属分别为镉（54.45%）、钴（27.14%）和铯（9.23%）（附加文件1：图S3）。
敏感性分析中，BKMR模型也显示了混合金属与骨关节炎风险的显著正相关（附加文件1：图S4）。在进一步调整参与者职业、糖尿病、高血压、心血管疾病、癌症和骨关节炎相关药物使用、调整调查周期，排除尿肌酐异常或排除孕妇参与者等敏感性分析中，结果没有实质性改变。
3.金属浓度与生物衰老标志物之间的关联
图3展示了基于线性回归的金属与生物老化标志物之间的关联。我们发现镉（Cd）、钴（Co）、铯（Cs）、铅（Pb）和铊（Tl）的最高四分位数（相对于第一四分位数）与生物年龄增长相关（所有趋势P值<0.001）。随着镉（Cd）、钴（Co）、铯（Cs）和铅（Pb）的分位数增加，表型年龄也增加（所有趋势P值<0.001）。镉（Cd）和铯（Cs）与端粒长度呈负相关（所有趋势P值<0.05）。此外，混合金属与生物年龄呈正相关（β 4.91，95%CI 4.52至5.31），与表型年龄呈正相关（β 5.90，95%CI 5.45至6.34），与端粒长度呈负相关（β -0.04，95% CI -0.05至-0.02）。
4.生物老化标志物与骨关节炎风险的关联
表2显示了基于逻辑回归的老化标志物与骨关节炎风险之间的关联。每增加1年的生物年龄，骨关节炎风险增加6%（95%CI 1.05至1.07）。同样，每增加1年的表型年龄与骨关节炎风险增加相关（OR 1.04，95%CI 1.04至1.05）。此外，每单位的端粒长度（平均T/S比值）增加，骨关节炎的OR减少74%（95%CI 0.07至0.97），与分位数分析结果一致（Q4对比Q1：OR 0.40，95%CI 0.17至0.90）。
5.中介分析
此外，进行了并行中介分析以评估生物老化对金属与骨关节炎风险关联的潜在中介效应。生物年龄在镉（Cd）、钴（Co）和铯（Cs）与骨关节炎风险关联中具有显著的中介效应，其中介比例分别为69.27%、18.43%和30.50%（所有P值<0.05）。表型年龄对镉（Cd）、钴（Co）和铯（Cs）与骨关节炎风险关联的中介比例分别为61.01%、31.83%和17.00%。此外，端粒长度也在铯（Cs）与骨关节炎风险关联中起到中介作用，中介比例为9.81%（附加文件1：表S7）。此外，生物年龄、表型年龄和端粒长度分别以57.60%、48.30%和9.50%的中介比例并行中介了混合金属与骨关节炎风险关联（所有P值<0.05）（图4）。
由于细胞衰老是全身老化的主要原因，我们进一步通过串联中介模型探讨了金属与骨关节炎风险关联的潜在途径。（附加文件1：表S8）显示了铯（Cs）与骨关节炎风险之间通过端粒长度-生物年龄途径和端粒长度-表型年龄途径的串联中介效应，中介比例分别为4.55%和4.17%。此外，混合金属与骨关节炎风险关联的端粒长度-生物年龄途径和端粒长度-表型年龄途径的串联中介效应分别为0.70%和0.60%。中介比例分别为11.67%和9.84%（图4）。
设计与统计学方法
一、研究设计
P：1999-2016年美国国家健康和营养检查调查（NHANES）的12584名参与者。
I：未进行分组，暴露因素为尿液金属化学物质。
O：结局：骨关节炎。
S：横断面研究。
二、统计方法
1.介绍以及差异性分析，所有分析都使用SAS（版本9.4）或R（版本3.6.3）进行，并考虑了NHANES的复杂抽样设计，附加文件2中提供了分析的代码。使用卡方检验和t检验评估骨关节炎状况下参与者的人口特征。
2.建立多变量逻辑回归模型，金属浓度被Ln转换为近似正态分布（连续变量），或者被分为四个四分位数（Q1、Q2、Q3和Q4）作为分类变量。应用多变量logistic回归估计金属和生物老化标志物与骨关节炎风险的关联的比值比（OR）及其95%置信区间（CI）。使用多变量线性回归探索金属与生物老化标志物之间的关联。在所有分析中，校正了年龄（<60岁或≥60岁）、性别（男性或女性）、种族/民族（墨西哥美国人、其他西班牙裔、非西班牙裔白人、非西班牙裔黑人或其他人）、婚姻状况（已婚/同居、丧偶/离婚/分居或从未结婚）、体力活动（中等或剧烈）、饮酒状况（曾经或从未）、体重指数（BMI）、家庭收入与贫困之比（PIR）和血清尼古丁浓度的混杂变量。对于混杂变量的缺失数据，对于分类变量，将其编码为缺失指示类别；对于连续变量，用中位数进行填补。在回归模型中使用了虚假发现率（FDR）校正来调整多重检验。使用整数值（1、2、3和4）计算不断增加的暴露组之间的趋势检验。
3.Pearson相关分析，使用Pearson相关分析评估Ln转换后的金属之间的相关性。有向无环图（使用R包“dagitty”和“ggdag”）显示了暴露、结果、混杂变量和中介变量之间的主要关系（附加文件1：图S1）。
4.加权分位数和回归，应用加权分位数和（Weighted Quantile Sum，WQS）回归来探索金属对骨关节炎的总体效应，因为它在描述环境混合物方面表现良好。R包（“gWQS”）可以根据个体金属浓度的加权和来计算WQS指数。WQS指数（范围从0到1）表示金属的混合暴露水平，重要成分通过非零权重来确定。最终结果被解释为混合金属一分位数增加对骨关节炎的同时效应。
5.敏感性分析，还进行了几项敏感性分析。首先，鉴于金属之间存在潜在的非线性和非加性关系，使用贝叶斯核机器回归（Bayesian Kernel Machine Regression，BKMR）评估了所有金属的联合效应（使用"BKMR"包），以及在固定其他金属浓度时单一金属与骨关节炎风险之间的剂量-反应关系。其次，我们进一步调整了参与者的职业（工业和农业、服务和交通等）。第三，我们进一步调整了与年龄相关的疾病，包括糖尿病（是或否）、高血压（是或否）、心血管疾病（是或否）和癌症（是或否），以及与骨关节炎相关的药物使用（如美洛昔康、塞来昔布）（是或否）。第四，我们进一步调整了调查周期（1、2、3、4、5、6、7、8和9）。第五，排除了尿肌酐异常的参与者。最后，孕妇参与者被排除在外。
6.中介分析，通过两个中介模型（平行中介和串行中介分析），估计了生物老化标志物对单一金属和混合金属（以WQS指数表示）与骨关节炎风险的关联的潜在中介效应。平行中介模型使用个体指标作为中介变量，而串行中介模型使用路径作为中介变量。中介分析采用了基于正态近似的千次模拟的准贝叶斯蒙特卡罗方法。直接效应（DE）表示金属暴露对骨关节炎的影响，没有中介变量。间接效应（IE）表示金属暴露通过中介变量对骨关节炎的影响。中介比例通过将IE除以TE（总效应）来计算。
小感悟
本期是回归三板斧系列填坑最终作，结合文章解析贝叶斯核机回归的R语言代码。
BKMR (R)
#首先清理缓存。
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.96*se,
ymax = est + 1.96*se)) +
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.96*sd,
ymax = est + 1.96*sd)) +
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的发生风险也上升。
一个专门做公共数据库的公众号，关注我们
696
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
