nhanesA实战：高效整理NHANES数据并完成加权分析

简介nhanesA是一款面向R语言用户的NHANES数据访问工具包旨在帮助医学科研人员、流行病学研究者与数据分析师快速浏览、检索并导入国家健康与营养检查调查数据。包内函数支持自定义选择变量和数据子集需联网使用适合需要高效获取公开健康调查数据的实证研究场景。压缩包共23个文件约40KB包含10个Rd格式函数帮助文档、3个R源码文件构成核心逻辑、包描述与命名空间等元数据文件以及vignettes示例教程、GitHub Actions工作流配置等目录结构完整符合R包开发规范。目前已有1893人学习。借助该包用户可以节省自行解析NHANES数据格式的时间通过现成函数完成数据查询、变量映射与翻译快速搭建基于调查问卷和体检数据的研究数据集也可对照帮助文档逐项学习包内函数用法适用于科研教学与数据挖掘入门实践。 前阵子接了一个活儿用美国国家健康与营养调查NHANES分析血压和血脂的关系。按理说这类公开数据库直接下载就行可真上手才发现光搞定数据文件就很折腾——一堆XPT格式文件散落在不同网页文件名还不规律变量说明又是另一个PDF真等把所有数据拼好半天已经过去了。后来换成nhanesA这个R包十几行代码就把人口学、身体测量、实验室检验这些数据全部拉齐整个分析流程一下子顺了。这篇就把我实际使用nhanesA整理NHANES数据的完整经验写出来从包的思路到实操代码再到踩坑记录给同样被这个数据库折磨过的朋友一个参考。1. NHANES数据为什么难搞以及nhanesA想解决的问题1.1 公共卫生研究常用的NHANES数据结构NHANES是美国疾病控制与预防中心CDC下属国家卫生统计中心NCHS开展的连续性横断面调查主要评估美国成人和儿童的健康与营养状况。这个调查从1999-2000年周期开始之后每两年一个周期数据文件按模块拆得特别细人口学资料DEMO、身体测量BMX、血压BPX、总胆固醇TCHOL、重金属暴露PBCD、膳食问卷DRXTOT等等每个周期少说几十个数据表多的上百个。更麻烦的是文件命名规律。1999-2000周期直接用DEMO.XPT这类名称之后每两个周期追加一个字母后缀2001-2002是DEMO_B2003-2004是DEMO_C一直往后排。所以你想拉2015-2016周期的人口学数据得知道那个周期对应的是DEMO_I。这种命名规则虽然规律性强但需要反复确认尤其不同模块在同一周期的字母后缀是否对齐光靠肉眼检查非常容易出错。1.2 手动下载数据的几个痛点手动下载NHANES数据的流程大概是先去官网找到对应周期的“Questionnaires, Datasets, and Related Documentation”页面一个一个看文件名找到需要的表下载XPT文件再用R的foreign::read.xport或haven::read_xpt读进来。听起来不复杂但实际操作会遇到几个具体问题。一个是变量名分散。你想分析“总胆固醇”和“血压”的关联胆固醇可能藏在TCHOL文件里血压在BPX文件里年龄性别这些协变量又在DEMO文件里。每个文件都有自己的SEQN个体序列号你得先把三个文件按SEQN合并起来。合并本身不难难的是搞清楚每个文件里到底有哪些变量、缺失编码是什么、需不需要剔除不符合条件的样本。另一个痛点是变量说明不在数据文件里。XPT文件读进来就是纯数值比如RIAGENDR性别里1和2代表什么RIAGENDR这个变量在变量说明文档里写着“男性1女性2”。每次都要打开文档对照分析十几个变量就得来回切换页面效率很低。nhanesA这个包就是冲着这些痛点来的它在R里封装了NHANES所有周期数据文件的下载、读取、变量查询、标签翻译、死亡状态数据获取等功能相当于把散落的官方数据源整理成了一个在R内直接调用的接口。对我这种经常要跑数据的人来说它最大的价值不是省了下载这一步而是把“变量发现-数据读取-标签翻译-多表连接”整条链路的效率提上来了。2. nhanesA核心API解析与选型2.1 nhanes()与nhanesFromYear()按周期拉取数据nhanesA最基础的函数是nhanes()直接传入文件名就能拉数据。比如nhanes(DEMO)会读取1999-2000周期的人口学数据nhanes(BPX_C)读取2003-2004周期的血压数据。这个函数背后是直接向CDC服务器请求对应的XPT文件读取后返回data.frame不需要你手动下载和导入。如果你不想记字母后缀可以用nhanesFromYear()按年份取数。比如nhanesFromYear(2015, DEMO)会返回2015-2016周期的人口学文件nhanesFromYear(2017, BPX)返回2017-2018周期的血压数据。我在实际操作中更喜欢用这个函数因为写代码时心里想的是“2015年”不用再去翻一个字母对应哪个周期减少一个容易出错的环节。这两个函数还支持一个重要的参数nhanes(DEMO, translate TRUE)可以在读取数据时自动把编码值转换成带标签的因子。比如RIAGENDR这个变量translate FALSE时是1和2translate TRUE时直接变成Male和Female。这个功能在做探索性分析和出图时非常方便但要注意如果后面要用survey包做加权分析因子变量可能会在某些计算中出问题所以正式建模阶段我一般还是会用原始数值编码自己手动维护一份变量标签。2.2 nhanesVar()与nhanesSearch()变量级检索与抽取NHANES数据文件很大比如实验室检验文件动辄上百个变量全读进来占内存不说还要花时间浏览哪些用得上。nhanesVar()函数可以只读取指定文件中的部分变量相当于在服务器端先做列筛选再返回给你省流量也省内存。比如我想从DEMO_I文件里只取SEQN、RIAGENDR、RIDAGEYR、RIDRETH1这4个变量可以这样写demo_sub - nhanesVar(DEMO_I, c(SEQN, RIAGENDR, RIDAGEYR, RIDRETH1))这个函数在数据量大的场景下特别实用。我有一次拉2017-2018周期的全身营养素数据P_NUTR整个文件几百MB用nhanes()全量读取直接把R跑崩了后来换nhanesVar()只取要用的几十个变量几秒钟就搞定。和变量抽取配合使用的是nhanesSearch()它能在所有周期所有文件中按变量名或描述搜索。比如我想知道“血清铁蛋白”在哪个文件里、叫什么名字可以这样搜search_result - nhanesSearch(c(ferritin, serum)) head(search_result)返回结果里会带变量名、变量描述、所在文件等信息。这个功能对不熟悉NHANES变量命名规则的新手来说特别友好比去官网一个文件一个文件翻要省事得多。2.3 nhanesMortality()直接拿死亡状态数据NHANES数据的一大价值在于它能和一些死亡登记数据关联用于队列分析。以前获取死亡状态数据需要单独去NCHS官网申请下载限制版数据流程比较麻烦。nhanesA包里内置了一个公开的死亡状态数据接口直接传年份就能返回那一年周期的受试者死亡状态、随访时间等信息。mort - nhanesMortality(2015) head(mort[, c(seqn, eligstat, permth_int, mortstat)])这里的变量含义是eligstat表示是否达到死亡状态追踪条件permth_int是从访谈日期到死亡日期或随访截止日期的整数月份数mortstat表示是否死亡。这些数据是做生存分析的核心输入以前要单独处理现在几行代码直接拿到省了不少事。2.4 nhanesQuery()与nhanesTables()数据目录和自定义查询nhanesA还提供了nhanesTables()来查看某个周期的文件清单以及nhanesQuery()来在本地SQLite缓存数据库上执行SQL查询。后者用处在于如果多个变量来自不同文件又想快速筛选可以直接写SQL联表查询。使用nhanesQuery()之前通常先调用一次nhanes()或相关函数让包把数据缓存在本地SQLite库里。然后就能执行类似这样的查询result - nhanesQuery(SELECT seqn, bpxsy1, bpxdi1 FROM nhanes.bpx_i LIMIT 20)需要留个心眼的是nhanesQuery的表名规则和原始文件名不完全一致通常是文件名小写并去掉特殊字符。我自己用的时候会先跑一次nhanesManifest()看看当前缓存里有哪些表再写SQL不然容易报“no such table”的错误。整个包的API大概就是这么几类文件级读取、变量级检索、变量值标签翻译、死亡数据、SQL查询。实际使用时按需组合就行没必要每个都用上。3. 完整实操从取数到合并到加权分析3.1 确定周期与文件清单我用一个具体例子演示完整流程分析2015-2016周期中BMI和收缩压的关系调整年龄、性别、种族并考虑复杂抽样设计。先列清楚需要哪些数据人口学变量年龄RIDAGEYR、性别RIAGENDR、种族RIDRETH1、抽样单元SDMVPSU、抽样层SDMVSTRA、权重WTMEC2YR来自DEMO_I身体测量BMIBMXBMI来自BMX_I血压收缩压BPXSY1、舒张压BPXDI1来自BPX_I因为这些文件都在2015-2016周期统一用nhanesFromYear()拉取library(nhanesA) demo - nhanesFromYear(2015, DEMO) bpx - nhanesFromYear(2015, BPX) bmx - nhanesFromYear(2015, BMX)拉完之后先逐个看维度dim(demo) dim(bpx) dim(bmx)正常情况下三个文件的行数应该一致都是这一周期的受试者总数大概一万人左右。如果行数不一致说明有重复记录或者某个文件包含子表需要进一步检查。3.2 变量筛选与多表合并拿到三个data.frame之后用dplyr做列筛选和按SEQN合并。SEQN是NHANES所有文件共有的个体识别号相当于数据库里的主键library(dplyr) demo_sub - demo %% select(SEQN, RIAGENDR, RIDAGEYR, RIDRETH1, SDMVPSU, SDMVSTRA, WTMEC2YR) bpx_sub - bpx %% select(SEQN, BPXSY1, BPXDI1) bmx_sub - bmx %% select(SEQN, BMXBMI) merged - demo_sub %% inner_join(bpx_sub, by SEQN) %% inner_join(bmx_sub, by SEQN)这里用inner_join是合理的因为做关联分析时只有三个文件都存在的样本才能入组。如果某个问卷文件存在整行缺失就会把样本排掉这也是NHANES分析中常见的有效样本量减少的原因之一。合并后可以快速检查一下数据范围和缺失情况summary(merged[, c(RIDAGEYR, BMXBMI, BPXSY1)])如果发现BPXSY1有大量缺失需要去核对原始血压文件里是否还有额外的测量记录或排除标记而不是直接带着缺失往下走。3.3 变量标签翻译与复杂抽样设计在进入建模前最好给分类变量加上标签方便后续出结果时阅读。nhanesA提供了nhanesTranslate()函数来做这件事merged_labelled - nhanesTranslate(merged, c(RIAGENDR, RIDRETH1)) table(merged_labelled$RIAGENDR)这个函数会在原数据框的基础上增加标签列比如RIAGENDR_f值变成“Male”、“Female”。做交叉表时直接看标签列比记忆1和2要直观得多。不过我用的时候发现偶尔某些变量的标签翻译会失败因为nhanesA内置的变量标签字典未必覆盖所有周期所有变量。遇到这种情况就对单个变量单独处理nhanesTranslate(merged, RIAGENDR)接下来是NHANES分析最重要的部分复杂抽样设计。NHANES不是简单随机抽样需要指定抽样单元、抽样层和抽样权重才能得到全国代表性的估计值。标准写法是library(survey) nhanes_design - svydesign( id ~SDMVPSU, strata ~SDMVSTRA, weights ~WTMEC2YR, nest TRUE, data merged )这段代码里id对应初级抽样单元PSUstrata对应分层变量weights对应两年周期内的抽样权重。nest TRUE是因为在部分周期中PSU编号在层内是复用的必须告诉survey包这个嵌套结构否则会高估标准误。建立设计对象后就可以直接计算加权描述统计了svymean(~BMXBMI, nhanes_design) svyby(~BPXSY1, ~RIAGENDR, nhanes_design, svymean)第一行会输出全人口加权平均BMI第二行会按性别分组输出加权平均收缩压。这些结果加上置信区间就是论文里常见的那张描述性统计表。对一个做了两周期数据的分析来说到这里基本就进入模型阶段了。3.4 多周期数据合并时的权重处理上面的例子只用了一个周期权重直接用WTMEC2YR即可。但如果想增加样本量、做跨周期趋势分析比如合并2015-2016和2017-2018两个周期就需要手动调整权重。NHANES每个周期都有约一万人跨度越长样本量越大但合并时权重不能简单相加。我常用的做法是如果合并k个周期就用每个周期的WTMEC2YR除以k作为新的合并权重。两年周期数据用原权重四个周期合并就除以2六个周期合并就除以3。原因在于每个两年周期的样本量大致相当直接相加会让某个周期因为抽样设计原因占比过大扭曲全国代表性。具体代码是在合并前把权重列除以周期数merged_1718 - nhanesFromYear(2017, DEMO) %% select(SEQN, RIAGENDR, RIDAGEYR, SDMVPSU, SDMVSTRA, WTMEC2YR) merged_combined - bind_rows( merged %% mutate(weight_adj WTMEC2YR / 2), merged_1718 %% mutate(weight_adj WTMEC2YR / 2) )然后把survey设计里的weights参数换成weight_adj。这里有一个容易忽略的点多周期合并时如果要分析的是长期趋势或合并患病率这个近似是够用的但如果对估计精度要求极高建议参考NCHS官方发布的趋势分析权重建议不同分析目标对应的权重处理方式略有差异。4. 常见问题与排查技巧4.1 下载慢或失败怎么办nhanesA本质上是联网请求CDC服务器上的XPT文件所以网络状况直接影响速度。有些实验室文件特别大比如P_NUTR这种几百MB的膳食营养数据下载可能需要几分钟到十几分钟不等。遇到下载中断或者超时可以先手动用浏览器下载对应XPT文件再用haven::read_xpt()读取本地文件不一定非要卡在nhanesA这一个下载通道上。另外我建议把常用数据缓存到本地。nhanesA内部有一个本地SQLite缓存机制你第一次拉过的文件会缓存在临时目录。继续跑分析时如果缓存还在就不会重复从服务器下载。R的临时目录一般重启后会清理如果你长期做NHANES项目可以考虑用选项设置来调整缓存目录避免每次从头下载。4.2 变量名、变量标签和文件名关系千万别搞混nhanesA返回的是data.frame变量名和原始XPT文件里的变量名一致都是大写字母加数字。但nhanesSearch()搜索出来的结果里变量描述和变量名是两回事有时候描述里有关键词不代表变量名也含有关键词反之亦然。我习惯搜索时同时看variable name和description两列避免漏掉自己想要的变量。还有一个隐藏坑同一个变量在不同周期可能含义有细微变化。比如某些问卷变量在2005-2006周期和2011-2012周期的取值编码不一样nhanesA只负责把数据拉下来不会自动帮你校验一致性。多周期合并前务必先对每个周期单独跑一下table()看看分类变量的取值分布是否一致再决定能不能直接合并。4.3 内存不够时的处理思路整表读取在大文件上很容易把内存占满。遇到这种情况优先用nhanesVar()只抽需要的变量而不是全量读取后再用select()筛选。如果变量特别多且分散在多个文件可以用循环配合nhanesVar()逐个文件抽取最后再合并内存占用会小很多。另外R的data.frame对字符型变量默认存储为因子或字符如果变量值本身是纯数字却以文本形式存在会额外占内存。读取后可以检查一下所有列的类型能转数值的尽量转成numeric能降为integer的不要用double这在几百万行的场景下差距很大。4.4 抽样设计里最常见的三个错误第一个是忘记加nest TRUE。很多教程里给的例子直接抄svydesign的模板没有仔细思考嵌套结构结果标准误会偏大。第二个是用错权重变量。不同周期、不同分析目标的权重不一样比如机械采样权重WTMEC2YR用来做整个样本的代表性估计有的子样本需要单独权重如果分析对象是某个亚组必须用对应的子样本权重。第三个是合并周期后直接用原权重不调整这样得到的全国代表性和标准误可能都不对。遇到这三个问题先回看自己的设计对象和权重处理方式通常能定位到原因。nhanesA这个包本身挺好上手的真正让我觉得有价值的是它把NHANES从“手工下载手工查文档”变成“代码调用代码查文档”整个过程可复现、可记录。我后来做其他周期或其他主题的分析基本都是同一套流程套上去只改文件名和变量名分析前的数据准备从半天压缩到半小时以内。如果你也经常跟NHANES打交道建议花一个下午把这几个核心函数过一遍之后每次跑数据都能省下不少时间。本文还有配套的精品资源点击获取