如何用R语言surveil包计算研究全时段的平均发病率？

全时段平均发病率的核心计算逻辑是全时段总发病数除以总暴露人口数，再换算为每10万人口的发病数。结合surveil的贝叶斯模型，我们可以基于后验样本计算该指标的分布，得到更稳健的点估计和置信区间，也可以直接用观测数据计算粗平均发病率。

方法一：基于贝叶斯模型后验的平均发病率 #

利用stan_rw输出的后验发病数（lambda），按分组汇总全时段的发病数，再结合总人口数计算：

1. 先计算每组全时段的总人口数：

total_pop <- tx.msa %>%
  group_by(Race) %>%
  summarise(total_pop = sum(Population), .groups = "drop")

2. 提取模型后验发病数并汇总全时段总发病数：

# lambda是三维数组：迭代数 × 时间点 × 分组
lambda_post <- fit$lambda
# 按迭代和分组汇总总发病数
total_cases_post <- apply(lambda_post, c(1, 3), sum)

3. 计算每个后验样本的平均发病率（/10万），并生成摘要：

# 计算每组的平均发病率后验分布
mean_inc_post <- t(apply(total_cases_post, 1, function(x) {
  x / total_pop$total_pop * 100000
}))
colnames(mean_inc_post) <- total_pop$Race

# 计算点估计（均值）和95%置信区间
mean_inc_summary <- data.frame(
  Race = total_pop$Race,
  mean_inc = apply(mean_inc_post, 2, mean),
  lower_95ci = apply(mean_inc_post, 2, quantile, 0.025),
  upper_95ci = apply(mean_inc_post, 2, quantile, 0.975)
)

# 输出格式化结果
kable(mean_inc_summary, digits = 2)

方法二：基于观测数据的粗平均发病率 #

如果不需要贝叶斯推断，直接用原始汇总数据计算：

crude_mean_inc <- tx.msa %>%
  group_by(Race) %>%
  summarise(
    mean_inc = sum(Count) / sum(Population) * 100000
  ) %>%
  select(Race, mean_inc)

kable(crude_mean_inc, digits = 2)

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