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=有）

数据清洗要点：

1. 将生存时间从月转换为人年（除以12）
2. 过滤缺失值
3. 检查极端值（比如生存时间超过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. 模型诊断与注意事项

泊松回归的关键假设检查：

1. 等离散性（均值≈方差）：用dispersiontest()检验
2. 异常值检测：检查标准化残差
3. 模型拟合优度：比较残差偏差与自由度

# 过度离散检验 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计数），零膨胀泊松模型可能更合适。