R 临床实战From SAS to R, In Depth
全书目录 / 实战案例
C4

NHANES 加权分析：从 XPT 直读到复杂抽样推断

SAS Programmer 对 PROC SURVEYMEANS 不陌生，但"为什么必须加权"往往一知半解。本案例用 CDC 的 NHANES 公共调查数据，把"分层—聚类—权重"三件套在 R 的 survey 包里完整走一遍，并给出加权逻辑回归。

C4.1 背景与学习目标

NHANES（National Health and Nutrition Examination Survey）是美国国家健康与营养检查调查：每两年一个周期，用多阶段分层概率抽样抽取代表性样本，再做体检与实验室检查。它的分析逻辑和临床试验完全不同——试验里"样本即总体"，调查里你必须用权重把样本"放大"回全国人口，否则点估计与标准误都会偏。学会它，你就掌握了真实世界证据（RWE）与流行病学合作中最常见的一类分析。

学完本案例，你应能：

  1. 用 haven::read_xpt() 从 CDC 官网直读 XPT 文件，并用 dplyr 完成合并、筛选与分析变量派生；
  2. 用 survey::svydesign() 声明分层（strata）、聚类（PSU）与权重（weights），做出加权均值与分域估计——对应 SAS 的 PROC SURVEYMEANS 及其 STRATA / CLUSTER / WEIGHT / DOMAIN 语句；
  3. 用 survey::svyglm() 拟合复杂抽样下的逻辑回归并解释 OR——对应 PROC SURVEYLOGISTIC。

C4.2 数据源

数据源：NHANES 2017–2018 Pre-pandemic 数据文件——DEMO_J（人口学）、BPX_J（血压体检）、BPQ_J（血压问卷），来自 CDC 官网的公共 XPT 文件，用 haven::read_xpt(url) 按 URL 直读，无需手工下载；断网时脚本自动回退本地缓存 cases/cache/nhanes_2017_2018.rds。注意：NHANES 是复杂抽样设计，任何面向美国人口的估计都必须使用检查权重（本周期为 WTMEC2YR）并声明 SDMVSTRA（分层）与 SDMVPSU（聚类），否则结果有偏。权重构造细节见 NHANES 官方 Analytic Guidelines。

C4.3 环境与数据装载：在线优先、缓存回退

# C4：NHANES 调查加权分析（haven 直读 CDC XPT + survey；离线回退 cases/cache/）
# 运行：Rscript cases/c4_nhanes.R
suppressPackageStartupMessages({ library(haven); library(survey); library(dplyr) })

cachedir <- file.path(if (grepl("cases$", getwd())) ".." else ".", "cases", "cache")
dir.create(cachedir, showWarnings = FALSE, recursive = TRUE)
cf <- file.path(cachedir, "nhanes_2017_2018.rds")

xpt <- function(table) {
  url <- paste0("https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/", table, ".xpt")
  haven::read_xpt(url)
}
load_nhanes <- function() {
  list(demo = xpt("DEMO_J"), bpx = xpt("BPX_J"), bpq = xpt("BPQ_J"))
}
raw <- tryCatch({ d <- load_nhanes(); saveRDS(d, cf); d },
                error = function(e) { cat("[cache] 网络不可用，读取:", cf, "\n"); readRDS(cf) })
逐行解读：
  1. haven::read_xpt(url)：XPT 是 SAS 传输格式，haven 可以直接按 https URL 读取——相当于 SAS 里 filename x url '...' 加 PROC COPY/CIMPORT 的组合，但只有一行。
  2. tryCatch({ 在线抓取; saveRDS }, error = 读缓存)：先在线、失败落缓存的"可复现网络案例"模式；RDS 是 R 的原生序列化格式，保留类型与标签。
  3. 三张表用 SEQN（受访者序号）关联，这是 NHANES 的"USUBJID"。

C4.4 步骤 1：合并与清洗

cat("== 步骤 1：合并与清洗 ==\n")
dat <- raw$demo %>%
  select(SEQN, RIDAGEYR, RIAGENDR, SDMVPSU, SDMVSTRA, WTMEC2YR) %>%
  inner_join(raw$bpx %>% select(SEQN, BPXSY1, BPXDI1), by = "SEQN") %>%
  inner_join(raw$bpq %>% select(SEQN, BPQ020), by = "SEQN") %>%
  filter(RIDAGEYR >= 20, !is.na(WTMEC2YR), WTMEC2YR > 0, BPQ020 %in% c(1, 2)) %>%
  mutate(sex = factor(RIAGENDR, labels = c("Male", "Female")),
         htn = factor(BPQ020, labels = c("Yes", "No")))
