NHANES复现教学 7：加权描述统计，复现论文 Table 1

上一节我们把数据集清洗到了论文口径：3,863 人。但先别急着算结果，NHANES 分析还有一道绕不过去的关：权重。

一句话立论：NHANES 不是"人海采样"，而是一份"设计出来的样本"——每个受访者身上都挂着一个系数（权重），告诉你"这个人背后站着多少个美国人"。不做加权的 NHANES 分析，算出来的患病率、均值、关联都不能代表全美国。

打个比方：NHANES 抽到的 11,933 个人，就像一场"每人票数不一样"的投票——权重就是每张票代表多少美国人：被刻意多抽的人群，一张票代表的人少；没被多抽的人群，一张票代表的人多。加权就是按票的代表人数计票，这样结果才等于"全美国"，而不是"这 11,933 个人"。

本篇要做两件事：讲清权重怎么选、怎么用，然后立刻用它把论文的第一张结果表（Table 1 基线特征表）复现出来，逐行与原文对照。

一、为什么要加权？用我们自己的数据看

第 02 篇讲过 NHANES 是多阶段分层抽样，样本里的人群结构和真实美国长得不一样。失真来自两处：其一，超采样——多数周期会刻意多抽特定人群（历史周期多抽老人、黑人、亚裔；2021–2023 疫情周期取消了种族超采，改为全抽 0–19 岁和 60 岁以上的家庭成员）；其二，即便不刻意超采，多阶段分层抽样本身也会让样本结构偏离人口（不同群体被抽中的基础概率与配合率不同）。

我们用 DEMO_L 实测了一组对比（图1）：样本里 60 岁以上老人占 29.4%，加权还原到全美人口只占 23.3%；样本里 0–17 岁儿童占 31.7%，真实只占 22.6%。2021–2023 周期正好做了官方说的"全抽 0–19 和 60+ 家庭成员"——老人和孩子的比例被人为拉高，所以差距格外明显。

图1：浅蓝是"样本里的人"，深蓝是"加权后代表的人"——差距就是超采样的痕迹。

直接用未加权样本算"全美 60 岁以上人口占比"，你会凭空多算 6 个百分点；算儿童占比会多算 9 个百分点。权重的作用就是给每个受访者一个放大倍数，把被人为扭曲的样本扶正回真实人口。实测权重总和约 3.27 亿，正好对应全美人口总量。

二、权重有哪几种？怎么选

打开 DEMO_L，你会看到好几个 WT 开头的变量（第 02 篇已经预告过：这些是"设计信息"）。它们不是重复，而是对应调查的不同环节——你在哪个环节取数，就用哪个环节的权重：

选择规则一句话：往下取数取到哪一层，就用哪一层的权重——用最"深"环节的权重。 只用问卷 → WTINT2YR；用了体检或化验 → WTMEC2YR；核心变量来自抽血化验 → WTPH2YR。

拿案例论文走一遍判断：

1. 用到入户访谈数据吗？ 用了（PHQ-9 问卷、人口学）；
2. 用到 MEC 体检/化验数据吗？ 用了（血清维生素 D 是 MEC 采血化验）；
3. 血液检测是不是核心暴露？ 是（LBXVIDMS 就是暴露变量）。

→ 结论：落在采血层，首选 WTPH2YR。

这里有一个诚实交代：论文原文只说"按 NHANES Analytical Guidelines 做了适当加权"，没写明用的是哪个权重。所以我们的做法不是猜，而是两种都跑一遍做对照（这就是下一节表格里那两行 WTPH2YR / WTMEC2YR 的意义）。

2.1 一个不做就会算错的细节：权重为 0 的人必须剔除

NHANES 的权重是分环节给的，所以会出现这种情况：没参加某个环节的人，在那个环节的权重上就是 0 或缺失。我们用 WTPH2YR（采血权重）实测成人样本 8,153 人的权重覆盖：

