基于独立负二项分布构建Dirichlet多项分布的死亡预测问题

问题背景 #

现有数据集包含ISO3国家代码、2010-2019年年份、24类死亡原因及对应死亡数，每年各死因死亡数总和固定，目标是预测2020-2023年各国各死因的预期死亡数。数据存在过度离散，原计划用R语言mgcv包的gam函数为每个国家-死因组合拟合独立负二项回归模型，但预测结果的总和无法匹配固定值；尝试手动缩放预测值会破坏模型拟合效果。参考资料提到“Dirichlet多项分布等价于具有相同尺度参数的独立负二项分布在总和条件下的集合”，需验证该思路并寻求可行实现方案。

样本数据 #

单国4类死因的样本数据如下：

# A tibble: 40 × 3
    Year Cause                      Deaths
   <dbl> <chr>                       <dbl>
 1  2010 Communicable diseaes       97203.
 2  2010 Non-communicable diseases   2743.
 3  2010 Injuries                    7128.
 4  2010 Other                      12712.
 5  2011 Communicable diseaes       96536.
 6  2011 Non-communicable diseases   2716.
 7  2011 Injuries                    7264.
 8  2011 Other                      13088.
 9  2012 Communicable diseaes       95786.
10  2012 Non-communicable diseases   2642.
11  2012 Injuries                    7140.
12  2012 Other                      13556.
13  2013 Communicable diseaes       94153.
14  2013 Non-communicable diseases   2652.
15  2013 Injuries                    7191.
16  2013 Other                      14058.
17  2014 Communicable diseaes       92192.
18  2014 Non-communicable diseases   2828.
19  2014 Injuries                    7288.
20  2014 Other                      14825.
21  2015 Communicable diseaes       96254.
22  2015 Non-communicable diseases   2617.
23  2015 Injuries                    8128.
24  2015 Other                      16339.
25  2016 Communicable diseaes      101976.
26  2016 Non-communicable diseases   2470.
27  2016 Injuries                    8366.
28  2016 Other                      17007.
29  2017 Communicable diseaes       97219 
30  2017 Non-communicable diseases   2365 
31  2017 Injuries                    7675 
32  2017 Other                      15686 
33  2018 Communicable diseaes       95826 
34  2018 Non-communicable diseases   2180 
35  2018 Injuries                    7937 
36  2018 Other                      15807 
37  2019 Communicable diseaes       97264 
38  2019 Non-communicable diseases   2080 
39  2019 Injuries                    8031 
40  2019 Other                      15660  

现有模型代码 #

原计划使用的独立负二项GAM模型代码：

library(mgcv)
library(dplyr)

# original dataset
df <- df %>%
  mutate(expected_cause = NA,
         expected_cause_se = NA) %>%
  arrange(Year)

# create a dataset to stored predictions for 2020-2023
pred_pand <- expand.grid(Year = rep(c(2020, 2021,2022,2023), each = 1),
                               Cause = rep(unique(df$Cause))) %>%
  mutate(expected_cause = NA,
         expected_cause_se = NA) %>%
  arrange(Year)

for (i in unique(df$Cause)) {
  # create temp dataset
  whichs <- which(df$Cause == i)
  temp <- df[whichs, ]

  # run the model
  model <- gam(Deaths ~ s(Year), data = temp, family = nb(theta = NULL, link = "log"))
  pred <- predict(model, se.fit = TRUE, type = "response", newdata = 
  data.frame(Year = c(2020,2021,2022,2023)))

  # store predictions for 2020-2023
  whichs_pred <- which(pred_pand$Cause == i)
  pred_pand[whichs_pred, "expected_cause"] <- pred$fit
  pred_pand[whichs_pred, "expected_cause_se"] <- pred$se.fit
  
  # store predictions for years prior to 2019
  years <- unique(temp$Year)
  range_years <- min(years):max(years)
  pred_hist <- predict(model, se.fit = TRUE, type = "response",
                       newdata = data.frame(Year = range_years))
  df[whichs, "expected_cause"] <- pred_hist$fit
  df[whichs, "expected_cause_se"] <- pred_hist$se.fit

}

解决方案思路与实现 #

1. 验证Dirichlet-负二项关联思路 #

