我和大家一样，最开始拿到NHANES数据时，第一反应就是读进来、跑个t检验、画个图，完事了。后来有一次投稿，审稿人直接问我：“你的加权权重在哪里？调查周期怎么处理的？为什么结果和官方发布的患病率差了这么多？”当时真是被问住了。从那以后我才认认真真把NHANES的加权分析从头捋了一遍，也才意识到：不搞清楚权重，NHANES的分析结果基本就是自说自话，上不了台面。
这篇内容就是把我自己实际跑通的完整流程分享出来——从数据读取、权重选择、survey设计对象构建，到加权均值/占比计算、分组比较、回归建模，再到把结果可视化，全程附上可直接运行的R代码。5分钟跑通不敢说人人做到，但只要你照着走一遍，整个逻辑就在脑子里了。适合刚接触NHANES的医学研究生、做流行病学分析的同学，以及任何被“复杂抽样加权”这四个字折磨过的人。
1. 为什么NHANES分析必须谈加权：一个算偏的真实场景
先不急着上代码，我想先讲一个真实案例，因为这个案例直接决定了你后续每一步怎么写。
我朋友之前做某生化指标与代谢综合征的关联分析，用的就是NHANES某一周期的数据。她当时为了省事，直接把所有样本当成随机样本，简单算了个患病率，结果算出来比美国CDC官网公布的数值低了差不多3个百分点。她又复查了一遍数据清洗，没发现问题，最后才发现是权重没加。
NHANES不是简单随机抽样，而是分层多阶段概率抽样。它会有意过抽某些亚群，比如老年人、墨西哥裔美国人、孕妇、青少年等，以确保这些亚群的样本量足够做统计推断。这就意味着，你手里的每一行并不是“代表一个人”，而是“代表一定数量的美国平民人群”。如果你不给每行赋予对应的“代表人数”，那过抽的亚群就会在你的结果里话语权过大，欠抽的亚群就被边缘化。
权重的作用就是把这个“代表数量”还原回去。加权后算出来的均值、患病率、回归系数，才能代表全美平民人口，而不是代表你这个样本本身。
NHANES的复杂抽样设计，主要包含三块：分层（Stratum）、整群（Cluster/PSU）、权重（Weight）。
分层变量在公开数据里通常是SDMVSTRA，PSU变量是SDMVPSU，权重变量则根据你要分析的数据类型不同而不同。很多人第一次接触这些变量名会懵，其实没关系，你只要记住：凡是做描述性统计和回归，必须声明它们；凡是涉及多周期合并，权重还需要做调整。
所以，不要嫌麻烦。你想让你的结果能进审稿人的法眼，第一步必须把survey设计对象建对。
1.1 NHANES权重类型：你应该选哪个
NHANES公开数据里有几个比较常见的权重变量，最常用的是这两个：
WTMEC2YR：MEC（Mobile Exam Center）检查权重，适合使用体检、实验室检测数据的分析。WTINT2YR：访谈权重，适合只使用问卷调查数据的分析。
如果你同时用了问卷数据和实验室数据，通常建议使用WTMEC2YR，因为样本量会收缩到那些真正做了体检的人，这样分母才一致。
如果合并了多个周期（比如2015-2016、2017-2018两个周期合并），那权重需要做如下调整：
- 2个周期合并：权重 =
WTMEC2YR / 2或WTINT2YR / 2。 - 3个周期合并：权重 =
权重 / 3。 - 4个周期合并：权重 =
权重 / 4。
原因是原始权重代表的是“该个体代表2年内的多少人”，你合并了多个2年周期后，如果不除以周期数，总人数就会虚高，方差就会被低估，p值也会偏小。这一点很多人会忽略，但审稿人很爱问。
2. 环境准备：R包选型与数据导入细节
工欲善其事，必先利其器。NHANES加权分析在R里最核心的包是survey包，这个包几乎是此类分析的标配。配套的还有tidyverse（数据清洗）、janitor（快速生成四格表/百分比）、gtsummary（生成发表级表格）、ggplot2（可视化）。
下面是我实际敲过的加载代码，你可以直接用：
library(tidyverse) library(survey) library(janitor) library(gtsummary) library(ggplot2)有人会问，nhanesA这个包不是能直接下载NHANES数据吗？为什么不用？我的习惯是：nhanesA的问题在于它偶尔会因为网络或服务器端接口变动出现下载失败，而且你还会被变量名搞晕。我更推荐直接去CDC官网按周期下载XPT或CSV格式的数据文件，存到本地，用read_xpt或read_csv读入。这样整个流程更稳，也不会受制于网络。
2.1 从CSV到survey设计对象：完整链路
假设你已经下载好了某周期的 demographics（人口学）数据和 examination（实验室/体检）数据，并且用SEQN（序列号）做了合并。读取和合并的代码大致是：
# 读取 demog <- read_xpt("DEMO_I.XPT") exam <- read_xpt("BIOPRO_I.XPT") # 按SEQN合并 nhanes_df <- demog %>% left_join(exam, by = "SEQN")合并之后，需要做一件非常重要的事：构建survey设计对象。这是后面所有加权计算的基石。
nhanes_design <- svydesign( id = ~SDMVPSU, # PSU，整群标识 strat = ~SDMVSTRA, # 分层标识 weight = ~WTMEC2YR, # MEC检查权重 nest = TRUE, # 一定要设为TRUE，因为PSU编号在层内是重复的 data = nhanes_df )nest = TRUE这个参数是新手最容易忽略的。NHANES的PSU编号在每个层内是独立编号的，换句话说，不同层的PSU编号会有重复。如果不加nest = TRUE，survey包会把不同层的相同PSU编号误当成同一个集群，标准误就会算错。
2.2 构建survey对象前先做数据清洗
构建survey对象前，建议先做几件琐碎但必要的事情：
- 年龄变量
RIDAGEYR已经是数值型了，但如果你要用年龄段分组，建议提前切好：
nhanes_df <- nhanes_df %>% mutate(age_group = cut(RIDAGEYR, breaks = c(0, 18, 40, 60, 80), right = FALSE, labels = c("<18", "18-39", "40-59", "60-79", "80+")))- 性别变量
RIAGENDR是1/2编码，建议转成因子：
nhanes_df <- nhanes_df %>% mutate(gender = factor(RIAGENDR, levels = c(1, 2), labels = c("Male", "Female")))- 别忘了排除关键变量的缺失值。抽样设计对象建好之后，缺失值会像滚雪球一样影响后续计算，所以最好在进
svydesign之前就把你分析所需的变量缺失行筛掉。
3. 加权描述统计：均值、患病率、分组对比
survey对象建好之后，加权描述就很简单了。通用函数是svymean、svytotal、svyby、svyglm。
3.1 加权连续变量汇总
假设你要计算某个血清指标的加权均值，代码长这样：
# 单变量加权均值 svymean(~LBXSCR, nhanes_design, na.rm = TRUE) # 按性别分组的加权中位数 svyby(~LBXSCR, ~gender, nhanes_design, svyquantile, quantiles = 0.5, na.rm = TRUE)svyby是个非常实用的函数，它的逻辑是：按某个分类变量分组，再对每组调用你指定的统计函数。不用写循环，一行就出分组结果。
3.2 加权患病率/百分比
分类变量的话，比如算一下代谢综合征的加权患病率：
# 假设disease是0/1变量，1代表患病 nhanes_df <- nhanes_df %>% mutate(disease = factor(disease, levels = c(0, 1), labels = c("No", "Yes"))) nhanes_design <- svydesign( id = ~SDMVPSU, strat = ~SDMVSTRA, weight = ~WTMEC2YR, nest = TRUE, data = nhanes_df ) svymean(~disease, nhanes_design, na.rm = TRUE)如果你想快速得到一个带95%置信区间的表格，用svyby配合confint：
props <- svyby(~disease, ~gender, nhanes_design, svymean, na.rm = TRUE) props confint(props)这里插一句经验：手动计算加权比率的置信区间时，很多人直接在均值上加1.96倍标准误，这在样本量比较大的时候没问题，但复杂抽样设计下的自由度可能小于传统正态近似的要求，特别是某些亚组样本量小时，建议用confint生成的基于t分布的置信区间，会更稳妥。
3.3 加权卡方检验与组间比较
组间比较时，你不能直接用chisq.test，因为那个函数完全不认识复杂抽样设计。你需要用svychisq：
svychisq(~gender + disease, nhanes_design)如果是比较两个加权连续变量的组间差异，可以用svyttest，它是加权版的t检验：
svyttest(LBXSCR ~ gender, nhanes_design)我实际用下来的体会是，加权t检验和普通t检验在多数情况下结论方向一致，但p值和置信区间宽度有差异。当你样本量分布不均衡时，这个差异会被放大。所以养成习惯——只要是NHANES数据，就一律用svyttest，不要用t.test。
4. 加权回归模型：逻辑回归的实际写法
做关联分析时，加权逻辑回归是高频操作。我见过不少人用glm(disease ~ age + gender + BMI, data = nhanes_df, family = binomial())直接跑，这其实是错的，因为它忽略了抽样设计，得到的标准误会偏小，假阳性风险变高。
正确的做法是用svyglm：
model <- svyglm( disease ~ age_group + gender + bmi, design = nhanes_design, family = quasibinomial() # 注意，复杂抽样下一般用quasibinomial ) summary(model)为什么用quasibinomial()而不是binomial()？这是survey包的一个重要细节。由于复杂抽样下观察值并非完全独立，准二项分布族可以允许一定程度的过度离散，得到的p值更稳健。如果你用binomial()，有时候会看到类似“non-integer #successes”的报错，改为quasibinomial()通常就正常了。
如果想提取OR值和95%置信区间：
# 转成数据框 model_tidy <- broom::tidy(model, conf.int = TRUE, exponentiate = TRUE) model_tidyexponentiate = TRUE会直接输出OR值和置信区间，方便你写论文结果。
4.1 连续变量还是分类变量：回归变量编码经验
在实际分析中，年龄既可以用连续变量，也可以分段。如果你用连续年龄，输出的是每增加一岁对应的OR值，这可能不太好解释。我更建议把它分成年龄组，这样结果与流行病学文献的习惯更一致。
另一个技巧是：如果你需要做趋势性检验（比如看教育水平与患病风险的趋势），可以把有序分类变量转成数值型再放进模型：
nhanes_df <- nhanes_df %>% mutate(edu_num = as.numeric(DMDEDUC2))这样得到的是每上升一个教育等级对应的OR，能直接支撑“趋势性关联”的结论。
5. 加权结果的可视化：从表格到能发表的图
可视化部分是我最喜欢的环节。很多人担心survey对象不能直接接ggplot2，其实思路很简单：先用survey函数算出加权结果，提取到普通数据框里，再交给ggplot2画图。你在图中展示的不是原始样本，而是加权后的代表总体的估计值。
5.1 加权患病率的条形图
比如你想画不同性别/年龄组的患病率条形图，还带误差条，代码如下：
# 先算加权患病率和置信区间 prevalence <- svyby( ~disease, ~gender + age_group, nhanes_design, svymean, na.rm = TRUE ) %>% as_tibble() %>% rename(prevalence = diseaseYes) %>% mutate( lower = prevalence - 1.96 * se, upper = prevalence + 1.96 * se ) # 再画图 ggplot(prevalence, aes(x = age_group, y = prevalence, fill = gender)) + geom_col(position = position_dodge(0.9), width = 0.7) + geom_errorbar( aes(ymin = lower, ymax = upper), position = position_dodge(0.9), width = 0.2 ) + labs(x = "Age group", y = "Weighted prevalence", fill = NULL) + theme_minimal(base_size = 14)这里有个易错点：svyby返回的列名取决于你的结果变量水平名。比如disease是0/1编码时，返回的列名可能是disease0和disease1。建议先glimpse()看一眼再rename，不然很容易对不上。
5.2 加权连续变量分布图
如果你想看某个指标的加权分布，可以把它和调查权重一起传入ggplot的aes：
nhanes_df %>% filter(!is.na(LBXSCR), !is.na(WTMEC2YR)) %>% ggplot(aes(x = LBXSCR, weight = WTMEC2YR)) + geom_histogram(bins = 30, fill = "#377EB8", color = "white") + labs(x = "Serum creatinine", y = "Weighted count") + theme_minimal(base_size = 14)在geom_histogram里设置weight = WTMEC2YR，画出来的直方图就是加权过的人数分布，能反映总体的形态，而不是样本形态。这一个参数的变化，往往就能改变你对数据的直观认知。
5.3 回归结果的森林图
如果你做了回归模型，最直观的可视化方式是画一个森林图，展示各个变量的OR值和置信区间。这个图改编自forestplot包或ggplot2手绘都可以，我习惯用ggplot2做，因为和上面图风格统一。
model_tidy %>% filter(term != "(Intercept)") %>% ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + geom_pointrange() + geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") + coord_flip() + scale_y_log10() + labs(x = NULL, y = "Odds ratio (95% CI)") + theme_minimal(base_size = 14)这里我把y轴取了对数，因为OR值是倍数关系，取对数后对称才更直观。
6. 绕开那些坑：我在实际运行中遇到的报错与处理
跑NHANES的这几年，我遇到了不少让人抓狂的报错。有几次卡了我几个小时，最后发现原因特别蠢。这里把最常见的几个写出来，希望你能绕过去。
6.1 “missing values in 'weight'”报错
如果你在构建svydesign时没筛缺失，survey包会直接罢工。它的逻辑是：如果权重缺失，整行都不能进入设计对象。解决办法是在svydesign之前的nhanes_df里把所有分析相关变量的缺失行筛掉：
nhanes_df <- nhanes_df %>% filter(!is.na(WTMEC2YR), !is.na(LBXSCR), !is.na(RIAGENDR))重要提醒：如果你分析时临时加了一个新变量，而这个变量有缺失，最好重新筛一遍，再重建survey对象。survey对象不像普通数据框那样会自动按新变量过滤。
6.2 “non-integer #successes in a binomial glm”警告
这个警告出现的原因很简单——你的因变量不是纯0/1，或者权重是小数。解决方式就是我前面提到的，在svyglm里用family = quasibinomial()。如果你坚持要整数的概念，可以把因变量转换成0/1因子，但还是建议直接用quasibinomial，这是survey数据建模的正常姿势。
6.3 自由度(degf)太低的问题
有时候亚组分析会报自由度太低，尤其当你按两个变量交叉分组后，某些层只剩一个PSU。此时survey没法稳定地估计方差。我的建议是：要么减少分组变量，要么换个分析粒度。千万不要硬算，因为标准误失真后，p值没有任何参考意义。
6.4 多周期合并时权重没除周期数
前面提过，合并两个周期后权重是WTMEC2YR / 2。如果你只合并了一个新变量但忘了做除法，方差会被低估，p值也会偏小。审稿人如果问“你怎么处理的调查周期”，你要能清晰说出来。
7. 我的写码习惯与结语前的一点建议
最后分享几个我自己的习惯，不一定适合所有人，但确实帮我少走了很多弯路。
第一，我会把分析流程拆成三个脚本：01_clean_data.R（读数据、清洗、构建survey设计对象）、02_descriptive.R（算加权描述统计和分组比较）、03_model_and_figures.R（跑回归出图）。这样每次重新分析时，不需要从头跑数据合并，很省时间。
第二，我习惯在分析前先把变量名存成一个向量，方便统一调整：
vars <- c("RIAGENDR", "RIDAGEYR", "RIDRETH1", "DMDEDUC2", "LBXSCR", "WTMEC2YR") nhanes_df <- nhanes_df %>% select(all_of(vars))这样如果后边要加变量，只需要改一处。
第三，每次跑完分析，我会顺手把R环境和包版本存下来，用sessionInfo()导出。这个习惯不仅方便自己复现，写补充材料时也可以贴给期刊。
NHANES的加权分析，本质上就是两件事：把survey设计对象建对，然后让后续所有函数都走survey这条路。只要这两点做到位，结果就是站得住脚的。希望这份代码能帮你少踩一些坑，把时间花在解读结果和打磨论文上，而不是耗在报错信息里。
