借助NHANES数据库发文的热度依然不减，对SCI有想法但手里头又没有合适的数据的老师或者同学，真的可以考虑好好筹备一番，抓住“福利的尾巴”。
从NHANES数据库中挖掘特异性的指标，然后将其与疾病相关联，是常用的思路。最近几年，“Mixture analysis”（混合物分析）发展壮大，随之而来的是几种统计分析方法，现在已经发展为大热或必备方法（或者，该称之为惯用套路），例如：
加权分位数和 （Weighted quantile Sum，WQS）回归；
贝叶斯核机回归 （Bayesian Kernel Machine Regression，BKMR）。
前段时间，我们对此做了简单的介绍，还不太了解的朋友可以先看看前面发布的文章：
接下来，笔者将以“gWQS”工具包作者提供的示例作为参照，为大家分享具体分析步骤。
# 查看内置的示例数据，yLBX是连续型的因变量,ybinLBX是二分类变量str(wqs_data)help(wqs_data)# 提取混合变量(要研究的自变量)PCBs = names(wqs_data)[1:34]
先做个探索性建模，看看混合变量对结局变量的影响是正向还是反向。
# 建模# yLBX是因变量，pwqs表示正向效应，nwqs表示负向，# mix_name是混合变量名# q指定分位数，validation指定验证集占比# b指定bootstrap次数（建议增加到最少100）# rh是重复次数（建议增加到最少100）# family指定分布类型（此处为高斯即正态分布）res1 = gwqs(yLBX ~ pwqs + nwqs, mix_name = PCBs, data = wqs_data,q = 10, validation = 0.6, b = 5, b1_pos = TRUE, rh = 2,family = "gaussian", seed = 123)summary(res1)
代码的解释在上述方框内部，结果如下：
因为此建模函数底层是glm，所以直接summary就可以出来结果。直接看pwqs和nwqs的置信区间，可见pwqs显著，nwqs不显著。
可视化结果：
# bar plotgwqs_barplot(res1)# scatter plot y vs wqsgwqs_scatterplot(res1)# scatter plot residuals vs fitted valuesgwqs_fitted_vs_resid(res1)# boxplot of the weights estimated at each repeated holdout stepgwqs_boxplot(results2i)
# 获取估计权重的确切值head(res1$final_weights)# WQS回归的权重表gwqs_weights_tab(res1)
由于负向效应不显著，所以可以直接拟合下面的模型。
res2 = gwqs(yLBX ~ wqs, mix_name = PCBs, data = wqs_data,q = 10, validation = 0.6, b = 5, b1_pos = TRUE, rh = 2,family = "gaussian", seed = 123)summary(res2)
res3 = gwqs(yLBX ~ wqs + sex, mix_name = PCBs, data = wqs_data,q = 10, validation = 0.6, b = 5, b1_pos = TRUE, rh = 2,family = "gaussian", seed = 123)summary(res3)
控制性别后，wqs还是显著的。绘图及其他步骤与上述的一样。
res4 = gwqs(ybinLBX ~ wqs, mix_name = PCBs, data = wqs_data,q = 10, validation = 0.6, b = 5, b1_pos = TRUE, rh = 2,family = "binomial", seed = 123)summary(res4)
gwqs_ROC(res4, wqs_data)参考文献：
[1]R包'gWQS'的帮助文档。
[2]https://renzetti.r-universe.dev/gWQS/doc/manual.html
封面和正文图片来自R软件输出结果，本文仅供学习、分享使用，如有侵权，请联系我们删除，谢谢。