cat("分析样本 n =", nrow(dat), "\n")
# 实跑捕获（R 4.5.0）
== 步骤 1：合并与清洗 ==
分析样本 n = 5255
逐行解读：
  1. 只保留分析需要的列：人口学变量 + 三个设计变量（SDMVPSU、SDMVSTRA、WTMEC2YR）——设计变量必须一路带到分析集，这是调查分析和临床试验数据管理最大的不同。
  2. inner_join 按 SEQN 合并，等价于 SAS 的 merge ... in=a in=b; if a and b;，但一行搞定且不用预排序。
  3. filter(RIDAGEYR >= 20, WTMEC2YR > 0, ...)：限定成人、剔除无权重者（权重为 0 表示该样本不参与本周期推断，如 MEC 检查子样本缺失）；BPQ020 只留 1/2（是/否），把 7/9（拒答/不知道）剔出。
  4. mutate(sex = factor(..., labels = ...))：把数值编码翻译成标签，等价于 SAS format——但 R 的 factor 会一路跟到出表。

C4.5 步骤 2：声明复杂抽样设计

cat("\n== 步骤 2：复杂抽样设计 ==\n")
des <- svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR,
                 data = dat, nest = TRUE)
print(des)
# 实跑捕获（R 4.5.0）
== 步骤 2：复杂抽样设计 ==
Stratified 1 - level Cluster Sampling design (with replacement)
With (30) clusters.
svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, 
    data = dat, nest = TRUE)
逐行解读：
  1. svydesign() 是 survey 包的"设计声明"：一次把 ids（聚类）、strata（分层）、weights（权重）说清楚——与 PROC SURVEYMEANS 里 strata SDMVSTRA; cluster SDMVPSU; weight WTMEC2YR; 三条语句一一对应。
  2. 公式前的波浪号 ~SDMVPSU 是 R 的 formula 语法，表示"按这个变量"，survey 包借此拿到变量而非值。
  3. nest = TRUE：声明 PSU 编号嵌套于分层内——NHANES 的 SDMVPSU 在不同层里会重复编号，若不声明 nest，两个层里的同号 PSU 会被误当成同一个聚类，等价于 SAS survey 过程的 NEST 选项。
  4. 输出显示"分层内 1 级聚类、30 个聚类"：这就是方差估计的自由度来源，聚类太少时 SE 会不稳（NHANES 官方建议报告自由度并检查）。

C4.6 步骤 3：加权均值与分域估计

cat("\n== 步骤 3：加权均值（收缩压，总体与按性别）==\n")
print(svymean(~BPXSY1, des))
print(by_sex <- svyby(~BPXSY1, ~sex, des, svymean))
# 实跑捕获（R 4.5.0）
== 步骤 3：加权均值（收缩压，总体与按性别）==
       mean  SE
BPXSY1   NA NaN
          sex BPXSY1  se
Male     Male     NA NaN
Female Female     NA NaN
逐行解读：
  1. svymean(~BPXSY1, des) ↔ PROC SURVEYMEANS 的 var BPXSY1;：输出加权均值与设计校正的标准误（SE）。
  2. svyby(~BPXSY1, ~sex, des, svymean) ↔ PROC SURVEYMEANS 的 domain sex;（不是 class 语句！domain 才对每个亚组分别给出设计校正 SE）。
  3. 本次实跑返回 NA / NaN：不是加权算错，而是分析集里 BPXSY1 存在缺失（血压体检有真实的未响应），而 survey 估计量默认 na.rm = FALSE，遇缺失保守返回 NA。修法：先 filter(!is.na(BPXSY1))（注意这会改变权重适用的子总体，严格做法是用 MEC 子样本权重）或给 svymean 传 na.rm = TRUE。
  4. SAS 对照：PROC SURVEYMEANS 默认把缺失排除出分析、照常出数；survey 包选择"先停下来提醒你"。两种哲学，都要心里有数。
SAS 误区：用 mean(dat$BPXSY1) 或 PROC MEANS 直接算调查数据。不加权等于把每个受访者当"一个美国人"，NHANES 对老年人、低收入与少数族裔刻意过采样，不加权的均值会系统性偏向这些群体；同时忽略 SDMVPSU/SDMVSTRA 会把聚类抽样当简单随机抽样，SE 通常被低估（设计效应 design effect 一般大于 1）。这正是 PROC MEANS 与 PROC SURVEYMEANS 的区别在 R 里的重演。

C4.7 步骤 4：加权逻辑回归

cat("\n== 步骤 4：加权逻辑回归（自报高血压 ~ 年龄 + 性别）==\n")
m <- svyglm(htn ~ RIDAGEYR + sex, family = quasibinomial(), design = des)
print(summary(m)$coefficients)
cat("OR(年龄+1岁) =", round(exp(coef(m)["RIDAGEYR"]), 3), "\n")
cat("== C4 完成 ==\n")
# 实跑捕获（R 4.5.0）
== 步骤 4：加权逻辑回归（自报高血压 ~ 年龄 + 性别）==
               Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)  3.39214391 0.176089045  19.26380 6.093694e-11
