用R语言dlnm包分析空气污染滞后效应：手把手教你复现芝加哥死亡数据案例-程序员充电站

用R语言dlnm包解析空气污染滞后效应：从数据清洗到结果可视化的完整指南

芝加哥的冬日总是灰蒙蒙的。作为一名公共卫生研究员，我常常盯着窗外被雾霾笼罩的天际线思考：这些悬浮颗粒究竟会在未来多少天内持续影响居民健康？直到发现了dlnm包，这个困扰我多年的问题终于有了科学的分析工具。本文将带您完整复现一个空气污染滞后效应分析案例，从数据导入到结果解读，手把手教您掌握这个强大的时间序列分析利器。

1. 环境准备与数据导入

在开始之前，我们需要确保所有必要的R包已经安装。dlnm包虽然是核心，但还需要几个辅助包来支持我们的分析：

# 安装必要包（如果尚未安装） install.packages(c("dlnm", "splines", "tsModel", "ggplot2")) # 加载包 library(dlnm) # 分布滞后非线性模型 library(splines) # 样条函数支持 library(tsModel) # 时间序列建模 library(ggplot2) # 高级绘图

芝加哥空气污染数据集包含1987-2000年间的每日记录，我们可以使用以下代码导入：

# 读取数据（假设文件位于工作目录） chicago_data <- read.csv("chicago_air.csv", header = TRUE, stringsAsFactors = FALSE) # 查看数据结构 str(chicago_data)

典型的环境健康数据集应包含以下关键变量：

变量名	类型	描述	单位
date	Date	记录日期	YYYY-MM-DD
death	int	每日死亡人数	人
pm10	num	PM10浓度	µg/m³
o3	num	臭氧浓度	ppb
temperature	num	日均温度	°F
humidity	num	相对湿度	%
dow	factor	星期几（因子变量）	1-7

提示：在实际分析前，务必检查数据的完整性和异常值。使用summary()和hist()函数快速了解数据分布特征。

2. 构建交叉基矩阵：模型的核心引擎

交叉基矩阵是dlnm包的核心创新，它同时捕捉暴露-反应关系和滞后效应。理解其参数设置对正确建模至关重要。

2.1 PM10的线性滞后模型

假设我们想研究PM10对死亡率的滞后影响，考虑15天的滞后窗口：

# 构建PM10的交叉基矩阵 cb_pm10 <- crossbasis( x = chicago_data$pm10, lag = 15, # 15天滞后窗口 argvar = list(fun = "lin"), # 线性暴露-反应关系 arglag = list(fun = "poly", degree = 4) # 4阶多项式滞后函数 )

参数解释：

argvar：定义暴露-反应关系
- fun="lin"表示线性关系
- 也可选择"ns"(自然样条)、"bs"(B样条)等非线性关系
arglag：定义滞后效应分布
- degree=4表示使用4阶多项式
- 也可选择"strata"(分层)等其他函数

2.2 温度的非线性控制模型

温度通常需要作为混杂因素控制，我们采用非线性建模：

# 构建温度的交叉基矩阵 cb_temp <- crossbasis( x = chicago_data$temperature, lag = 3, argvar = list(fun = "ns", df = 5, cen = 65), # 以65°F为参考 arglag = list(fun = "strata", breaks = 1) # 滞后分层 )

这里使用了几个关键技巧：

cen=65：设定65°F(约18°C)为参考温度
breaks=1：将0-1天和1-3天的滞后效应分开建模

3. 拟合统计模型与结果预测

有了交叉基矩阵后，我们可以构建广义线性模型：

# 拟合准泊松回归模型 model <- glm( death ~ cb_pm10 + cb_temp + ns(time, df = 7*14) + dow, family = quasipoisson(), data = chicago_data ) # 模型摘要 summary(model)

注意：ns(time, df=7*14)用于控制长期趋势和季节效应，其中df参数建议设为研究周期年数×14。

3.1 生成预测结果

使用crosspred()函数计算特定暴露水平的风险：