官方指南明确要求：使用某环节权重时，该环节权重为 0 或缺失的个体必须从分析样本中排除。 道理很直白——一个权重为 0 的人代表 0 个美国人，他对"全国估计"没有任何贡献；把他留在样本里只会虚增样本量、干扰部分统计量（比如某些软件会把权重 0 的人算进分母）。

# 建设计对象之前先检查一遍（养成习惯）
sum(is.na(ana$WTPH2YR))            # 缺失多少
sum(!is.na(ana$WTPH2YR) & ana$WTPH2YR == 0)   # 权重为 0 多少
# 官方做法：把无效权重的人剔掉再建设计对象
ana <- subset(ana, !is.na(WTPH2YR) & WTPH2YR > 0)

我们的案例里这一步的结果是"0 人需要剔除"——因为前面的样本漏斗已经把"没有维生素 D 结果的人"筛掉了（那批人恰好就是没采血的人）。实测确认：3,863 人的分析样本里，WTPH2YR 缺失 0 人、为 0 的 0 人。

但别因此跳过这一步：换一个用体检变量的研究（比如用 BMI、血压），若不做这层筛选，就会出现"十几% 的样本权重为 0 却仍在分析里"的情况——这类错误不会报错、只会让结果悄悄偏掉，是 NHANES 分析里最典型的"沉默的坑"。

2.2 官方分析指南（Analytic Guidelines）：NHANES 的"交规"该看哪几条

NHANES 官方有一份 Analytic Guidelines（分析指南），是这类分析的"交通规则"。新手不必通读，但下面这几条必须知道。本系列其实一直在按它做，这里集中交代出处：

建议的读法：把指南当"查手册"用。动手前扫一遍目录知道有什么，遇到具体问题时再翻对应条款；官方也会随周期更新（例如 2017–2020 疫情周期的权重规则就与常规周期不同）。写论文的方法学部分时，写一句"分析遵循 NHANES Analytic Guidelines"是常规做法，但你最好真的核对过上面七条。

三、svydesign：把"设计信息"交给 R（含一个新手必踩的姿势错误）

加权不是"乘以一个数"这么简单。正确的动作是：告诉 R 这份数据的抽样设计长什么样，之后所有统计都基于这个"设计对象"来算，方差、置信区间、P 值才会算对。

library(survey)
des <- svydesign(id = ~SDMVPSU,      # 主抽样单元：先抽县/片区
                 strata = ~SDMVSTRA, # 分层：抽样时的层
                 weights = ~WTPH2YR, # 权重：每个受访者代表多少人
                 nest = TRUE,        # PSU 编号在层内嵌套（官方数据必须设 TRUE）
                 data = ana)         # ana = 上一节清洗好的 3,863 人

四个参数逐个讲清楚：

• id（PSU）：NHANES 先抽"县/片区"再抽人，同一片区的人不是独立的——方差要按"片区"这个层级算；
• strata（层）：抽样是分层的，层内才可比较，方差估计要吃掉这层结构；
• weights：核心，就是上面选出来的那个权重变量；
• nest = TRUE：NHANES 的 PSU 编号在层内循环复用，不加这个参数 R 会把不同层的同号 PSU 当成一个，方差算错。官方数据一律设 TRUE。

下面是新手最容易踩的姿势错误，也是我们必须纠正的一处：想做"只看抑郁人群"的子分析时，不要先把数据筛好再建设计对象：

# ❌ 错误姿势：先筛数据，再建对象
sub <- ana[ana$depression == 1, ]
des_sub <- svydesign(id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTPH2YR, nest=TRUE, data=sub)

# ✅ 正确姿势：先建全样本设计对象，再对"设计对象"取子集
des     <- svydesign(id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTPH2YR, nest=TRUE, data=ana)
des_dep <- subset(des, depression == 1)      # subset 作用在设计对象上

加权之后，置信区间为什么会变宽？ 因为 R 的 survey 包默认用 Taylor 级数线性化（Taylor series linearization） 来估计方差。这名字听着唬人，意思其实很朴素：它不假设每个人互相独立，而是把"同在一个片区、同一层里的人"当成彼此相关的，按这个前提重新算误差。相关信息少了，方差自然更大。（另有一类做法是重复抽样法，如 BRR、jackknife，用在特定场景。）

