我有（过度分散的）计数数据，这些数据是在不同的设施收集的，我想使用负二项随机效应模型进行建模，其中设施具有随机效应。在进一步检查数据时，每个设施的过度分散变化很大（即，一些设施中的结果分布近似于泊松分布，其他设施高度过度分散，一些介于两者之间），即：

我想指定模型中每个组的离差。我尝试使用GLMMadaptive::mixed_model，因为它有参数n_phis和initial_values，文档指定允许您输入多少个参数及其初始值，所以我通过MASS:fitdistr估计了色散参数（szs）并尝试将其插入：

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, # also tried initial_values = c(phis = szs) given the function documentation
                          family = GLMMadaptive::negative.binomial())

然而，这是不成功的，因为它仍然使用单个分散参数建模（删除n_phis和initial_values命令返回完全相同的模型）。我在lme4::glmer.nb中没有看到这种类型的任何选项。
经过一些搜索，我看到这可以通过brms包中的贝叶斯方法来完成（交叉验证问题在这里和here），但是可以使用使用其他方法的包（即GLMMadaptive::mixed_model，lme4::glmer.nb或其他包）来完成吗？需要说明的是，我的问题是关于编码实现，而不是底层的统计概念。
数据

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)))