# 预测PM10增加0-20µg/m³的影响 pred_pm10 <- crosspred( cb_pm10, model, at = 0:20, # PM10浓度范围 bylag = 0.2, # 滞后步长 cumul = TRUE # 计算累积效应 )

关键参数：

at：指定暴露水平的计算点
bylag：控制滞后曲线的平滑度
cumul：是否计算累积效应

4. 结果可视化与专业解读

4.1 单日暴露的滞后效应曲线

# 绘制PM10增加10µg/m³的滞后效应 plot(pred_pm10, "slices", var = 10, col = "darkgreen", ylab = "相对风险(RR)", main = "PM10增加10µg/m³的日滞后效应", ci.arg = list(col = "lightgreen"))

曲线解读要点：

横轴：暴露后的天数(滞后天数)
纵轴：相对风险(RR)
初期高峰：通常出现在暴露后2-3天
后期下降：可能由于补偿效应或测量误差
置信区间：反映估计的不确定性

4.2 累积效应分析

# 绘制累积效应 plot(pred_pm10, "slices", var = 10, cumul = TRUE, col = "red", ylab = "累积相对风险", main = "PM10增加10µg/m³的累积效应")

累积效应图揭示了更完整的公共卫生影响：

短期累积：前3-5天的风险叠加
长期趋势：可能出现的风险逆转现象
关键值：峰值累积RR及其滞后天数

4.3 结果表格输出

对于学术报告，我们常需要数值结果：

# 提取关键结果 results <- data.frame( Lag = seq(0, 15, by = 1), RR = pred_pm10$matRRfit[11, seq(1, 151, by = 10)], # 10µg/m³的RR LowCI = pred_pm10$matRRlow[11, seq(1, 151, by = 10)], HighCI = pred_pm10$matRRhigh[11, seq(1, 151, by = 10)] ) # 打印表格 knitr::kable(results, caption = "PM10增加10µg/m³的滞后效应估计")

示例输出表格：

Lag	RR	LowCI	HighCI
0	1.012	1.008	1.016
1	1.018	1.013	1.023
...	...	...	...
15	0.998	0.994	1.002

5. 高级技巧与常见问题排查

5.1 模型诊断与敏感性分析

优秀的分析需要验证模型稳健性：

# 敏感性分析：改变滞后窗口 cb_pm10_sens <- crossbasis(chicago_data$pm10, lag = c(10, 20), # 测试不同窗口 argvar = list(fun = "lin"), arglag = list(fun = "poly", degree = 3)) # 比较模型拟合 AIC(update(model, . ~ . - cb_pm10 + cb_pm10_sens))

常见敏感性分析维度：

滞后窗口长度(如7天vs15天)
自由度选择(使用AIC/BIC准则)
参考值设定(特别是对非线性关系)

5.2 常见报错与解决方案

在实际操作中，您可能会遇到：

数据格式错误
- 症状：Error in crossbasis(): x must be a numeric vector
- 解决：检查变量类型class(chicago_data$pm10)
内存不足
- 症状：cannot allocate vector of size...
- 解决：减少at参数的分辨率或简化模型
收敛问题
- 症状：glm()算法未收敛
- 解决：调整glm.control(maxit=100)

5.3 使用ggplot2增强可视化

基础图形有时不能满足出版要求，我们可以用ggplot2增强：

library(ggplot2) # 准备绘图数据 plot_data <- data.frame( Lag = rep(seq(0, 15, length.out = 50), 2), RR = c(pred_pm10$matRRfit[11, ], pred_pm10$cumRRfit[11, ]), Type = rep(c("单日效应", "累积效应"), each = 50) ) # 创建高级图形 ggplot(plot_data, aes(x = Lag, y = RR, color = Type)) + geom_line(size = 1.2) + geom_hline(yintercept = 1, linetype = "dashed") + labs(title = "PM10滞后效应分析", x = "滞后天数", y = "相对风险(RR)") + theme_minimal()