NHANES 采用了分层四阶段抽样设计,因此在统计分析时,多采用 R 软件 survey 包,进行复杂抽样加权分析。
有不少同学在进行多因素回归分析时,当自变量较多时,R 软件给出的回归结果,所有 P 值均为 NaN,举例如下:
library(survey)
load("./data_frame.RData")
data_frame <- data_to_save
# Constructing a complex survey design object
svy_design <- survey::svydesign(
strata = ~SDMVSTRA,
id = ~SDMVPSU,
weights = ~WTINT2YR,
nest = TRUE,
data = data_frame
)
# Fitting a survey-weighted logistic regression model
mod <- survey::svyglm(
RIAGENDR ~ DMDHRGND + DMDCITZN + SIALANG +
DMDHRAGZ + RIDRETH3 + DMDHRMAZ + INDFMPIR,
family = binomial(),
design = svy_design
运行结果如下:
软件提示 Zero or nagative residual df; p-values not defined
到这里,很多同学懵圈了,不知道下一步怎么办了。
这里我们查阅了survey包的说明书,也和survey包的作者在github里做了交流,我们避开晦涩的统计术语,用临床医生听得懂的话给大家解释一下:
发生这个情况的原因是,模型残差自由度消耗为零。通常由于自变量个数比较多,但样本量又不够大的情况。如何处理呢?
首先是减少自变量的个数。但很多同学并不乐意这么做。那就采用其他理论来设置模型残差的自由度。具体原理如下:
可以在summary(model) 的时候,强行设置模型残差的自由度 df.resid。
1)当不指定 df.resid 时,软件使用 Korn 和 Graubard 推荐的方法计算空模型的自由度(df),通过减法计算残差自由度。这对于 PSU 水平的协变量是正确的,但对于个体水平的协变量可能非常保守。
结果如上图
2)使用 df.resid = degf(result$survey.design) 可计算 PSUs 减去 strata 的数目来获得残差自由度。
结果如下:
3)使用 df.resid = Inf 可基于正态分布进行检验
结果如下:
这样,P 值就计算出来了。
官方对上述方法的评价原文:
If df.resid is not specified the df for the null model is computed by degf and the residual df computed by subtraction. This is recommended by Korn and Graubard and is correct for PSU-level covariates but is potentially very conservative for individual-level covariates. To get tests based on a Normal distribution use df.resid=Inf, and to use number of PSUs-number of strata, specify df.resid=degf(design).
如果觉得以上表述太复杂,也可以用一个傻瓜软件 Mstata 来一键解决问题:
软件位于 www.mstata.com 进入软件后选 “复杂抽样加权(NHANES等)” 下拉菜单,选 “单因素多因素分析” 模块。
我们看到该软件在做复杂抽样加权回归的时候,有个这样的选项:
在选第一个默认设置的时候,有时候会出现残差自由度为零,多因素分析没有 P 值的情况:
软件自动给出了红色报警文字提示,一键修改设置后:
再点分析按钮,则顺利显示多因素分析 P 值:
真是方便得不要不要的。
同时,软件也内置了上述过程的源代码在附件里,需要留底给审稿人的,直接复制软件自动生成的 R 源代码保存即可。