1. 为什么NHANES数据库值得你花心思去啃
先说个我经常遇到的场景:身边不少做临床、公共卫生、营养研究的朋友,手头有一堆好问题,却苦于没有数据,只能发问卷、跑医院收集样本,折腾大半年数据量还不够。其实很多时候,一个免费、公开、持续更新了六十多年的数据库就摆在那里——NHANES数据库,全称National Health and Nutrition Examination Survey,美国国家健康与营养调查。它是一个横断面调查项目,由美国疾病控制与预防中心(CDC)旗下的国家卫生统计中心(NCHS)负责执行,每两年一个周期,收集全美人口的营养状况、健康指标、生活方式、环境暴露、遗传信息等海量数据。
我第一次用NHANES数据库的时候,第一反应是“这数据也太散了吧”。它不像很多现成数据集那样给你一张干干净净的宽表,而是把人口学信息、问卷结果、身体检查、实验室检验分门别类存成几十个独立文件。这意味着,不管你想做什么分析,几乎都绕不开数据合并这一步。而合并完之后,如果直接拿原始数据跑统计分析,结果很可能是有偏的——因为NHANES数据库的抽样设计是复杂的多阶段分层概率抽样,不是简单随机抽样。你必须用上它提供的权重变量和分层、聚类信息,才能得到真正能代表全美非机构化人口的估计值。
这篇内容就是冲着这两个核心痛点去的:一是把数据合并这件事讲透,二是把加权分析的门道说明白。适合正在用或准备用NHANES数据库做论文、做课题、做课程设计的医学生、研究生、临床医生和公卫从业者。我会用R语言做演示,因为R的survey包处理这类复杂抽样数据非常成熟,当然部分思路和原理同样适用于Stata和SAS。
2. NHANES数据库的底细:文件结构、周期与读取方式
2.1 两年一个周期,文件命名是有规律的
NHANES数据库从1999年开始采用连续两年为一个周期的模式,比如1999-2000、2001-2002,一直到最近的2017-2018、2021-2023(COVID期间特殊)。每个周期内,所有参与者会先接受家庭访谈问卷(Interview),然后其中一部分人被邀请到移动检查中心(MEC,Mobile Examination Center)做体格检查和实验室检测。这个设计产生了两个不同的“样本底子”:一个是完成家庭访谈的全样本,一个是进一步完成MEC检查的样本。这两类样本对应的权重变量不同,后面我会专门讲。
数据下载页面把文件按数据域分类,比如Demographics(人口学)、Dietary(膳食)、Examination(检查)、Laboratory(化验)、Questionnaire(问卷)。每个文件名的结构通常是字母缩写加上年份区间,例如DEMO_X代表2015-2016周期的人口学文件,P_ALB_CR代表2017-2018周期的尿白蛋白肌酐比文件。用多了你就会发现,文件名里的前缀(如DEMO、BMX、ALQ、DR1TOT)和调查板块之间有明确的对应关系,查一下数据文档里的Variable List就能对上号。
2.2 xpt格式怎么读取
NHANES数据库提供SAS transport格式(.XPT)文件下载,R语言里可以用两种方式读取:一个是foreign包里的read.xport函数,是老牌工具但偶尔会在变量标签上出问题;另一个是haven包里的read_xpt函数,它读取SAS格式的体验更顺滑,变量标签、格式信息保留得更好。我个人强烈推荐haven包。
library(haven) demo15 <- read_xpt("DEMO_I.XPT") alq15 <- read_xpt("ALQ_I.XPT")这里的文件后缀I表示这是2015-2016周期。拿到手之后先别急着做别的,第一件事永远是查看变量名和标签:
library(dplyr) glimpse(demo15)这一步能帮你快速了解文件里有哪些候选变量。NHANES数据库的变量名有时很反直觉,比如体重变量叫“BMXWT”,腰围叫“BMXWAIST”,如果你不看数据字典,单靠猜,很容易在合并之后发现变量根本不存在。我在早期使用时就吃过这个亏,拿着一个旧周期的变量名去匹配新周期,结果当然是找不到。
2.3 别忽略数据字典文档
每个周期的主页面都会附带一份数据字典(Codebook),以变量名、变量说明、值标签、缺失编码的形式逐条列出。下载数据文件时一定要顺手把对应字典存下来。NHANES数据库的变量命名有继承性,多数变量跨周期保持一致(比如年龄RIDAGEYR、性别RIAGENDR),但不同周期的编码方式可能有变化,尤其是一些问卷类变量。举个例子,吸烟状态变量SMQ020在不同周期里的值标签措辞有细微差异,如果不查字典直接合并,很容易把“现在每天吸”和“现在偶尔吸”当成同一个东西处理。
我的习惯是为每个项目单独建一个变量字典表,列出变量名、所属文件、周期、取值含义,分析的时候随时对照。这个看似笨办法的习惯,能让你在合并数据时少走很多弯路。
3. 数据合并的核心操作:从一维表拼出你的分析数据集
3.1 合并的主角:SEQN
NHANES数据库里的每一个被调查者都有一个唯一的受访者编号,变量名为SEQN(Respondent sequence number)。所有文件——无论是人口学、问卷、检查还是化验——都通过SEQN关联。这其实就是数据库设计里的主键概念。
合并操作按方向分就两类:横向合并(增加列)和纵向合并(增加行)。在NHANES数据库场景里,把不同文件按SEQN拼到一起是横向合并,也就是把多个文件的变量拼到同一个人的同一行记录上;而当你需要把多个周期放在一起分析时,则要纵向合并,也就是把不同周期的行堆叠起来。
处理NHANES数据库合并最舒服的方式是用dplyr包里的`full_join`或`left_join`系列函数。我通常会以人口学文件(DEMO)为主表,然后用left_join把其他文件依次并进来。
merged_data <- demo15 %>% left_join(alq15, by = "SEQN") %>% left_join(bmx15, by = "SEQN") %>% left_join(lab15, by = "SEQN")这里使用left_join而不是inner_join的原因后面会展开,但先记住一条核心原则:人口学文件里的每一行代表一个受访者,问卷和化验文件里的行数不一定是每个人一行,有些检验会有多次测量(比如不同时间点测了两遍血压),所以合并时一定要先检查一下对应文件里SEQN是否为唯一值。
3.2 合并前先查重复SEQN:多对多合并的隐患
这是新手最常踩的坑。如果某个化验文件对同一个SEQN出现了多行记录(例如某个指标测了两次),而你直接left_join上去,那么主表中这一行会被复制成多行,样本量被悄悄放大。这种错误在后续加权分析时很难被察觉,但结果已经严重偏差了。
所以每次合并新文件前,我都会跑一下这段代码确认:
lab15 %>% count(SEQN) %>% filter(n > 1) %>% nrow()如果返回0,说明这个文件里每个受访者只有一行,可以安全合并。如果返回大于0,就得回到文档里看重复记录代表什么。有些文件本身就设计成长格式(long format),比如膳食回忆数据,一个人一天可能吃了很多种食物,每种食物一行。这种文件在合并前需要先做数据重塑(例如把某种营养素的摄入总量按人聚合),再合并到主表。
3.3 只有访谈数据的样本和做完MEC检查的样本——合并范围要分清
前面提到,NHANES数据库把参与者分成两类:完成家庭访谈的(Interview sample)和进一步完成MEC检查的(MEC exam sample)。并不是每个参与访谈的人都会去做MEC检查,所以不同数据文件对应的样本底子不一样。
当你把一个问卷文件和一个化验文件合并到人口学文件上时,你会发现合并后某些变量的缺失值变多了,这是正常的。因为问卷是全样本的,而化验只有MEC子样本才有。后续分析如果同时用到问卷变量和化验变量,你的有效分析样本就自动收缩到MEC子样本,此时必须使用MEC权重(WTMEC2YR),而不是访谈权重(WTINT2YR)。这个选择直接关系到你的结果能不能代表美国全国人口,后面加权部分再细说。
3.4 跨周期合并:权重变量名冲突怎么办
如果你的研究需要更大的样本量,就可能要把多个周期的数据纵向拼接起来。这时候会遇到一个很实际的问题:每个周期都有自己的权重变量,名字都叫WTMEC2YR,直接合并会造成变量重复命名。
我的处理方法是先给每个周期数据加一列周期标识,然后重命名权重变量:
demo15 <- demo15 %>% mutate(cycle = "2015-2016") %>% rename(WTMEC = WTMEC2YR, WTINT = WTINT2YR, SDMVPSU = SDMVPSU, SDMVSTRA = SDMVSTRA) demo17 <- demo17 %>% mutate(cycle = "2017-2018") %>% rename(WTMEC = WTMEC2YR, WTINT = WTINT2YR, SDMVPSU = SDMVPSU, SDMVSTRA = SDMVSTRA) combined_demo <- bind_rows(demo15, demo17)这里的SDMVPSU是初级抽样单元(PSU)编号,SDMVSTRA是分层变量。如果你合并两个周期,后续计算权重时要记得把权重除以2;合并三个周期就除以3。原因是每个周期原本都代表全国人口,两个周期堆在一起后人口总数翻倍了,权重必须统一缩放到能代表两年平均人口的水平。
4. 加权分析的底层逻辑:为什么不能直接算平均值
4.1 NHANES数据库的抽样设计决定了你必须加权
NHANES数据库不是从美国人口里均匀随机抽人的。它的抽样策略是多阶段分层概率抽样:先在全国范围内抽取若干个县(初级抽样单元),再在县内抽取街区,然后在街区内抽取住户,最后在住户内抽取个人。同时,为了提高某些亚群(如老年人、墨西哥裔美国人、孕妇、低收入人群)的估计精度,调查设计会有意地超采样(oversampling)这些群体。
这就带来一个后果:数据里一个人的背后代表的人数并不一样。如果一个群体被超采样了,那么它的样本占比就比真实人口占比高。你不加权,算出来的任何率、任何均值都只代表“样本里的分布”,而不是“全美人口的分布”。举个例子,如果设计时超采了老年人,那么你没加权算出来的高血压患病率会高于真实值,因为样本里老年人太多了。
权重变量就是用来修正这个问题的。它本质上是一个“放大系数”,表示这个受访者代表了多少个实际人口。把每个受访者的权重加总,约等于该周期美国的总人口数。
4.2 权重变量到底怎么选:WTMEC2YR还是WTINT2YR
这是我在社区里见到被问得最多的问题。简单说:
- 如果你只用了问卷变量(只有访谈数据),用WTINT2YR(访谈权重)。
- 只要你用到了任何MEC检查或化验数据,就必须用WTMEC2YR(MEC权重)。
实际应用中,绝大多数分析都会用到化验指标(比如血糖、血脂、维生素D),所以WTMEC2YR是出场率最高的权重。但不要无脑用——如果你的研究只用年龄和收入这些问卷变量做分析,用WTMEC2YR虽然不会造成灾难性后果(因为权重结构相似),但并不是最规范的用法。
除了这两个全局权重,NHANES数据库还有一些特殊模块的专属权重。比如膳食数据里,第一天饮食回顾的权重是WTDRD1,第二天是WTDRD2;某些子研究(如肾功能研究)会有单独的权重变量。另外还有为合并多个周期设计的权重。记住一个原则:每个数据文件里面带Wt开头的变量,都是给这个特定模块用的,选择时要先看文档说明。
4.3 声明调查设计对象:R语言survey包的正确姿势
加权分析在R里最常用的实现方式是survey包。它的思路是先用`svydesign`函数声明“我的数据是怎么抽出来的”,之后所有统计分析函数(如svymean、svytotal、svyglm)都会自动考虑分层、聚类和权重。
library(survey) nhanes_design <- svydesign( id = ~SDMVPSU, # 初级抽样单元 strata = ~SDMVSTRA, # 分层变量 weights = ~WTMEC2YR, # MEC权重 nest = TRUE, data = merged_data )这里的nest = TRUE非常关键,因为NHANES数据库里的PSU编号是在每个层内独立编号的,PSU编号在不同层之间会重复。如果不设置nest = TRUE,survey包会误以为不同层用了同一个PSU,方差估计会出错。
我见过不少人在这一步偷懒,直接用某个现成代码模板,结果没注意PSU编号的问题,后面算置信区间时发现区间窄得离谱。方差估计是加权分析里水最深的地方之一,务必小心。
4.4 加权前后到底差多少:一个直观对比
我在给学生演示时经常做一件事:同样的数据,算一下未加权的血清维生素D平均水平,再算一下加权后的平均水平。结果通常相差不小——因为NHANES数据库超采样了维生素D水平普遍偏低的某些亚群(比如非西班牙裔黑人),未加权估计会被拉低,加权之后才接近真实人口水平。
再有就是患病率。很多人用未加权数据算出一个“全国患病率”发到文章里,审稿人一眼就能看出问题。凡是涉及NHANES数据库的描述性统计,都必须报告加权估计值,并在方法部分明确写出使用的权重变量和调查设计声明方式。现在很多期刊对这类报告规范(如STROBE声明)要求很严格,方法描述不完整会被退回。
5. 一个完整的实战示例:从原始文件到加权回归
5.1 研究场景设定
假设我想研究一个很实际的问题:2015-2016年美国成年人中,血清维生素D水平与肥胖指标(BMI)之间是什么关系。这个分析里我至少要合并三个文件——人口学文件(获取年龄、性别)、身体测量文件(获取BMI)、实验室文件(获取血清维生素D,25-hydroxyvitamin D)。
5.2 读取和合并数据
library(haven) library(dplyr) library(survey) demo15 <- read_xpt("DEMO_I.XPT") bmx15 <- read_xpt("BMX_I.XPT") vid15 <- read_xpt("VID_I.XPT") analysis_data <- demo15 %>% select(SEQN, RIDAGEYR, RIAGENDR, RIDRETH1, WTMEC2YR, SDMVPSU, SDMVSTRA) %>% left_join(bmx15 %>% select(SEQN, BMXBMI), by = "SEQN") %>% left_join(vid15 %>% select(SEQN, LBXVIDMS), by = "SEQN") %>% filter(RIDAGEYR >= 20) summary(analysis_data$BMXBMI) summary(analysis_data$LBXVIDMS)这里我用了`select`提前把需要的变量选出来,而不是把整个文件都并进来。这能让数据框保持精简,避免变量名冲突,也让后续操作更流畅。先用`summary`看一下缺失情况,如果某个变量缺失太多,就要查查是文件合并问题还是样本本身的缺失设计。
5.3 声明调查设计并跑加权统计
nhanes_design <- svydesign( id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, nest = TRUE, data = analysis_data ) # 加权平均BMI和维生素D svymean(~BMXBMI, nhanes_design, na.rm = TRUE) svymean(~LBXVIDMS, nhanes_design, na.rm = TRUE) # 加权线性回归 model <- svyglm(LBXVIDMS ~ BMXBMI + RIDAGEYR + factor(RIAGENDR), design = nhanes_design) summary(model)`na.rm = TRUE`在处理缺失值时很好用,但在做回归时要小心,`svyglm`的缺失值处理方式是整行删除,也就是说只要任何变量缺失,这一行就直接被排除。在NHANES数据库这种合并场景里,样本量的减少可能很可观。比较好的做法是在建模前先检查一下有效样本量,比如:
analysis_data_complete <- analysis_data %>% filter(!is.na(BMXBMI), !is.na(LBXVIDMS), !is.na(RIDAGEYR), !is.na(RIAGENDR)) nrow(analysis_data_complete)5.4 亚组分析时的加权注意事项
如果你只想分析女性,或者只想分析糖尿病人群,最直接的做法是filter之后重新声明svydesign。在这个子集上,权重仍然有效,因为子样本的抽取机制仍然服从于整个调查设计。但要注意,子集的样本量可能变得很小,这会导致估计方差不稳定,置信区间变得很宽。
有一种更精细的做法是用`subset`函数配合survey设计对象,这样能保留原始设计信息中的层与PSU结构,得到的方差估计在统计上更严谨:
female_design <- subset(nhanes_design, RIAGENDR == 2) svymean(~LBXVIDMS, female_design, na.rm = TRUE)我个人在实际使用中偏好这种做法,因为它的计算过程会自动处理某些只有一两个观测值的层(single-PSU stratum)问题,代码上也更干净。
6. 时间趋势分析和多周期合并的正确姿势
6.1 为什么要做多周期合并
单个周期的样本量在7000到10000人之间,对于常见疾病和指标的描述性分析通常够用,但如果你想做罕见指标的估计、某些小亚群的分析,或者想做时间趋势(比如“近十年肥胖率怎么变”),就必须把多个周期合在一起。
多周期合并时的核心变化是权重调整:如果你合并了n个周期,每个周期的MEC权重就要除以n。原因是每个周期的权重分别把样本放大到了该周期对应的全国人口,合并后如果不除以周期数,加权总人口就成了真实人口的n倍。虽然对点估计本身没有影响(分子分母同时放大),但对方差估计有直接影响,因为survey包要根据权重总和来调整估计精度。
6.2 多周期合并的完整代码模板
demo_1517 <- demo15 %>% mutate(cycle = "2015-2016") %>% rename(WTMEC = WTMEC2YR, WTINT = WTINT2YR) %>% bind_rows( demo17 %>% mutate(cycle = "2017-2018") %>% rename(WTMEC = WTMEC2YR, WTINT = WTINT2YR) ) %>% mutate( WTME_2cycle = WTMEC / 2, WTIN_2cycle = WTINT / 2 ) combined_design <- svydesign( id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTME_2cycle, nest = TRUE, data = demo_1517 )这段模板里我用新的变量名WTME_2cycle存调整后的权重,而不是直接覆盖原列。这样做的好处是万一计算或分析过程中出错,还能回去核对原始权重。
6.3 时间趋势分析中的连续变量处理
如果想比较2015-2016和2017-2018两个周期的某个指标变化,最直接的做法是把年份作为分类变量放进模型:
combined_data <- combined_data %>% mutate(cycle_factor = factor(cycle, levels = c("2015-2016", "2017-2018"))) model2 <- svyglm(BMXBMI ~ cycle_factor + RIDAGEYR + factor(RIAGENDR), design = combined_design) summary(model2)这里的回归系数反映的是两个周期之间的BMI差异(在控制了年龄和性别后)。如果你想更进一步,把调查周期当作连续时间变量(如1999-2000记为0,2001-2002记为1……)来处理,那就能直接得到“平均每个两年期BMI变化多少”的趋势估计。这类分析在NHANES数据库相关的顶刊文章里非常常见。
7. NHANES数据库实战中绕不开的坑与经验性建议
7.1 样本量陷阱:加权后有效样本量不等于行数
很多人看数据框有9800行,就觉得样本量很大。但加权后,由于设计效应(design effect)的存在,你的有效样本量可能只有实际行数的三分之一甚至更少。设计效应反映的是复杂抽样比简单随机抽样效率低的程度。分层通常会降低设计效应(提高精度),而聚类则会增加设计效应(降低精度)。NHANES数据库的聚类效应很明显,因为同一社区里的人往往更相似。
所以当你看到加权估计的置信区间比自己预想的宽很多时,不用惊讶,那可能是真实的设计效应在起作用。写文章时可以考虑报告加权的样本量或者设计效应,这能体现你对方法的理解。
7.2 缺失值处理要高调,不要假装不存在
NHANES数据库的缺失值编码有多种:有些变量用“.”表示缺失,有些用“77777”、“99999”、“Refused”(拒绝回答)、“Don't know”(不知道)表示。在读取xpt文件后,haven包会把一些特殊值转换成NA,但并非所有值都会被自动处理,尤其是一些问卷变量。
读取数据后建议立刻检查变量取值分布,看看最大值是否异常。比如年龄变量RIDAGEYR的取值上限是80(80岁及以上统一编码为80),这是NHANES数据库保护隐私的一种做法,分析时要知道这一点。再比如有些家庭收入变量会用到“$75,000 and over”这样的区间编码,处理时不能当成数值变量直接放进模型。
7.3 加权分析不等于每个步骤都加权
有些人在做倾向性评分匹配、做变量筛选时也用权重,这其实是个误区。权重的用途是让描述性统计和推断统计能够代表目标人口,而建模过程中的变量筛选和特征选择更多关注的是关联结构,过度加权反而可能让你选入仅在大权重个体上显著的噪声变量。我的一般做法是:描述性统计全部加权,回归模型加权,但中间的数据清洗、异常值排查、缺失模式分析不需要加权。
7.4 别忽略隐藏的重复测量
有些NHANES数据库变量是重复测量的:血压测了三次、膳食摄入测了两天、实验室检测有重复样本。如果你做的是“一个人”层面的分析,必须先把重复测量聚合成一个人一个值(比如取平均值),不然合并后会引入严重的相关性偏差。
以血压为例,文件里通常有BPXSY1、BPXSY2、BPXSY3三个收缩压测量值,标准做法是取前两个或三个的平均值(具体看文档),或者用MEC检查里推荐的平均值变量。有些问卷文件直接给出了平均变量,比如BMXBMI就是直接计算好的BMI,这种直接拿来用就可以。
7.5 下载数据时Procedures和Documentation一个都别落
NHANES数据库官网每个周期的数据页里,除了数据文件本身的下载链接,还有一系列配套文档:Analytic Notes(分析注意事项)、Codebook(字典)、Procedure(调查流程说明)。Analytic Notes里经常藏着重要信息,比如某个变量在本次周期改了统计口径、某个检验方法换了实验室、某个子样本权重需要特殊处理。
我强烈建议,每下载一个数据文件,就顺手把它对应的Codebook和Analytic Notes存成一个同名文件夹。不要攒到最后一口气看,因为分析到一半再回去翻文档会打断思路。我个人的整理习惯是:
- 原始数据文件夹按周期分,里面再按板块分子文件夹
- 每个子文件夹放一个README.txt,写明这个文件里的关键变量、样本量、权重变量、特殊编码
- 分析脚本和输出文件单独放,不污染原始数据
这个习惯帮过我很大忙。NHANES数据库的数据文件非常大,变量维度动辄几百上千列,没有良好的文件管理,后期找人复核代码或者自己复盘时会非常痛苦。
8. 把加权分析结果写进论文时的报告要点
很多审稿人对使用公共数据库的研究有额外要求,尤其是方法学部分必须写清楚这几点:
- 数据来源:明确写出使用NHANES数据库哪个周期,数据公开下载的日期。
- 抽样设计声明:说明数据是基于多阶段分层概率抽样,分析时考虑了抽样权重、分层变量和聚类变量。
- 权重选择:明确说明使用WTMEC2YR访谈权重或MEC权重,多周期合并时说明权重调整方式。
- 缺失值处理方式:报告排除了哪些缺失记录,有效分析样本量是多少。
- 统计软件和包:例如R版本、survey包版本。
我见过一些文章因为权重说明模糊而被要求大改甚至拒稿。另一个常见的低级错误是把未加权样本量写进正文,但表格里报告的却是加权后的百分比,弄得读者对不上号。建议在做表时,把所有加权结果和未加权样本量在表格下方注明清楚。
9. 写在最后:一点个人的使用心得
我在处理NHANES数据库的过程中,最大的体会是:这个数据库的价值不只在数据量大,更在于它的设计是开放的、可追溯的、有完整文档的。每一份数据文件背后都有详细的方法学说明,你可以清楚地知道你的分析结果是在什么样的抽样框架下得到的。这比很多封闭的临床数据要强得多。
如果你刚开始接触NHANES数据库,建议不要急着去分析一个大课题。先挑一个单一周期,选一个明确的小问题,走通“下载数据-查看字典-合并文件-声明调查设计-跑加权描述统计-写一个简单回归”这条完整链路。把这条链路跑顺了,再去处理多周期合并、复杂亚组分析、趋势检验这些进阶问题都会顺利很多。
最后再分享一个小技巧:每次分析之前,我会把数据准备阶段做的所有操作写成一个R脚本,从读取原始文件开始,到最终生成分析数据集结束,全程不手工干预。这样一旦发现某个环节出错,只需要改脚本再重跑一遍,不需要从头再来。NHANES数据库的文件结构相对稳定,这套脚本在下一个周期数据出来时往往只需改一下文件名和变量名就能复用,也算是给自己攒下的一笔时间资产。