负二项混合模型中组特异性离散参数的编码实现问询

我有不同机构的过度离散计数数据，计划用带机构随机效应的负二项随机效应模型建模。经检查，各机构的过度离散程度差异显著：部分机构结果分布接近Poisson，部分高度过度离散，其余介于两者之间。

我尝试用GLMMadaptive::mixed_model通过n_phis和initial_values参数为每个组指定离散参数，先通过MASS::fitdistr估计各机构的离散参数szs再代入，但模型仍仅使用单一离散参数（移除相关参数后模型结果完全一致）。lme4::glmer.nb也无对应功能选项。已知brms的贝叶斯方法可实现该需求，现寻求非贝叶斯R包的编码实现方案，无需讲解统计概念。

附数据生成代码：

set.seed(123)
n <- 100
dta <- data.frame(facility = rep(LETTERS[1:4], each = n),
                  phase = rep(c("Intervention", "Control"), each = n / 2),
                  outcome = c(rnbinom(n, mu = 5, size = 0.02), 
                              rpois(n, lambda = 5),
                              rpois(n, lambda = 8), 
                              rnbinom(n, mu = 8, size = 0.5)))

附之前尝试的代码：

szs <- sapply(unique(dta$facility), \(i) MASS::fitdistr(dta$outcome[dta$facility == i], "negative binomial")[[1]][1])

GLMMadaptive::mixed_model(fixed = outcome ~ phase,
                          random = ~ 1|facility,
                          data = dta,
                          n_phis = length(unique(dta$facility)),
                          initial_values = szs, # 也试过initial_values = c(phis = szs)
                          family = GLMMadaptive::negative.binomial())

方案1：使用glmmTMB包 #

glmmTMB支持通过dispformula参数为分组指定异质离散参数，选择nbinom2负二项模型（方差为mu + mu^2/theta，theta对应负二项的size参数），通过dispformula = ~ 0 + facility让每个机构拥有独立的离散参数：

# 安装加载包
install.packages("glmmTMB")
library(glmmTMB)

# 拟合模型
model_tmb <- glmmTMB(outcome ~ phase + (1|facility),
                     dispformula = ~ 0 + facility,  # 每个facility对应独立离散参数
                     data = dta,
                     family = nbinom2())

# 查看结果
summary(model_tmb)

# 提取各机构的离散参数theta（转换为原始尺度）
exp(fixef(model_tmb)$disp)

注：glmmTMB对离散参数的建模在对数尺度上，因此用exp()转换回原始的theta（即负二项分布的size参数）。

方案2：使用flexmix包 #

flexmix可拟合分组混合模型，每个机构作为一个独立组件，每个组件拥有专属的离散参数，同时支持加入随机效应：

# 安装加载包
install.packages("flexmix")
library(flexmix)

# 拟合模型：按facility分组，每个组对应一个负二项组件，包含随机截距
model_flex <- flexmix(outcome ~ phase + (1|facility) | facility,
                      data = dta,
                      model = FLXMRglmnegbin())

# 查看结果
summary(model_flex)

# 提取各组件（对应facility）的离散参数
sapply(model_flex@components, function(x) x@family$getTheta())

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