借助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 plot
gwqs_barplot(res1)
# scatter plot y vs wqs
gwqs_scatterplot(res1)
# scatter plot residuals vs fitted values
gwqs_fitted_vs_resid(res1)
# boxplot of the weights estimated at each repeated holdout step
gwqs_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软件输出结果,本文仅供学习、分享使用,如有侵权,请联系我们删除,谢谢。