上个月有个做营养流行病学的读者私信我说拿NHANES跑出来的BMI均值和CDC官方报告差了十万八千里代码检查来检查去没发现问题最后把权重往survey design里一塞数字才正常。这种事我见过太多次了甚至我自己刚接触NHANES时也栽过同样的跟头。NHANESNational Health and Nutrition Examination Survey是美国国家健康与营养调查一个公开的、免费的、数据质量极高的医学数据库但它的抽样设计是复杂的多阶段分层整群抽样不是我们从统计学课本上学的简单随机抽样。也就是说每个样本并不代表“一个人”而是代表一大群人权重就是用来还原这层代表性的关键。这篇就当保姆级教程我用R语言带你从数据下载、加权描述、加权回归一直做到可视化全文代码都能直接跑跑完你就掌握了NHANES分析最核心的一套流程。适合医学研究生、流行病学研究者以及所有想用公开数据库发文章又卡在加权分析上的朋友。1. NHANES数据库和普通数据集的区别在哪1.1 复杂抽样设计为什么不能把NHANES当普通数据算很多人在学校学统计时默认数据来自简单随机抽样也就是每个人的入选概率都一样所以算均值、算回归都不需要考虑权重。但NHANES不是这么玩的。它为了保证少数族裔、老年人、低收入人群等子群体也有足够的样本量会有意地“过采样”这些人。打个比方如果老年人只占总人口的15%但在NHANES样本里占了25%那你直接对样本人群求平均老年人多的数据点就会把整体均值带偏。加权就是给每个样本一个“放大系数”让样本结构重新匹配美国总人口结构。具体到NHANES它的抽样设计至少包含三个层次分层strata、初级抽样单元PSU、个体individual。所谓分层就是把全国按地理位置、人口密度等分成很多层每层内再抽取若干个县/社区作为PSU然后在PSU内部再抽取居民。这个设计在统计分析里的后果是样本点之间不是独立的同一个PSU里的人往往更相似误差方差的计算方式也和简单随机样本完全不同。所以如果你直接对NHANES数据跑t.test()、lm()、glm()结果肯定有偏。偏多少取决于你关注的变量和过采样的程度。很多时候可不只是几个百分点的差别我见过有人把糖尿病患病率算低了3到4个百分点结论完全反过来。这就是NHANES和普通数据集最大的不同抽样权重和抽样结构必须进入到每一行分析代码里。1.2 权重不加上去偏差到底有多大有人觉得加权只是让均值稍微变一变无伤大雅。真不是。给大家一个直观场景NHANES过采样了老年人而老年人本身更可能患高血压、高血脂、肥胖等问题。如果你不加权就等于把一群“特别容易生病人群”的样本当成了全国人口的缩影患病率自然偏高某些种族亚组的研究里偏差会更明显。再举个我实际跑过的例子用2015-2016年NHANES数据估计20岁以上成人的BMI均值。不加权重时算出来可能接近28.1加权重后可能在27.6到27.8之间。看起来只差0.3、0.4但在大样本下这种差值已经是高度显著了而且如果你按年龄分层画图差异可能集中在某个年龄段被放大。回归分析里权重的加入还会改变回归系数的标准误进而影响P值和置信区间甚至改变是否显著的结论。到写论文时审稿人一看“未加权”大概率直接拒稿因为NHANES官网和绝大多数高质量期刊都要求分析时考虑权重和复杂抽样设计。所以别嫌麻烦加权不是可选项是必选项。2. 动手前必须搞懂的3个基础概念2.1 MEC权重和访谈权重到底用哪个NHANES数据里最常见的权重变量是WTINT2YR和WTMEC2YR。这两个都是两年周期使用的权重但适用场景不同。WTINT2YR对应的是“访谈权重”用于你在家通过问卷获得的信息比如收入、吸烟史、饮食频率这类问卷变量WTMEC2YR对应的是“MEC权重”MEC是Mobile Examination Center的缩写也就是流动体检中心用于你在体检中心测得的体格检查数据、实验室生化指标等。举个例子body mass indexBMXBMI是体检时测出来的所以分析BMI时用WTMEC2YR而你问一个人“每周锻炼几次”这种自报数据则用WTINT2YR。如果你同时分析问卷变量和体检变量通常也是以MEC权重为准因为最终纳入的模型里包含了体检结果。最简单的记忆方式是凡是要进模型的变量里有体格检查或实验室指标就优先用WTMEC2YR纯问卷数据则用WTINT2YR。另外不同周期还有对应的权重版本。比如2015-2016年单周期数据里有WTMEC2YR如果你把两个连续周期合并成四年数据使用就要用WTMEC4YR如果自己手动合并两个周期还需要对权重做重新计算公式我会在后续常见问题部分详细说。总之先记住一个原则你用了哪些周期的数据就要匹配对应周期的权重变量。2.2 SDMVPSU和SDMVSTRA分层与聚类的身份信息除了权重NHANES数据里还有两个特别关键的变量SDMVPSU和SDMVSTRA。SDMVPSU是初级抽样单元的编号SDMVSTRA是分层的编号。这两个变量合在一起唯一标识了每个样本在抽样设计里的“位置”。在R的survey包里构建设计对象时需要用id ~SDMVPSU指定PSU用strata ~SDMVSTRA指定分层并且设置nest TRUE。为什么要指定它们因为标准误的计算必须考虑“层内抽样”的方差结构。你可以把PSU理解成一个“班级”层理解成“学校”同一班级里的学生成绩往往比不同班级的成绩更接近。如果你忽视这个结构计算标准误时等于是把所有个体当成独立样本这样算出来的标准误会偏小置信区间就会太窄结论显得过于“精确”这在统计分析里是很危险的事情。这也是我看到很多拿NHANES做分析的人犯得最隐蔽的一个错误权重加了但survey design里的id和strata没写对导致方差估计仍然不对。具体到R代码标准写法是svydesign(id ~SDMVPSU, strata ~SDMVSTRA, weights ~WTMEC2YR, nest TRUE, data xxx)。关于nest TRUE的含义简单说就是告诉R“PSU编号在各层之间是重复出现的”这是NHANES数据的标准设置。这些细节不需要理解得很深你只要照做统计结果就和主流论文对齐了。2.3 数据下载一条R命令把三张表拉下来以前下载NHANES数据挺痛苦的需要在官网一个个点文件、下载XPT、再读入R。现在有了nhanesA这个包事情简单多了。安装并加载nhanesA后用nhanes(表名)就能拉取指定表。比较常用的表名有DEMO人口学信息、BMX体格检查、BPX血压、DIQ糖尿病问卷、HDL血脂等等。我这次演示拉三张表DEMO、BMX、BPX分别代表人口学信息、身体测量指标和血压测量结果。这三张表基本够我们跑一个完整的“BMI和高血压关系”的演示。代码很简单核心就一行一个表。需要提醒的是nhanesA的接口在不同版本上略有差异有的版本要带year 2015-2016参数有的版本默认拉最新周期。如果你不想纠结版本问题也可以直接从官网下载XPT文件再用haven::read_xpt()读取效果完全一样。下载部分我把两种方式都写在代码里你根据自己的网络和包版本选一种就行。3. 5分钟核心流程完整可跑代码3.1 环境准备与数据加载推荐先新建一个RStudio项目把代码写在脚本文件里不要直接在控制台里临时敲。下面这段代码完成了安装、加载、数据下载和合并标注清楚每一步在干什么。# 1. 安装包如果已经装过可以跳过 install.packages(c(tidyverse, survey, srvyr, nhanesA, haven)) # 2. 加载包 library(tidyverse) library(survey) library(srvyr) library(nhanesA) # 3. 下载NHANES 2015-2016周期的三张表 # 如果版本支持 year 参数就带 year不支持就直接 nhanes(DEMO) demo - nhanes(DEMO, year 2015-2016) bmx - nhanes(BMX, year 2015-2016) bpx - nhanes(BPX, year 2015-2016) # 4. 按SEQN样本唯一ID合并 nhanes_data - demo %% left_join(bmx, by SEQN) %% left_join(bpx, by SEQN)如果nhanesA在你机器上拉数据太慢或者一直报错我提供一个备选方案直接从官网把XPT文件下到本地再读入。这里以2015-2016年为例文件名中的I代表该周期代码# 备选方案直接下载XPT文件 download.file( https://wwwn.cdc.gov/Nchs/Nhanes/2015-2016/DEMO_I.XPT, DEMO_I.XPT, mode wb ) demo - haven::read_xpt(DEMO_I.XPT) # BMX_I.XPT 和 BPX_I.XPT 同理读完再按SEQN合并第一次下载数据确实会花几分钟但这是值得的。我的习惯是把合并后的数据保存成本地RDS文件之后所有分析直接readRDS()读取就不用再等了。代码是saveRDS(nhanes_data, nhanes_2015_2016.rds)下次load进来就行。3.2 构建survey.design让R知道数据是“加权”的数据合并好之后下一步就是构建survey design对象。这一步是整个NHANES分析的核心入口权重、分层、聚类都在这一个对象里声明。代码如下# 5. 关键设置处理稀有层/独立PSU的方差估计问题 options(survey.lonely.psu adjust) # 6. 构建抽样设计 nhanes_design - svydesign( id ~SDMVPSU, # 初级抽样单元 strata ~SDMVSTRA, # 分层变量 weights ~WTMEC2YR, # MEC权重 nest TRUE, data nhanes_data )为什么要加options(survey.lonely.psu adjust)这一句因为NHANES在某些层里可能只有一个PSU这种“孤独PSU”会导致方差估计出问题。survey包遇到这种情况时默认会报错或者给一个奇怪的标准误。设置成adjust后它会用中心化调整的思路处理这在实际分析里是比较常见的做法很多官方教程也推荐。如果不设置你跑某些亚组分析时可能会看到warning甚至直接报错所以强烈建议在任何svydesign构建之前就把这个全局选项设好。另外要注意svydesign里的data参数必须是data.frame而且weights参数要写成~WTMEC2YR这种公式形式不是直接写列名。这是survey包的固定语法很多刚上手的人会在这一步卡住。3.3 加权描述统计均值、比例、分组对比构建好survey design之后描述统计就很简单了。先看加权均值是多少# 7. 加权BMI均值 svymean(~BMXBMI, nhanes_design, na.rm TRUE)svymean返回的不仅有点估计还有标准误这个标准误已经包含了抽样设计的方差结构。如果你想看加权比例比如高血压的患病率先把血压变量转换成二分类再用svymean()# 8. 生成“是否高血压”变量 nhanes_design - update(nhanes_design, hypertension ifelse(BPXSAR 130 | BPXSAR 80, 1, 0) ) # 加权高血压患病率 svymean(~hypertension, nhanes_design, na.rm TRUE)注意这里用了一个不太显眼但很重要的函数update()它是survey包专用的修改变量方式能够让新变量被纳入已经定义好的design对象里。如果你直接往nhanes_data这个原始data.frame里加列design对象并不知道你加了什么进而可能报错。分组描述则用svyby()。假设想看不同性别的BMI均值# 9. 按性别分组的加权BMI均值 svyby(~BMXBMI, ~RIAGENDR, nhanes_design, svymean, na.rm TRUE)RIAGENDR是NHANES里的性别变量1代表男性2代表女性。svyby的输出会同时给出每组的均值和标准误。这里有一个非常推荐的对比技巧先跑一遍不加权版本再跑一遍加权版本专门看看差异有多大。不加权版本就是mean(nhanes_data$BMXBMI, na.rm TRUE)或者tapply一个循环。如果两个版本差异很小说明变量在权重上的敏感性不高但写论文时仍然要报告加权结果如果差异很大恭喜你这个发现本身就可以写进讨论部分说明样本结构对结果的影响显著。3.4 加权回归svyglm做线性和逻辑回归描述统计只是开胃菜回归分析才是NHANES产出结论的重头戏。survey包里的回归函数叫svyglm()用法和glm()非常接近但会自动把抽样设计考虑进去。线性回归示例研究年龄、性别和BMI的关系。# 10. 加权线性回归 model_lm - svyglm(BMXBMI ~ RIDAGEYR factor(RIAGENDR), design nhanes_design, na.rm TRUE) summary(model_lm)RIDAGEYR是年龄连续变量factor(RIAGENDR)把性别转为分类变量。summary的输出里每一个回归系数都带上了加权标准误。如果你是做流行病学研究的通常还想看二分类结局比如高血压是否受年龄和性别影响那就用逻辑回归。这块有一个特别臭名昭著的坑在svyglm中做logistic回归family参数不要写binomial而要写quasibinomial。# 11. 加权逻辑回归 model_glm - svyglm(hypertension ~ RIDAGEYR factor(RIAGENDR), design nhanes_design, family quasibinomial(), na.rm TRUE) summary(model_glm)为什么用quasibinomial因为受复杂抽样设计影响观测不一定完全满足二项分布的独立性假设survey包用准二项分布来处理过离散的问题实际估计的系数和binomial几乎一致但标准误会更稳健。如果你用family binomial()某些版本下会弹出warning甚至直接报错说“non-integer #successes”原因和迭代加权最小二乘的算法有关。新手最容易在这个坎上停住直接把binomial改成quasibinomial就好。拿到系数之后想报告OR值和95%置信区间也很简单# 12. 输出OR和95%CI exp(cbind(OR coef(model_glm), confint(model_glm)))这里把logist回归的log-odds系数指数化就得到了OR值。需要注意NA存在时svyglm对象本身可能已经通过na.action处理掉了缺失样本confint()能自动给出基于校正自由度的置信区间。看到区间不包括1就说明该变量在0.05水平上显著。3.5 可视化把权重画进图里很多人在可视化阶段犯错用加权结果构建了模型画图时直接拿原始数据丢给ggplot2。这样画出来的描述图柱状图和箱线图全都是未加权的跟前面表格的加权结果对不上审稿人一眼就能看出来。正确的思路是先算加权统计量再画图。这里我用srvyr包配合ggplot2来完成步骤很顺。# 13. 转成srvyr设计对象语法更接近dplyr nhanes_svy - nhanes_data %% as_survey_design(ids SDMVPSU, strata SDMVSTRA, weights WTMEC2YR, nest TRUE) # 14. 按年龄分组计算加权BMI均值和标准误 bmi_by_age - nhanes_svy %% filter(!is.na(BMXBMI)) %% mutate(age_grp cut(RIDAGEYR, breaks c(20, 40, 60, 80), include.lowest TRUE)) %% group_by(age_grp) %% summarise(mean_bmi survey_mean(BMXBMI, na.rm TRUE)) # 15. 用ggplot2画柱状图误差线 library(ggplot2) ggplot(bmi_by_age, aes(x age_grp, y mean_bmi)) geom_col(fill #4E79A7) geom_errorbar(aes(ymin mean_bmi - 1.96 * mean_bmi_se, ymax mean_bmi 1.96 * mean_bmi_se), width 0.2) labs(x 年龄组, y 加权BMI均值) theme_minimal(base_size 14)srvyr包的神奇之处在于它的函数名和dplyr几乎一样group_by之后用survey_mean()就是在分组计算加权均值。结果里的mean_bmi_se就是标准误画误差线时用1.96乘上它得到近似的95%置信区间。如果你觉得柱状图太简单也可以画分组箱线图但要注意ggplot2的geom_boxplot不认权重这时最好改用加权分位数比如用svyquantile算出第25、50、75百分位数再画。否则你看到的箱线图又是未加权的分布属于经典错误。另外一个好用的可视化是加权直方图或密度图用来展示BMI的加权分布。这在survey包里没有特别漂亮的包装函数但原理很简单画图前把每条记录的权重当作geom_histogram的weight参数。# 16. 加权BMI直方图 ggplot(nhanes_svy$variables, aes(x BMXBMI)) geom_histogram(aes(weight WTMEC2YR), bins 40, fill #4E79A7, color white) labs(x BMI, y 加权频数) theme_minimal(base_size 14)nhanes_svy$variables就是原始data.frame直接取出来画图加上weight参数后直方图反映的就是加权后的分布。注意这里没有考虑复杂抽样设计的方差只是描述性展示但比未加权版本更有意义。4. 常见问题与排查技巧实录4.1 多个周期合并时权重怎么重新计算新手最常见的需求是合并两个连续周期比如2015-2016和2017-2018这样可以提高亚组分析的检验效能。合并之后权重不能直接用WTMEC2YR。正确做法是如果使用的是连续两个周期且二者都是2年数据新权重 WTMEC2YR / 2。原因很简单个体代表的人口总量翻倍了每个样本的相对权重减半这样才能保证合并后总人数接近真实人口规模。如果是2017-2020这种特殊性周期因为本身覆盖了多个年份官方已经给出了专门的4年权重变量。你可以查一下对应周期的文档说明确认它自带的就是WTMEC4YR。凡是合并周期的数据构建survey design时weights列一定要写调整后的新权重列而不是原始的WTMEC2YR。这一步错掉后面的所有统计量都会失去代表性。4.2 报错lonely PSU / 自由度太小运行svyby或svyglm时如果出现“lonely PSU”警告说明某些层里只有一个PSU。这个问题的经典解法就是前面提到的options(survey.lonely.psu adjust)一定要写在svydesign之前。如果你在自定义函数里跑survey分析也要在函数里再次设置这个选项因为R启动后默认可能被重置。如果你发现置信区间特别宽甚至可以宽到没有意义那多半是自由度太小。NHANES每个周期的设计自由度大概是PSU数减去层数通常在几十左右。自由度越小t分布尾部越厚置信区间越宽。这不是代码错误是复杂抽样的内在特性。稳妥的做法是报告结果时主动指出自由度并解释这种保守的估计是为了覆盖抽样设计的不确定性。4.3 亚组分析时权重要不要重新算有人做男性亚组分析时只把男性数据筛选出来再重新构建一个svydesign并继续使用原来的权重这是不对的但又不能简单地把权重加起来归一化。正确的做法是用subset()函数在原design对象上做子集化male_design - subset(nhanes_design, RIAGENDR 1)subset()是survey包专为处理设计对象而生的方法它会保留原始的层和PSU结构同时更新权重使子集内的估计仍然有代表性。如果你用dplyr的filter()先过滤再重新svydesign可能会导致部分层的PSU数量为0进而出现方差估计异常。所以凡是需要亚组分析一律用subset()而不是filter重建设计。这是一个能让R帮你省去大量麻烦的习惯。4.4 缺失值处理顺序和其他细节svymean和svyglm里都看到na.rm TRUE或na.action参数但要注意描述统计里的na.rm TRUE只对当前变量生效回归分析里缺失值是按整行删除的。因此如果你在做回归前发现样本量骤减先检查是不是合并变量时把很多行都变成NA了特别是BPX血压变量可能因为检测失败而缺失。我通常会先构建一个只含模型所需变量的data.frame同时检查缺失数量再决定是直接丢弃还是用多重插补。多重插补在复杂抽样设计下比较麻烦常规论文里一般用complete case analysis就够你只要在方法部分写清楚即可。另外一个容易被忽略的坑是变量合并时的列名冲突。BMX和BPX里都有SEQN这是对的但某些变量名可能在不同表之间重复比如都叫“XXX”left_join时会出现.x和.y后缀后续引用变量名就容易搞混。建议每次合并后用一个简单的names()检查一下或者在下一次赋值前把不需要的列删干净。关于可视化再补一句心里话不要一上来就追求花哨的可视化大屏NHANES这种复杂抽样数据的可视化最重要的是“统计量算对”尤其是误差线必须来自加权标准误。图丑一点没关系画错了才致命。等加权收入、加权患病率这些核心数字都准确了再考虑配色、字体、大屏展示这些锦上添花的东西。个人实操下来最冤的情况就是辛辛苦苦下载数据、清洗合并结果因为少设置一个survey.lonely.psu选项全部分析不能跑。我现在的习惯是每次新建R脚本第一件事就是把options(survey.lonely.psu adjust)和library(survey)写在最前面再开始写业务逻辑。这个微不足道的小习惯能帮你省下大量查报错的时间。NHANES分析并不难难点在于你愿不愿意把抽样设计当成第一公民来看待而一旦你掌握了survey包这套工作流以后处理任何复杂抽样数据库都是同一个套路一劳永逸。
上个月有个做营养流行病学的读者私信我说拿NHANES跑出来的BMI均值和CDC官方报告差了十万八千里代码检查来检查去没发现问题最后把权重往survey design里一塞数字才正常。这种事我见过太多次了甚至我自己刚接触NHANES时也栽过同样的跟头。NHANESNational Health and Nutrition Examination Survey是美国国家健康与营养调查一个公开的、免费的、数据质量极高的医学数据库但它的抽样设计是复杂的多阶段分层整群抽样不是我们从统计学课本上学的简单随机抽样。也就是说每个样本并不代表“一个人”而是代表一大群人权重就是用来还原这层代表性的关键。这篇就当保姆级教程我用R语言带你从数据下载、加权描述、加权回归一直做到可视化全文代码都能直接跑跑完你就掌握了NHANES分析最核心的一套流程。适合医学研究生、流行病学研究者以及所有想用公开数据库发文章又卡在加权分析上的朋友。1. NHANES数据库和普通数据集的区别在哪1.1 复杂抽样设计为什么不能把NHANES当普通数据算很多人在学校学统计时默认数据来自简单随机抽样也就是每个人的入选概率都一样所以算均值、算回归都不需要考虑权重。但NHANES不是这么玩的。它为了保证少数族裔、老年人、低收入人群等子群体也有足够的样本量会有意地“过采样”这些人。打个比方如果老年人只占总人口的15%但在NHANES样本里占了25%那你直接对样本人群求平均老年人多的数据点就会把整体均值带偏。加权就是给每个样本一个“放大系数”让样本结构重新匹配美国总人口结构。具体到NHANES它的抽样设计至少包含三个层次分层strata、初级抽样单元PSU、个体individual。所谓分层就是把全国按地理位置、人口密度等分成很多层每层内再抽取若干个县/社区作为PSU然后在PSU内部再抽取居民。这个设计在统计分析里的后果是样本点之间不是独立的同一个PSU里的人往往更相似误差方差的计算方式也和简单随机样本完全不同。所以如果你直接对NHANES数据跑t.test()、lm()、glm()结果肯定有偏。偏多少取决于你关注的变量和过采样的程度。很多时候可不只是几个百分点的差别我见过有人把糖尿病患病率算低了3到4个百分点结论完全反过来。这就是NHANES和普通数据集最大的不同抽样权重和抽样结构必须进入到每一行分析代码里。1.2 权重不加上去偏差到底有多大有人觉得加权只是让均值稍微变一变无伤大雅。真不是。给大家一个直观场景NHANES过采样了老年人而老年人本身更可能患高血压、高血脂、肥胖等问题。如果你不加权就等于把一群“特别容易生病人群”的样本当成了全国人口的缩影患病率自然偏高某些种族亚组的研究里偏差会更明显。再举个我实际跑过的例子用2015-2016年NHANES数据估计20岁以上成人的BMI均值。不加权重时算出来可能接近28.1加权重后可能在27.6到27.8之间。看起来只差0.3、0.4但在大样本下这种差值已经是高度显著了而且如果你按年龄分层画图差异可能集中在某个年龄段被放大。回归分析里权重的加入还会改变回归系数的标准误进而影响P值和置信区间甚至改变是否显著的结论。到写论文时审稿人一看“未加权”大概率直接拒稿因为NHANES官网和绝大多数高质量期刊都要求分析时考虑权重和复杂抽样设计。所以别嫌麻烦加权不是可选项是必选项。2. 动手前必须搞懂的3个基础概念2.1 MEC权重和访谈权重到底用哪个NHANES数据里最常见的权重变量是WTINT2YR和WTMEC2YR。这两个都是两年周期使用的权重但适用场景不同。WTINT2YR对应的是“访谈权重”用于你在家通过问卷获得的信息比如收入、吸烟史、饮食频率这类问卷变量WTMEC2YR对应的是“MEC权重”MEC是Mobile Examination Center的缩写也就是流动体检中心用于你在体检中心测得的体格检查数据、实验室生化指标等。举个例子body mass indexBMXBMI是体检时测出来的所以分析BMI时用WTMEC2YR而你问一个人“每周锻炼几次”这种自报数据则用WTINT2YR。如果你同时分析问卷变量和体检变量通常也是以MEC权重为准因为最终纳入的模型里包含了体检结果。最简单的记忆方式是凡是要进模型的变量里有体格检查或实验室指标就优先用WTMEC2YR纯问卷数据则用WTINT2YR。另外不同周期还有对应的权重版本。比如2015-2016年单周期数据里有WTMEC2YR如果你把两个连续周期合并成四年数据使用就要用WTMEC4YR如果自己手动合并两个周期还需要对权重做重新计算公式我会在后续常见问题部分详细说。总之先记住一个原则你用了哪些周期的数据就要匹配对应周期的权重变量。2.2 SDMVPSU和SDMVSTRA分层与聚类的身份信息除了权重NHANES数据里还有两个特别关键的变量SDMVPSU和SDMVSTRA。SDMVPSU是初级抽样单元的编号SDMVSTRA是分层的编号。这两个变量合在一起唯一标识了每个样本在抽样设计里的“位置”。在R的survey包里构建设计对象时需要用id ~SDMVPSU指定PSU用strata ~SDMVSTRA指定分层并且设置nest TRUE。为什么要指定它们因为标准误的计算必须考虑“层内抽样”的方差结构。你可以把PSU理解成一个“班级”层理解成“学校”同一班级里的学生成绩往往比不同班级的成绩更接近。如果你忽视这个结构计算标准误时等于是把所有个体当成独立样本这样算出来的标准误会偏小置信区间就会太窄结论显得过于“精确”这在统计分析里是很危险的事情。这也是我看到很多拿NHANES做分析的人犯得最隐蔽的一个错误权重加了但survey design里的id和strata没写对导致方差估计仍然不对。具体到R代码标准写法是svydesign(id ~SDMVPSU, strata ~SDMVSTRA, weights ~WTMEC2YR, nest TRUE, data xxx)。关于nest TRUE的含义简单说就是告诉R“PSU编号在各层之间是重复出现的”这是NHANES数据的标准设置。这些细节不需要理解得很深你只要照做统计结果就和主流论文对齐了。2.3 数据下载一条R命令把三张表拉下来以前下载NHANES数据挺痛苦的需要在官网一个个点文件、下载XPT、再读入R。现在有了nhanesA这个包事情简单多了。安装并加载nhanesA后用nhanes(表名)就能拉取指定表。比较常用的表名有DEMO人口学信息、BMX体格检查、BPX血压、DIQ糖尿病问卷、HDL血脂等等。我这次演示拉三张表DEMO、BMX、BPX分别代表人口学信息、身体测量指标和血压测量结果。这三张表基本够我们跑一个完整的“BMI和高血压关系”的演示。代码很简单核心就一行一个表。需要提醒的是nhanesA的接口在不同版本上略有差异有的版本要带year 2015-2016参数有的版本默认拉最新周期。如果你不想纠结版本问题也可以直接从官网下载XPT文件再用haven::read_xpt()读取效果完全一样。下载部分我把两种方式都写在代码里你根据自己的网络和包版本选一种就行。3. 5分钟核心流程完整可跑代码3.1 环境准备与数据加载推荐先新建一个RStudio项目把代码写在脚本文件里不要直接在控制台里临时敲。下面这段代码完成了安装、加载、数据下载和合并标注清楚每一步在干什么。# 1. 安装包如果已经装过可以跳过 install.packages(c(tidyverse, survey, srvyr, nhanesA, haven)) # 2. 加载包 library(tidyverse) library(survey) library(srvyr) library(nhanesA) # 3. 下载NHANES 2015-2016周期的三张表 # 如果版本支持 year 参数就带 year不支持就直接 nhanes(DEMO) demo - nhanes(DEMO, year 2015-2016) bmx - nhanes(BMX, year 2015-2016) bpx - nhanes(BPX, year 2015-2016) # 4. 按SEQN样本唯一ID合并 nhanes_data - demo %% left_join(bmx, by SEQN) %% left_join(bpx, by SEQN)如果nhanesA在你机器上拉数据太慢或者一直报错我提供一个备选方案直接从官网把XPT文件下到本地再读入。这里以2015-2016年为例文件名中的I代表该周期代码# 备选方案直接下载XPT文件 download.file( https://wwwn.cdc.gov/Nchs/Nhanes/2015-2016/DEMO_I.XPT, DEMO_I.XPT, mode wb ) demo - haven::read_xpt(DEMO_I.XPT) # BMX_I.XPT 和 BPX_I.XPT 同理读完再按SEQN合并第一次下载数据确实会花几分钟但这是值得的。我的习惯是把合并后的数据保存成本地RDS文件之后所有分析直接readRDS()读取就不用再等了。代码是saveRDS(nhanes_data, nhanes_2015_2016.rds)下次load进来就行。3.2 构建survey.design让R知道数据是“加权”的数据合并好之后下一步就是构建survey design对象。这一步是整个NHANES分析的核心入口权重、分层、聚类都在这一个对象里声明。代码如下# 5. 关键设置处理稀有层/独立PSU的方差估计问题 options(survey.lonely.psu adjust) # 6. 构建抽样设计 nhanes_design - svydesign( id ~SDMVPSU, # 初级抽样单元 strata ~SDMVSTRA, # 分层变量 weights ~WTMEC2YR, # MEC权重 nest TRUE, data nhanes_data )为什么要加options(survey.lonely.psu adjust)这一句因为NHANES在某些层里可能只有一个PSU这种“孤独PSU”会导致方差估计出问题。survey包遇到这种情况时默认会报错或者给一个奇怪的标准误。设置成adjust后它会用中心化调整的思路处理这在实际分析里是比较常见的做法很多官方教程也推荐。如果不设置你跑某些亚组分析时可能会看到warning甚至直接报错所以强烈建议在任何svydesign构建之前就把这个全局选项设好。另外要注意svydesign里的data参数必须是data.frame而且weights参数要写成~WTMEC2YR这种公式形式不是直接写列名。这是survey包的固定语法很多刚上手的人会在这一步卡住。3.3 加权描述统计均值、比例、分组对比构建好survey design之后描述统计就很简单了。先看加权均值是多少# 7. 加权BMI均值 svymean(~BMXBMI, nhanes_design, na.rm TRUE)svymean返回的不仅有点估计还有标准误这个标准误已经包含了抽样设计的方差结构。如果你想看加权比例比如高血压的患病率先把血压变量转换成二分类再用svymean()# 8. 生成“是否高血压”变量 nhanes_design - update(nhanes_design, hypertension ifelse(BPXSAR 130 | BPXSAR 80, 1, 0) ) # 加权高血压患病率 svymean(~hypertension, nhanes_design, na.rm TRUE)注意这里用了一个不太显眼但很重要的函数update()它是survey包专用的修改变量方式能够让新变量被纳入已经定义好的design对象里。如果你直接往nhanes_data这个原始data.frame里加列design对象并不知道你加了什么进而可能报错。分组描述则用svyby()。假设想看不同性别的BMI均值# 9. 按性别分组的加权BMI均值 svyby(~BMXBMI, ~RIAGENDR, nhanes_design, svymean, na.rm TRUE)RIAGENDR是NHANES里的性别变量1代表男性2代表女性。svyby的输出会同时给出每组的均值和标准误。这里有一个非常推荐的对比技巧先跑一遍不加权版本再跑一遍加权版本专门看看差异有多大。不加权版本就是mean(nhanes_data$BMXBMI, na.rm TRUE)或者tapply一个循环。如果两个版本差异很小说明变量在权重上的敏感性不高但写论文时仍然要报告加权结果如果差异很大恭喜你这个发现本身就可以写进讨论部分说明样本结构对结果的影响显著。3.4 加权回归svyglm做线性和逻辑回归描述统计只是开胃菜回归分析才是NHANES产出结论的重头戏。survey包里的回归函数叫svyglm()用法和glm()非常接近但会自动把抽样设计考虑进去。线性回归示例研究年龄、性别和BMI的关系。# 10. 加权线性回归 model_lm - svyglm(BMXBMI ~ RIDAGEYR factor(RIAGENDR), design nhanes_design, na.rm TRUE) summary(model_lm)RIDAGEYR是年龄连续变量factor(RIAGENDR)把性别转为分类变量。summary的输出里每一个回归系数都带上了加权标准误。如果你是做流行病学研究的通常还想看二分类结局比如高血压是否受年龄和性别影响那就用逻辑回归。这块有一个特别臭名昭著的坑在svyglm中做logistic回归family参数不要写binomial而要写quasibinomial。# 11. 加权逻辑回归 model_glm - svyglm(hypertension ~ RIDAGEYR factor(RIAGENDR), design nhanes_design, family quasibinomial(), na.rm TRUE) summary(model_glm)为什么用quasibinomial因为受复杂抽样设计影响观测不一定完全满足二项分布的独立性假设survey包用准二项分布来处理过离散的问题实际估计的系数和binomial几乎一致但标准误会更稳健。如果你用family binomial()某些版本下会弹出warning甚至直接报错说“non-integer #successes”原因和迭代加权最小二乘的算法有关。新手最容易在这个坎上停住直接把binomial改成quasibinomial就好。拿到系数之后想报告OR值和95%置信区间也很简单# 12. 输出OR和95%CI exp(cbind(OR coef(model_glm), confint(model_glm)))这里把logist回归的log-odds系数指数化就得到了OR值。需要注意NA存在时svyglm对象本身可能已经通过na.action处理掉了缺失样本confint()能自动给出基于校正自由度的置信区间。看到区间不包括1就说明该变量在0.05水平上显著。3.5 可视化把权重画进图里很多人在可视化阶段犯错用加权结果构建了模型画图时直接拿原始数据丢给ggplot2。这样画出来的描述图柱状图和箱线图全都是未加权的跟前面表格的加权结果对不上审稿人一眼就能看出来。正确的思路是先算加权统计量再画图。这里我用srvyr包配合ggplot2来完成步骤很顺。# 13. 转成srvyr设计对象语法更接近dplyr nhanes_svy - nhanes_data %% as_survey_design(ids SDMVPSU, strata SDMVSTRA, weights WTMEC2YR, nest TRUE) # 14. 按年龄分组计算加权BMI均值和标准误 bmi_by_age - nhanes_svy %% filter(!is.na(BMXBMI)) %% mutate(age_grp cut(RIDAGEYR, breaks c(20, 40, 60, 80), include.lowest TRUE)) %% group_by(age_grp) %% summarise(mean_bmi survey_mean(BMXBMI, na.rm TRUE)) # 15. 用ggplot2画柱状图误差线 library(ggplot2) ggplot(bmi_by_age, aes(x age_grp, y mean_bmi)) geom_col(fill #4E79A7) geom_errorbar(aes(ymin mean_bmi - 1.96 * mean_bmi_se, ymax mean_bmi 1.96 * mean_bmi_se), width 0.2) labs(x 年龄组, y 加权BMI均值) theme_minimal(base_size 14)srvyr包的神奇之处在于它的函数名和dplyr几乎一样group_by之后用survey_mean()就是在分组计算加权均值。结果里的mean_bmi_se就是标准误画误差线时用1.96乘上它得到近似的95%置信区间。如果你觉得柱状图太简单也可以画分组箱线图但要注意ggplot2的geom_boxplot不认权重这时最好改用加权分位数比如用svyquantile算出第25、50、75百分位数再画。否则你看到的箱线图又是未加权的分布属于经典错误。另外一个好用的可视化是加权直方图或密度图用来展示BMI的加权分布。这在survey包里没有特别漂亮的包装函数但原理很简单画图前把每条记录的权重当作geom_histogram的weight参数。# 16. 加权BMI直方图 ggplot(nhanes_svy$variables, aes(x BMXBMI)) geom_histogram(aes(weight WTMEC2YR), bins 40, fill #4E79A7, color white) labs(x BMI, y 加权频数) theme_minimal(base_size 14)nhanes_svy$variables就是原始data.frame直接取出来画图加上weight参数后直方图反映的就是加权后的分布。注意这里没有考虑复杂抽样设计的方差只是描述性展示但比未加权版本更有意义。4. 常见问题与排查技巧实录4.1 多个周期合并时权重怎么重新计算新手最常见的需求是合并两个连续周期比如2015-2016和2017-2018这样可以提高亚组分析的检验效能。合并之后权重不能直接用WTMEC2YR。正确做法是如果使用的是连续两个周期且二者都是2年数据新权重 WTMEC2YR / 2。原因很简单个体代表的人口总量翻倍了每个样本的相对权重减半这样才能保证合并后总人数接近真实人口规模。如果是2017-2020这种特殊性周期因为本身覆盖了多个年份官方已经给出了专门的4年权重变量。你可以查一下对应周期的文档说明确认它自带的就是WTMEC4YR。凡是合并周期的数据构建survey design时weights列一定要写调整后的新权重列而不是原始的WTMEC2YR。这一步错掉后面的所有统计量都会失去代表性。4.2 报错lonely PSU / 自由度太小运行svyby或svyglm时如果出现“lonely PSU”警告说明某些层里只有一个PSU。这个问题的经典解法就是前面提到的options(survey.lonely.psu adjust)一定要写在svydesign之前。如果你在自定义函数里跑survey分析也要在函数里再次设置这个选项因为R启动后默认可能被重置。如果你发现置信区间特别宽甚至可以宽到没有意义那多半是自由度太小。NHANES每个周期的设计自由度大概是PSU数减去层数通常在几十左右。自由度越小t分布尾部越厚置信区间越宽。这不是代码错误是复杂抽样的内在特性。稳妥的做法是报告结果时主动指出自由度并解释这种保守的估计是为了覆盖抽样设计的不确定性。4.3 亚组分析时权重要不要重新算有人做男性亚组分析时只把男性数据筛选出来再重新构建一个svydesign并继续使用原来的权重这是不对的但又不能简单地把权重加起来归一化。正确的做法是用subset()函数在原design对象上做子集化male_design - subset(nhanes_design, RIAGENDR 1)subset()是survey包专为处理设计对象而生的方法它会保留原始的层和PSU结构同时更新权重使子集内的估计仍然有代表性。如果你用dplyr的filter()先过滤再重新svydesign可能会导致部分层的PSU数量为0进而出现方差估计异常。所以凡是需要亚组分析一律用subset()而不是filter重建设计。这是一个能让R帮你省去大量麻烦的习惯。4.4 缺失值处理顺序和其他细节svymean和svyglm里都看到na.rm TRUE或na.action参数但要注意描述统计里的na.rm TRUE只对当前变量生效回归分析里缺失值是按整行删除的。因此如果你在做回归前发现样本量骤减先检查是不是合并变量时把很多行都变成NA了特别是BPX血压变量可能因为检测失败而缺失。我通常会先构建一个只含模型所需变量的data.frame同时检查缺失数量再决定是直接丢弃还是用多重插补。多重插补在复杂抽样设计下比较麻烦常规论文里一般用complete case analysis就够你只要在方法部分写清楚即可。另外一个容易被忽略的坑是变量合并时的列名冲突。BMX和BPX里都有SEQN这是对的但某些变量名可能在不同表之间重复比如都叫“XXX”left_join时会出现.x和.y后缀后续引用变量名就容易搞混。建议每次合并后用一个简单的names()检查一下或者在下一次赋值前把不需要的列删干净。关于可视化再补一句心里话不要一上来就追求花哨的可视化大屏NHANES这种复杂抽样数据的可视化最重要的是“统计量算对”尤其是误差线必须来自加权标准误。图丑一点没关系画错了才致命。等加权收入、加权患病率这些核心数字都准确了再考虑配色、字体、大屏展示这些锦上添花的东西。个人实操下来最冤的情况就是辛辛苦苦下载数据、清洗合并结果因为少设置一个survey.lonely.psu选项全部分析不能跑。我现在的习惯是每次新建R脚本第一件事就是把options(survey.lonely.psu adjust)和library(survey)写在最前面再开始写业务逻辑。这个微不足道的小习惯能帮你省下大量查报错的时间。NHANES分析并不难难点在于你愿不愿意把抽样设计当成第一公民来看待而一旦你掌握了survey包这套工作流以后处理任何复杂抽样数据库都是同一个套路一劳永逸。