衡量"复杂抽样 vs 简单随机抽样"差距的指标叫设计效应（deff，design effect）：

deff = 复杂抽样设计下的方差 ÷ 同样本量简单随机抽样的方差

我们实测（分析样本 3,863 人，加权患病率 11.88%）：

但代价是：同一份数据、同一个估计，只要忽略抽样设计，标准误就少算近一半，置信区间跟着变窄、P 值跟着变小，结论看起来比实际更"确定"。这不是小毛病，是结论对不对的问题（第 10 篇做权重敏感性分析时，我们会看到它足以让一个"边缘显著"的结果翻转）。

为什么？ 先把人剔掉，等于把"这些人原来所在的抽样层和片区"信息一起扔掉：某个片区可能只剩 1 个人了，R 只能报错或给出错误的方差（自由度、置信区间全跟着错）。对设计对象取子集，抽样结构才能完整保留。这是 NHANES 官方分析指南点名强调的做法，也是 survey 包设计的本意。记住这条，后面 Table 1 的分组比较、第 09 篇的亚组分析全靠它。

四、读结果的四个基本功：从患病率到 OR，每个数字怎么读

这一节给整个结果部分（本篇与第 08、09、10 篇）打地基。这四个概念在论文里到处出现，却几乎没人解释。先花十分钟讲清，后面所有表格你都能自己读。

基本功 ①：患病率 ≠ 发病率（横断面只能算前者）

NHANES 是横断面的（第 01 篇第四节讲过）：暴露与结局在同一时点测量，所以你只能算患病率。所以你的说法只有两种：

1. 不能说"发病率"。论文从头到尾用的都是 prevalence，中文就是患病率；
2. 不能说"导致"或"预防"。没有时间先后，就谈不上因果。你只能说"相关"或"同时存在"。

案例里这句话该怎么写："血清维生素 D 水平与抑郁症状患病率呈负相关"（正确），而不是"维生素 D 缺乏导致抑郁风险升高"（错误）。

基本功 ②：95% 置信区间的正确读法

先说它字面上的意思：把同一个研究重复做无数次，每次都算一个区间，其中 95% 的区间会套住真实值。

最常见的误读是"真值有 95% 的概率落在这个区间里"。严格讲这句话不对：真实值是固定的，会变的是区间。

实用读法只有两句话：

1. 区间多宽 = 这个估计有多不确定（样本越小、设计越复杂 → 越宽）；
2. 区间有没有跨过"无效值"：OR 看是否包含 1，均值差看是否包含 0。跨过去了，就说明证据不足以断定有关联。

拿我们的结果读一遍：加权抑郁患病率 11.88%（95%CI 10.03–13.72） → "美国成人抑郁症状患病率估计约 11.9%，考虑到抽样误差，大致在 10%–13.7% 之间"。

基本功 ③：P 值是什么、不是什么

先把它的意思说白：假设真实效应其实为零（统计上叫零假设），那么出现"当前这么大、甚至更大的差异"的概率有多大？ 这个概率就是 P 值。

它不是下面任何一件事：

• ❌ 不是"效应为零的概率"；
• ❌ 不是"结果由偶然造成的概率"；
• ❌ 不是"效应很大的概率"。

两条实用规则：

1. P 值和置信区间是同一个检验的两种表达：对 OR 来说，"P < 0.05"与"95%CI 不含 1"是一回事；
2. P 值不是"有/无"的开关。我们的 Table 1 里，吸烟的组间 P = 0.0572，论文是 0.0475，恰好一个在 0.05 之上、一个在之下。P = 0.057 的正确读法是"证据不足以断定有差异"，而不是"没有差异"；把 0.05 当成一条非黑即白的线，是医学论文里最常见的误用之一。

基本功 ④：odds 与 OR（理解 Table 2 的地基）

