从零开始复现一篇6.2分CHARLS纵向研究：烹饪燃料与呼吸健康的关联分析
1. 为什么选择CHARLS数据做纵向研究
CHARLS（中国健康与养老追踪调查）作为国内最具代表性的中老年健康数据库之一，其价值在于长达十年的追踪设计。我最初接触这个数据集时，最惊讶的是它完整记录了2011-2018年间同一批受访者的健康指标、生活习惯和环境暴露数据。这种"时间维度"正是纵向研究的黄金矿脉——比如我们可以观察同一个体从使用煤炉到改用天然气的过程中肺功能变化，这种证据强度远超单次横断面调查。
具体到烹饪燃料研究，CHARLS的三大优势特别突出：
- 暴露评估全面：不仅记录主要燃料类型（煤气/电力/生物质等），还包含使用频率、厨房通风等细节
- 结局指标丰富：肺功能检测、呼吸道症状问卷、医疗诊断记录形成多维度验证
- 混杂因素详尽：吸烟史、职业暴露、居住环境等关键协变量均有覆盖
不过真实操作时会遇到两个典型挑战：一是变量分散在数十个数据文件中需要拼接，二是同一变量在不同调查波次的命名规则可能变化。有次我为了找"厨房通风情况"这个变量，花了三天时间比对2013和2015年的问卷差异。
2. 数据获取与清洗实战
2.1 数据下载与初步探索
从CHARLS官网下载数据时，建议直接获取STATA格式的"Harmonized CHARLS"版本，这个经过预处理的版本已经统一了变量名和编码规则。我通常按这个顺序下载文件：
- 核心问卷（Demographics、Health Status）
- 生物测量（Physical Measurements）
- 家庭环境（Household Environment）
- 社区数据（Community Data）
初次清洗时必做的四步操作：
// 示例：处理2011基线数据 use "wave1.dta", clear keep ID householdID wave // 保留关键ID和时间标记 merge 1:1 ID using "health_wave1.dta" // 合并健康数据 drop if _merge != 3 // 剔除匹配失败样本 label values COPD yesno // 统一编码标签2.2 构建纵向数据结构
这是最考验耐心的环节。我的经验是先用Excel绘制"变量-波次"映射表，标出每个目标变量在各年份的原始名称。比如肺功能指标FEV1在：
- 2011年：pef
- 2015年：lung_test1
- 2018年：resp_exam3
推荐使用STATA的reshape命令转换数据结构：
// 将宽数据转为长数据 reshape long @fuel @lung @smoke, i(ID) j(wave 2011 2013 2015 2018)这个步骤常会遇到的问题包括：时间点缺失导致样本量骤减、变量单位不一致（比如2013年前燃料使用量记录为"小时/天"，之后改为"是/否"）。有个取巧的办法是查阅原文献的补充材料，通常作者会披露处理细节。
3. 统计建模的关键选择
3.1 模型选型思路
当看到"纵向"二字时，新手常犯的错误是直接套用混合效应模型。实际上需要先回答三个问题：
- 研究问题本质：是观察暴露的累积效应？还是追踪个体变化轨迹？
- 数据结构特征：时间点是否等距？缺失模式是否随机？
- 结局变量类型：连续型（肺功能值）还是二分类（ COPD诊断）？
原文献使用的广义相加混合模型（GAMM）是个灵活选择，既能处理非线性关系（燃料使用年数与肺功能下降的曲线形态），又考虑了个体随机效应。R语言实现示例：
library(mgcv) model <- gam(FEV1 ~ s(fuel_years) + s(age) + smoke + s(ID, bs="re"), data=lung_data)3.2 结果解读要点
纵向分析最大的价值在于揭示"变化中的关联"。比如原文发现：
- 短期（3年内）改用清洁燃料对已有COPD患者症状改善不明显
- 长期（5年以上）坚持使用天然气可显著减缓肺功能年下降率
这类发现必须配合敏感性分析才有说服力。我通常会做：
- 不同缺失数据处理方法的比较（完整病例分析 vs 多重填补）
- 关键参数的不同设定（如平滑函数自由度）
- 亚组分析（城乡差异、性别差异）
4. 从复现到创新的跨越
完全复现文献只是起点。基于现有数据可以做的拓展分析包括：
- 暴露累积效应：构建燃料暴露指数（FEI），结合使用年限和类型加权
- 临界点分析：使用阈值模型识别肺功能急剧下降的暴露阈值
- 联合暴露影响：空气污染与烹饪燃料的交互作用分析
有个实用的建议：在Results部分先复现原文核心表格，然后在Discussion部分加入自己的拓展分析。这样既展示了对前人工作的尊重，又体现了学术增量。比如我曾在复现某篇CHARLS研究后，增加了"燃料转换时间窗"分析，最终帮助论文影响因子提升了1.5分。
最后提醒几个容易踩的坑：
- CHARLS的抽样权重必须纳入分析，否则可能高估城市人群效应
- 不同波次的肺功能检测协议可能有细微差异，需要查阅技术文档
- 广义相加模型的平滑曲线解释需要配合效应值表格，避免视觉误导