R语言实战：基于泊松回归模型计算人年发病率及其置信区间
1. 泊松回归模型基础
泊松回归是广义线性模型的一种，专门用于分析计数数据或事件发生率。在流行病学和临床研究中，我们经常需要计算人年发病率（incidence rate per person-year），比如每1000人年的疾病发生率。泊松回归的核心假设是：事件发生的次数服从泊松分布，且其方差等于均值。
举个生活中的例子：假设你经营一家奶茶店，想预测每天顾客投诉的次数。投诉次数是典型的计数数据，泊松回归可以帮你分析影响因素（比如天气、店员数量）。在医学场景中，我们可能关注的是“每1000人年随访中癌症复发的次数”。
泊松回归的数学形式为：
log(λ) = β₀ + β₁X₁ + ... + βₖXₖ其中λ是事件发生率，X是自变量。这个对数链接函数（log link）确保预测值始终为正数。
2. 数据准备与清洗
我们使用乳腺癌生存数据作为示例（数据可通过公开渠道获取）。首先导入并检查数据：
library(dplyr) bc <- read.csv("Breast_cancer_survival_agec.csv", header=TRUE) head(bc) # 重命名乱码列 names(bc)[1] <- "id" # 关键变量说明： # status: 死亡状态（0=存活，1=死亡） # time: 生存时间（月） # ln_yesno: 淋巴结肿大（0=无，1=有）数据清洗要点：
- 将生存时间从月转换为人年（除以12）
- 过滤缺失值
- 检查极端值（比如生存时间超过10年的记录）
bc <- bc %>% mutate(pyears = time/12) %>% filter(!is.na(status) & !is.na(pyears))3. 人年发病率计算原理
人年发病率的计算公式为：
发病率 = 事件发生数 / 总人年数 × 1000以淋巴结肿大分组为例：
# 无淋巴结肿大组 group0 <- bc %>% filter(ln_yesno==0) events0 <- sum(group0$status) pyears0 <- sum(group0$pyears) rate0 <- events0/pyears0*1000 # 有淋巴结肿大组 group1 <- bc %>% filter(ln_yesno==1) events1 <- sum(group1$status) pyears1 <- sum(group1$pyears) rate1 <- events1/pyears1*1000手工计算结果可能如下：
- 无淋巴结肿大组：11例/1000人年
- 有淋巴结肿大组：28例/1000人年
4. 使用glm构建泊松回归模型
R中的glm()函数可以方便地拟合泊松模型：
model <- glm(status ~ ln_yesno + offset(log(pyears)), family=poisson, data=bc) summary(model)关键参数说明：
offset(log(pyears))：将人年数作为偏移项，相当于在方程右边加上log(人年)family=poisson：指定泊松分布
模型输出会包含：
- 系数估计（Estimate）
- 标准误（Std. Error）
- z值（z value）
- p值（Pr(>|z|)）
5. 计算发病率比与置信区间
泊松回归的系数经过指数变换后就是发病率比（Rate Ratio, RR）：
# 提取系数 coef <- summary(model)$coefficients rr <- exp(coef[2,1]) # 第二行是ln_yesno的系数 # 计算95% CI se <- coef[2,2] lower <- exp(coef[2,1] - 1.96*se) upper <- exp(coef[2,1] + 1.96*se)假设输出显示RR=0.40（95%CI:0.25-0.64），意味着：
- 有淋巴结肿大的患者死亡率是无淋巴结肿大患者的0.4倍
- 或者说，无淋巴结肿大组的死亡率是有淋巴结肿大组的2.5倍（1/0.4）
6. 使用popEpi包简化计算
对于不想手动计算的研究者，popEpi包提供了便捷函数：
library(popEpi) rate_ratio(x = c(events0, pyears0), y = c(events1, pyears1), SE.method = FALSE)这个函数会直接输出：
- 两组的人年发病率
- 发病率比（RR）
- 95%置信区间
7. 结果可视化
用ggplot2展示结果更直观：
library(ggplot2) result <- data.frame( Group = c("No Lymph Node", "Lymph Node"), Rate = c(rate0, rate1), Lower = c(rate0_lower, rate1_lower), Upper = c(rate0_upper, rate1_upper) ) ggplot(result, aes(x=Group, y=Rate)) + geom_col(width=0.5, fill="#4E84C4") + geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.2) + labs(title="Incidence Rate per 1000 Person-Years", y="Rate (95% CI)") + theme_minimal()8. 模型诊断与注意事项
泊松回归的关键假设检查：
- 等离散性（均值≈方差）：用
dispersiontest()检验 - 异常值检测：检查标准化残差
- 模型拟合优度：比较残差偏差与自由度
# 过度离散检验 library(AER) dispersiontest(model) # 残差分析 res <- residuals(model, type="pearson") plot(res, ylab="Pearson Residuals")如果存在过度离散（overdispersion），考虑改用负二项回归：
library(MASS) model_nb <- glm.nb(status ~ ln_yesno + offset(log(pyears)), data=bc)9. 扩展应用
9.1 调整混杂因素
可以加入年龄、肿瘤分级等协变量：
model_adj <- glm(status ~ ln_yesno + age + histgrad + offset(log(pyears)), family=poisson, data=bc)9.2 分层分析
按年龄分层计算发病率：
bc %>% group_by(age_group, ln_yesno) %>% summarise(events=sum(status), pyears=sum(pyears), rate=events/pyears*1000)9.3 时间依赖性分析
如果需要考虑发病率随时间变化，可以加入时间交互项：
bc$time_cat <- cut(bc$pyears, breaks=c(0,2,5,10)) model_time <- glm(status ~ ln_yesno*time_cat + offset(log(pyears)), family=poisson, data=bc)10. 常见问题解决
问题1：R报错“NA/NaN/Inf in foreign function call”
- 检查数据是否有缺失值
- 确保人年数没有0值
问题2：置信区间过宽
- 可能是事件数太少，考虑合并分类
- 检查是否有异常值影响
问题3：模型不收敛
- 尝试增加最大迭代次数：
glm(..., control=list(maxit=100))- 检查自变量间是否存在完全分离
实际分析中，我发现当样本量较小时，泊松回归的置信区间可能不够准确。这时可以考虑使用精确泊松检验或贝叶斯方法。另外，对于零膨胀数据（很多0计数），零膨胀泊松模型可能更合适。