这一条最关键，因为第 09 篇整篇都在报 OR。

odds（比值）≠ 概率：

odds = 发生概率 ÷ 不发生概率

抑郁患病率 12% → odds = 0.12 ÷ 0.88 ≈ 0.136。注意 0.136 不是 12%，两者是不同的量。

OR（odds ratio，比值比） = 两组 odds 的比值。例如"暴露每升高 1 ng/mL，抑郁的 odds 变为原来的 0.992 倍"（OR = 0.992）。

为什么论文都用 OR 而不是 RR（风险比）？ 因为 logistic 回归天然给出的是 OR。两者的关系是：

• 结局罕见（<10%）时，OR ≈ RR，可以直接近似解读为"风险"；
• 结局不罕见时，OR 会把效应放大，这时把"odds 降低"说成"风险降低"会夸大结论。

我们的案例：抑郁患病率约 12%（不算罕见），所以最高四分位那行 OR = 0.508 应表述为"odds 降低约 49%"，不能直接写成"风险降低 49%"。（第 10 篇讲结果解读时还会用到这条。）

小结：患病率/发病率决定你能说什么，CI 决定估得准不准，P 值决定证据够不够，odds/OR 决定效应怎么说。四个概念各管一段，混了任何一个，结论都会走偏。

五、第一次加权描述：加权 vs 未加权

先做一次最简单的加权描述，算抑郁患病率：

des <- svydesign(id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTPH2YR, nest=TRUE, data=ana)
svymean(~depression, des)     # 加权患病率
confint(svymean(~depression, des))   # 95% 置信区间

实测结果（3,863 人分析样本，2026 年 9 月）：

图2：（R ggplot 实算图，不是示意图）。同一个人群、三种算法下的患病率几乎重叠——误差线也几乎重叠。

两个结论：

1. 加权后患病率从样本数字变成全国数字：11.88% 才是"美国成年人的抑郁症状患病率"估计；
2. 两种权重结果几乎一致（11.88% vs 11.76%，差 0.12 个百分点）——这印证了"权重敏感性低"。所以案例论文没写明用哪个权重，并不影响结论，我们在报告里也如实记录这一点（第 10 篇）。

还要记住"乘法规则"：如果将来合并多个周期（上一节的 append 用上了），权重不能直接混用，每个周期各有一套权重，合并时要按官方规则换算（例如除以周期数）。跨周期研究务必查官方 Analytic Guidelines 的对应条款。

六、复现论文 Table 1（上）：连续变量

论文 Table 1 是"按有无抑郁分组的加权基线特征表"。连续变量报告加权均值（95%CI），并给组间 P 值。

# 一个变量跑一遍：加权均值 + 95%CI + 组间 P
for (v in c("RIDAGEYR", "INDFMPIR", "BMXBMI", "vitD_ngml", "d2_ngml", "d3_ngml")) {
  d0 <- svymean(as.formula(paste0("~", v)), subset(des, depression == 0))   # 无抑郁组
  d1 <- svymean(as.formula(paste0("~", v)), subset(des, depression == 1))   # 有抑郁组
  p  <- summary(svyglm(as.formula(paste0(v, " ~ depression")), design = des))$coef[2, 4]
  print(rbind(confint(d0), confint(d1))); print(p)
}

说明两点：暴露变量按论文口径换算成 ng/mL（官方单位 nmol/L ÷ 2.5），所以上面的 vitD_ngml 是除以 2.5 之后的列；组间 P 值用加权回归的系数检验（等价于加权 t 检验）。

逐行对照（本表左边是我们的实跑值，方括号里是论文原值）：

6 行全部对上，最大偏差 0.17（维生素 D3：30.86 vs 30.96），多数偏差 ≤0.02。 这个量级的差异来自软件实现与迭代算法，不是数据处理错误，判定复现成功。

另外核对一个数：论文正文报告了分析样本的未加权描述，平均年龄 53.56（SD 16.71）、男性 45.74%。我们实跑：53.56（SD 16.71）、男性 45.74%，逐位一致。这也再次确认：我们和论文切出的是同一批 3,863 人。