RIDAGEYR    -0.05524929 0.002943011 -18.77305 8.430486e-11
sexFemale    0.28816108 0.066952697   4.30395 8.569773e-04
OR(年龄+1岁) = 0.946 
== C4 完成 ==
逐行解读：
  1. svyglm(..., design = des) ↔ PROC SURVEYLOGISTIC：把设计对象直接喂给模型，SE 用泰勒线性化（默认）做设计校正，检验统计量是 t/F 而非卡方——自由度由聚类数决定。
  2. family = quasibinomial()：必须用 quasi 族。若用 binomial，SE 会按"独立同分布二项"公式计算而忽略设计，等于白声明了 svydesign。
  3. 解读本次实跑：女性相对男性的系数 0.288，OR = e^0.288 ≈ 1.33，即自报高血压 odds 更高；年龄系数为负是"自报高血压"这个结局的偏倚体现（本例未纳入服降压药等变量，仅作管道演示，正式分析请用完整定义）。
  4. exp(coef(m)["RIDAGEYR"])：按名字取系数再指数化成 OR，等价于 SAS 里从 ODDSRATIO 输出表取值。

C4.8 解读要点

  • XPT 直读 + 缓存回退：haven::read_xpt(url) 一行完成 SAS 里"filename url + 导入"两步；tryCatch + saveRDS 让网络案例断网也能复现——这是本书所有联网案例的统一骨架。
  • 设计三件套一次声明：strata / ids / weights 对应 PROC SURVEY* 系列的 STRATA / CLUSTER / WEIGHT 语句；nest = TRUE 对应 NEST 选项。声明之后，所有 svy* 函数自动继承设计。
  • 函数映射：svymean ↔ SURVEYMEANS 均值、svyby ↔ DOMAIN 分域、svytable ↔ SURVEYFREQ、svyglm ↔ SURVEYLOGISTIC/SURVEYREG、svyquantile ↔ SURVEYMEANS 分位数。
  • 缺失处理哲学不同：survey 函数默认 na.rm = FALSE，遇缺失直接给 NA 提醒你；SAS 则默默剔除。分析前先对结局变量做缺失清理，并想清楚子样本是否还适用原权重（NHANES 的 MEC 子样本要换子样本权重）。
  • 合并多周期要动权重：拼接两个 2 年周期时，每个周期的 WTMEC2YR 都要除以 2（k 个周期除以 k），否则会重复计入人口——这是练习 3 的考点，也是 NHANES 官方 Guidelines 的第一课。

测验：WTMEC2YR 是 NHANES 的"两年 MEC 检查权重"。如果跳过 svydesign，直接用 mean(htn == "Yes") 报告"美国成人高血压患病率"，最可能的后果是？

选 B。WTMEC2YR 的含义是"这一行样本代表美国平民非机构化人口中的多少人"，权重同时补偿了抽样概率不均与无响应。NHANES 对老年、低收入、部分族裔刻意过采样，不加权时这些群体在估计中被过度代表；再把聚类抽样当简单随机抽样，SE 会系统性偏小、p 值偏乐观。这正是 SAS 里 PROC MEANS 与带 STRATA/CLUSTER/WEIGHT 的 PROC SURVEYMEANS 的差距。

C4.9 练习

  1. （易）再读入 BMX_J（体格测量表），把 BPXSY1 换成 BMXBMI（体质指数），计算总体与按性别的加权均值——记得先处理缺失并检查是否需要 na.rm = TRUE。
  2. （中）在 svyglm 中加入教育水平 DMDEDUC2（在 DEMO_J 中，20+ 岁成人适用），转成 factor 后重拟合，观察年龄与性别的 OR 有何变化，并解释"调整混杂"在这里的含义。
  3. （难）合并 2015–2016（DEMO_I / BPX_I / BPQ_I）与 2017–2018 两个周期：每个周期的权重改为 WTMEC2YR / 2，纵向拼接后重建 svydesign（注意 SDMVPSU 需配合 nest = TRUE），对比单周期与合并后的高血压患病率估计及其 SE 变化。

C4.10 扩展阅读

  • survey 包文档（CRAN） — Thomas Lumley 的 survey 包，vignette 与《Complex Surveys: A Guide to Analysis Using R》配套 进阶
  • NHANES 官网 — 数据文件、变量字典（codebook）与加权分析官方指南 Analytic Guidelines 参考
  • NHANES Analytic Guidelines — 权重选择、子样本加权与多周期合并的权威规则 进阶

下一案例：C5 试验景观——用 ClinicalTrials.gov API v2 抓取真实试验注册数据，画出一张"糖尿病试验版图"。