如何用modelsummary并排展示mice填补的glmmTMB条件与零膨胀模型？

我用以下R代码完成缺失值填补和零膨胀混合模型拟合：

library(glmmTMB)
library(mice)

Owls <- Owls

Owls[sample(1:nrow(Owls),10), "ArrivalTime"] <- NA

predictorMatrix <- quickpred(
  Owls,
  include = colnames(Owls),
  exclude = c(),
  mincor = 0.1
)

imp <- mice(
  data = Owls,
  m = 5,
  maxit = 10,
  method = "pmm",
  predictorMatrix = predictorMatrix,
  seed = 123
)

fit <- with(imp, 
             glmmTMB(SiblingNegotiation~(FoodTreatment+ArrivalTime)*SexParent+
                       offset(log(BroodSize))+(1|Nest),
                     data=data.frame(mget(ls())),
                     ziformula=~FoodTreatment+ArrivalTime,
                     family=poisson))

我希望通过modelsummary把模型的条件部分和零膨胀部分并排展示，但试了Stack Overflow上的方案都无效。排查后发现broom.mixed::tidy和broom.mixed::glance无法为mice合并后的模型生成component列。

我自己写了一段代码，但不确定是否正确，同时有个疑问：零膨胀模型的glance统计量在哪里可以获取？

get_modsum_list <- function(model, effects = "fixed", comp = "cond"){
  tidy_df <- broom.mixed::tidy(pool(model), effects = effects, component = comp)
  tidy_df <- tidy_df %>% mutate(component=pool(model)$pooled$component)
  tidy_df <- tidy_df %>% filter(component==!!comp)
  glance_df <- pool(model)$glanced %>% summarise(across(everything(), mean))
  mod <- list(
    tidy = tidy_df,
    glance = glance_df
  )
  class(mod) <- "modelsummary_list"
  mod
}

list(
  'conditional' = get_modsum_list(fit),
  'zero-inflated' = get_modsum_list(fit, comp="zi")
) |>
  modelsummary(stars=TRUE)

一、自定义代码的正确性与优化 #

你的代码逻辑是可行的，但有可以简化和更稳妥的地方：

1. broom.mixed::tidy(pool(model), component = comp)已经会自动过滤对应component的系数，不需要再手动filter
2. 直接用传入的comp参数赋值component列，避免从pool(model)$pooled提取可能出现的长度不匹配问题

优化后的函数：

library(dplyr)
library(modelsummary)

get_modsum_list <- function(model, effects = "fixed", comp = "cond"){
  # 提取对应component的系数表并补充component标记
  tidy_df <- broom.mixed::tidy(pool(model), effects = effects, component = comp) %>%
    mutate(component = comp)
  
  # 对每个填补模型的统计量取均值
  glance_df <- pool(model)$glanced %>%
    summarise(across(everything(), mean, na.rm = TRUE))
  
  # 包装为modelsummary兼容格式
  mod <- list(tidy = tidy_df, glance = glance_df)
  class(mod) <- "modelsummary_list"
  mod
}

二、零膨胀部分的glance统计量获取 #

pool(model)$glanced包含的是每个填补模型的整体统计量（如AIC、BIC），如果需要单独提取零膨胀部分的统计量，需要遍历每个填补后的模型，分别提取对应component的信息再合并：

# 提取每个填补模型的零膨胀部分统计量
zi_glance_list <- lapply(fit$analyses, function(x) {
  broom.mixed::glance(x, component = "zi")
})

# 合并所有填补模型的结果并取均值
zi_glance_df <- bind_rows(zi_glance_list) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

如果要让modelsummary分别展示条件和零膨胀部分的专属统计量，可以修改自定义函数，让它返回对应component的glance结果：

get_modsum_list <- function(model, effects = "fixed", comp = "cond"){
  tidy_df <- broom.mixed::tidy(pool(model), effects = effects, component = comp) %>%
    mutate(component = comp)
  
  # 提取对应component的专属统计量
  glance_list <- lapply(model$analyses, function(x) {
    broom.mixed::glance(x, component = comp)
  })
  glance_df <- bind_rows(glance_list) %>%
    summarise(across(everything(), mean, na.rm = TRUE))
  
  mod <- list(tidy = tidy_df, glance = glance_df)
  class(mod) <- "modelsummary_list"
  mod
}

三、最终调用代码 #

用优化后的函数调用，即可得到条件与零膨胀部分并排的结果：

list(
  "条件模型" = get_modsum_list(fit, comp = "cond"),
  "零膨胀模型" = get_modsum_list(fit, comp = "zi")
) |>
  modelsummary(
    stars = TRUE,
    caption = "缺失值填补后的零膨胀泊松混合模型",
    coef_map = function(x) gsub("\\(Intercept\\)", "截距", x) # 可选：将英文系数名转为中文
  )

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