七、复现论文 Table 1（下）：分类变量

分类变量报告加权百分比（95%CI），同样给组间 P 值。三个函数分工：

tab  <- svytable(~ male + depression, des)   # 加权列联表
prop <- prop.table(tab, 2) * 100             # 转成每组内的百分比
p    <- svychisq(~ male + depression, des)$p.value   # 组间差异检验

• svytable() 算的是"加权后的人数"，prop.table(..., 2) 把它转成组内百分比；
• 比例的 95%CI 用 svyciprop()（默认 logit 法）；
• svychisq() 做组间比较。注意它默认用 Rao-Scott F 检验（考虑抽样设计的卡方），不是普通 Pearson 卡方。

逐行对照：

百分比全部对上，最大偏差 0.45 个百分点（非西裔黑人：9.75% vs 9.37%）。

一个必须如实交代的差异：P 值不对齐，而且有一处跨过了 0.05

比例对得这么齐，P 值却有肉眼可见的差别（性别 0.0098 vs 0.0034、饮酒 0.0223 vs 0.0048），吸烟更是 0.0572 对 0.0475，一个在 0.05 之上、一个在之下。

这不是算错了，原因是检验方法不同：论文用的软件（Empower）未披露具体卡方实现，而 R 的 svychisq() 默认用 Rao-Scott F 检验（一种专门为复杂抽样设计校正的卡方，带设计自由度的修正）。同一份数据、同一种"卡方检验"的名义下，不同软件的统计量与自由度约定不同，P 值就会有差异——样本量越大、越接近 0.05 边界，这种差异越容易被放大。

所以复现时要分清两类数字：点估计（均值、比例、OR）应当逐位对齐；P 值只要求"同量级、结论方向一致"。如果你的 P 值也算到了 0.05 边缘，别急着改数据——先换一种检验实现（例如试着用设计自由度校正的版本）看看它是不是本来就站在边界上。

八、Table 1 复现验收

里程碑 M4 达成：论文的第一张结果表已经在你手里，而且不是"抄来的"，是你自己算出来的。

【暗线进度条 · 里程碑 M4 ✅】 数据集清洗到论文口径（第 06 篇，N=3,863）→ 权重选好、svydesign 姿势正确、Table 1 加权基线表完整复现 → 下一站：用 svyglm 跑三个递进模型，复现论文 Table 2 与 Figure 2/3（第 08 篇） （对应论文 2.4 权重 + 3.1 结果）

本篇对应的复现脚本段落（完整代码）

总脚本第 6 步就是这张表，代码一字未改。

##### 第 6 步：Table 1——加权基线特征表（按抑郁分组，完整版）#################
#
# 论文 Table 1：把分析样本按"有/无抑郁症状"分成两组，各算一遍加权均值（连续变量）
# 或加权百分比（分类变量），并做组间比较（p 值）。

# 连续变量（加权均值 + 95%CI + 组间比较），含论文 Table 1 的 D2/D3 行
cat("\n===== Table 1 完整版（加权，按抑郁分组）=====\n")
for (v in c("RIDAGEYR", "INDFMPIR", "BMXBMI", "vitD_ngml")) {
  d0 <- svymean(as.formula(paste0("~", v)), subset(des, depression == 0))
  d1 <- svymean(as.formula(paste0("~", v)), subset(des, depression == 1))
  pval <- summary(svyglm(as.formula(paste0(v, " ~ factor(depression)")),
                         design = des))$coef[2, 4]
  cat(sprintf("%-10s 无抑郁 %.2f (%.2f–%.2f) | 有抑郁 %.2f (%.2f–%.2f) | P=%.4f\n",
              v, coef(d0), confint(d0)[1], confint(d0)[2],
              coef(d1), confint(d1)[1], confint(d1)[2], pval))
}
# D2/D3 分量（ng/mL 口径）
for (v in c("LBXVD2MS", "LBXVD3MS")) {
  ana[[paste0(v, "_ngml")]] <- ana[[v]] / 2.5
  des <- update(des, tmp = ana[[paste0(v, "_ngml")]])
  d0 <- svymean(~tmp, subset(des, depression == 0)); d1 <- svymean(~tmp, subset(des, depression == 1))
  pval <- summary(svyglm(tmp ~ factor(depression), design = des))$coef[2, 4]
  cat(sprintf("%-14s 无抑郁 %.2f | 有抑郁 %.2f | P=%.4f\n", paste0(v, "(ng/mL)"), coef(d0), coef(d1), pval))
}

