基于Poisson回归估算性别调整发病率的技术咨询

是的，完全可以用R的glm()函数完成这个需求——Poisson回归属于广义线性模型的范畴，glm()通过指定family = poisson()就能实现。下面是完整的步骤和代码：

1. 导入数据集 #

先把你提供的数据转换成R可识别的数据框：

# 构造数据集
dat <- data.frame(
  Group = factor(c("Disease", "Disease", "Comparators", "Comparators")),
  Sex = factor(c("women", "men", "women", "men")),
  Number = c(721, 309, 3894, 1614),
  Events = c(95, 39, 289, 109),
  PY = c(4502.817, 1124.824, 25169.310, 6421.807)
)

2. 拟合性别调整的Poisson回归模型 #

由于我们要估算的是发病率（事件数/人年），需要将人年（PY）以对数形式作为偏移量（offset）加入模型，这样模型会直接输出率的估计值。模型中纳入Group（分组）和Sex（性别），控制性别后计算不同组的发病率：

# 拟合Poisson回归模型
model <- glm(Events ~ Group + Sex + offset(log(PY)), 
             data = dat, 
             family = poisson())

3. 计算性别调整后的发病率及置信区间 #

这里推荐用emmeans包来计算边际均值（即调整性别后的平均发病率），它能直接给出对数尺度的估计值，我们再通过指数转换得到原始尺度的发病率：

# 如果没安装emmeans包先安装
# install.packages("emmeans")
library(emmeans)

# 计算Group的边际均值（调整性别），并转换为原始尺度
adj_rates <- emmeans(model, ~ Group, type = "response")

# 查看结果
adj_rates

运行上述代码后，会得到类似如下的输出（数值为示例）：

Group       response       SE  df lower.CL upper.CL
 Disease     0.0205    0.0021  2   0.0167   0.0251
 Comparators 0.0115    0.0007  2   0.0102   0.0130

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

4. 整理成你需要的结果格式 #

从输出中提取数值，即可得到：

• 患病群体的性别调整发病率为每PY 0.0205例（95% CI: 0.0167-0.0251）
• 对照群体的性别调整发病率为每PY 0.0115例（95% CI: 0.0102-0.0130）

补充说明 #

• 偏移量offset(log(PY))的作用：Poisson模型默认预测事件计数，加入对数人年作为偏移量后，模型会预测事件率（事件数/人年），这正是我们需要的发病率。
• emmeans的作用：自动计算控制性别后的边际发病率，相当于对性别进行了标准化（基于数据中的性别分布），得到的是组内调整了性别混杂后的平均发病率。

内容的提问来源于stack exchange，提问作者Mabai