该思路是正确的：当多个独立负二项分布共享相同的尺度参数θ时，在总和固定的条件下，各变量的联合分布等价于Dirichlet多项分布。核心是让所有死因的模型共享过度离散参数，并通过条件约束保证总和匹配。

2. 三种可行实现方案 #

方案A：后缩放校正（简单易行）

基于现有独立GAM模型的预测结果，按年份将预测值缩放至固定总死亡数，同时校正标准误以保留不确定性：

# 假设已知2020-2023年的固定总死亡数（替换为实际值）
total_deaths <- c("2020" = 120000, "2021" = 122000, "2022" = 121000, "2023" = 123000)

# 按年份缩放预测值
pred_pand <- pred_pand %>%
  group_by(Year) %>%
  mutate(
    pred_total = sum(expected_cause),
    scale_factor = total_deaths[as.character(Year)] / pred_total,
    # 缩放预测均值
    expected_cause_scaled = expected_cause * scale_factor,
    # 缩放标准误（近似校正）
    expected_cause_se_scaled = expected_cause_se * scale_factor
  ) %>%
  ungroup()

注意：该方法不会扭曲各死因的比例趋势，仅调整量级，若原模型的比例预测合理，是高效的解决方案。

方案B：共享θ的联合贝叶斯建模（贴合理论）

使用brms包构建联合负二项模型，让所有死因共享过度离散参数，再通过后处理缩放匹配固定总和：

library(brms)
library(tidyr)

# 转换数据为宽格式（每行对应一个年份，各死因列为独立变量）
df_wide <- df %>%
  pivot_wider(names_from = Cause, values_from = Deaths)

# 拟合共享θ的联合负二项GAM模型
model <- brm(
  bf(mvbind(Communicable_diseaes, Non_communicable_diseases, Injuries, Other) ~ s(Year)),
  data = df_wide,
  family = negbinomial(link = "log"),
  # 设定共享shape参数（即θ的倒数）
  prior = prior(constant(1), class = shape),
  chains = 4, cores = 4
)

# 生成2020-2023年的预测
pred_new <- tibble(Year = 2020:2023)
pred_fit <- predict(model, newdata = pred_new, type = "response")

# 缩放预测值至固定总死亡数
pred_scaled <- t(apply(pred_fit, 1, function(x) {
  year_total <- total_deaths[as.character(pred_new$Year[which(apply(pred_fit,1,identical,x))])]
  x / sum(x) * year_total
}))

# 整理为长格式结果
pred_pand_brms <- expand.grid(Year=2020:2023, Cause=unique(df$Cause)) %>%
  arrange(Year, Cause) %>%
  mutate(expected_cause = as.vector(pred_scaled))

方案C：直接拟合Dirichlet多项模型（原生支持比例约束）

使用DirichletReg包直接建模各死因的比例，自然满足总和固定的约束，同时处理过度离散：

library(DirichletReg)
library(tidyr)

# 整理数据为DirichletReg要求的格式：每行对应一个年份，包含各死因计数和总数
df_prop <- df %>%
  group_by(Year) %>%
  summarise(
    total = sum(Deaths),
    counts = list(Deaths)
  ) %>%
  ungroup() %>%
  mutate(
    # 转换为矩阵格式的响应变量
    y = lapply(counts, function(x) matrix(x, nrow=1))
  )

# 拟合带Year平滑项的Dirichlet多项模型
model_dir <- DirichletReg(y ~ s(Year), data = df_prop, family = "multinomial")

# 预测2020-2023年的死因比例
pred_prop <- predict(model_dir, newdata = data.frame(Year=2020:2023), type="response")

# 转换为死亡数（乘以固定总死亡数）
pred_pand_dir <- expand.grid(Year=2020:2023, Cause=unique(df$Cause)) %>%
  arrange(Year, Cause) %>%
  mutate(
    expected_cause = as.vector(pred_prop) * rep(total_deaths, each=length(unique(df$Cause)))
  )

关键注意事项 #

• 若预测年份的总死亡数未知，需先单独拟合总死亡数的负二项GAM模型，得到预测总数后再进行比例分配。
• 共享过度离散参数是保证Dirichlet-负二项理论一致性的核心，避免独立模型导致的参数矛盾。
• 后缩放方法需注意校正标准误，否则会低估/高估预测的不确定性。

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