# 分类变量（加权百分比 + 卡方检验），含饮酒四组
for (v in c("male", "race", "education", "smoking", "drinking", "diabetes", "hypertension")) {
  tab <- svytable(as.formula(paste0("~", v, " + depression")), des)
  prop <- prop.table(tab, 2) * 100
  cs <- svychisq(as.formula(paste0("~", v, " + depression")), des)
  pv <- unname(cs$p.value)
  cat(sprintf("%-13s 组间 P=%.4f\n", v, pv))
  for (lv in rownames(prop)) cat(sprintf("   %-18s 无抑郁 %.1f%% | 有抑郁 %.1f%%\n", lv, prop[lv, 1], prop[lv, 2]))
}

本篇常见坑

1. 把患病率写成发病率：横断面数据只能算患病率，也不能说"导致/预防"，这是审稿人一眼就抓的硬伤。
2. 把权重当普通变量乘进数据里：权重要交给 svydesign，让 survey 包在估计和方差里统一处理，不要自己手算加权平均（否则置信区间和 P 值全错）。
3. 先筛子集再建设计对象：方差会算错。正确做法是先 svydesign 全样本，再 subset(des, ...) 对设计对象取子集。
4. nest = TRUE 忘了设：NHANES 的 PSU 编号在层内复用，不设会把不同层的同号 PSU 当同一个，方差偏低。
5. 拿 svytable 的原始输出当百分比：它是加权人数，要 prop.table(..., 2) 转成组内比例。
6. 用普通卡方理解 svychisq 的 P 值：它是抽样设计校正后的检验，P 值与论文不同属正常，只看结论方向与量级。
7. 以为权重能修一切：权重修"抽样设计造成的结构偏差"，修不了"某类人根本没数据"——第 06 篇讲的完整病例偏倚就是权重无解的那部分。
8. 不检查权重为 0 / 缺失的个体：这类人必须从分析样本中排除（本篇 2.1），漏掉会让结果悄悄偏掉且不报错。
9. 把 P 值当"有/无"的开关：P = 0.057 是"证据不足"，不是"没有差异"。
10. 把 OR 当 RR 讲：结局不罕见时 OR 会放大效应，要说"odds 降低"而不是"风险降低"。

下篇预告

第 08 篇：把关联算出来，svyglm() 加权 logistic 回归的完整讲法（系数怎么变成 OR、置信区间怎么来、P 值怎么读），三个递进模型的构建逻辑，连续变量与四分位的双口径结果，以及 Figure 2/3 剂量反应曲线的画法；同时讲清 OR 的量纲含义（"每 1 ng/mL"到底意味着什么），并如实处理论文未披露的四分位切点口径。

参考来源

1. CDC/NCHS. NHANES Tutorials — Weighting. https://wwwn.cdc.gov/nchs/nhanes/tutorials/weighting.aspx
2. CDC/NCHS. NHANES Analytic Guidelines. https://wwwn.cdc.gov/nchs/nhanes/analyticguidelines.aspx
3. Lumley T. survey: analysis of complex survey samples（R 包文档）. https://cran.r-project.org/web/packages/survey/
4. 案例论文：Front Nutr. 2025;12:1545443. PMID: 40497025（Table 1 原值取自 PMC 全文 XML）
5. 本篇加权均值/比例/CI/P 值均为 2026 年 9 月 R 4.5.2 + survey 实算；权重结构实测数据取自 DEMO_L