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, headerTRUE) 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_yesno0) events0 - sum(group0$status) pyears0 - sum(group0$pyears) rate0 - events0/pyears0*1000 # 有淋巴结肿大组 group1 - bc %% filter(ln_yesno1) 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)), familypoisson, databc) summary(model)关键参数说明offset(log(pyears))将人年数作为偏移项相当于在方程右边加上log(人年)familypoisson指定泊松分布模型输出会包含系数估计Estimate标准误Std. Errorz值z valuep值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)假设输出显示RR0.4095%CI:0.25-0.64意味着有淋巴结肿大的患者死亡率是无淋巴结肿大患者的0.4倍或者说无淋巴结肿大组的死亡率是有淋巴结肿大组的2.5倍1/0.46. 使用popEpi包简化计算对于不想手动计算的研究者popEpi包提供了便捷函数library(popEpi) rate_ratio(x c(events0, pyears0), y c(events1, pyears1), SE.method FALSE)这个函数会直接输出两组的人年发病率发病率比RR95%置信区间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(xGroup, yRate)) geom_col(width0.5, fill#4E84C4) geom_errorbar(aes(yminLower, ymaxUpper), width0.2) labs(titleIncidence Rate per 1000 Person-Years, yRate (95% CI)) theme_minimal()8. 模型诊断与注意事项泊松回归的关键假设检查等离散性均值≈方差用dispersiontest()检验异常值检测检查标准化残差模型拟合优度比较残差偏差与自由度# 过度离散检验 library(AER) dispersiontest(model) # 残差分析 res - residuals(model, typepearson) plot(res, ylabPearson Residuals)如果存在过度离散overdispersion考虑改用负二项回归library(MASS) model_nb - glm.nb(status ~ ln_yesno offset(log(pyears)), databc)9. 扩展应用9.1 调整混杂因素可以加入年龄、肿瘤分级等协变量model_adj - glm(status ~ ln_yesno age histgrad offset(log(pyears)), familypoisson, databc)9.2 分层分析按年龄分层计算发病率bc %% group_by(age_group, ln_yesno) %% summarise(eventssum(status), pyearssum(pyears), rateevents/pyears*1000)9.3 时间依赖性分析如果需要考虑发病率随时间变化可以加入时间交互项bc$time_cat - cut(bc$pyears, breaksc(0,2,5,10)) model_time - glm(status ~ ln_yesno*time_cat offset(log(pyears)), familypoisson, databc)10. 常见问题解决问题1R报错“NA/NaN/Inf in foreign function call”检查数据是否有缺失值确保人年数没有0值问题2置信区间过宽可能是事件数太少考虑合并分类检查是否有异常值影响问题3模型不收敛尝试增加最大迭代次数glm(..., controllist(maxit100))检查自变量间是否存在完全分离实际分析中我发现当样本量较小时泊松回归的置信区间可能不够准确。这时可以考虑使用精确泊松检验或贝叶斯方法。另外对于零膨胀数据很多0计数零膨胀泊松模型可能更合适。