r语言 - 为什么 'parameters::model_parameters' 会为负二项式 'rstanarm' 模型抛

Question

我正在尝试为stan_glmer.nb(rstanarm)输出创建一个表，但是model_parameters包parameters抛出一个奇怪的错误，我不确定如何解决。也许这是一个错误。

缩短了版本信息的sessionInfo()输出：

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
parameters_0.8.2
rstanarm_2.21.1

一个可重现的示例：

library(rstanarm)
library(parameters)
x<-rnorm(500)
dat<-data.frame(x=x,z=rep(c("A","B","C","D","E"),100), y=.2+x*.7)
mod1<-stan_glmer(y~x+(x|z),data=dat)
model_parameters(mod1, effects="all")

我将在这里省略输出，因为它并不重要，但函数有效。 现在负二项式模型：

dat.nb<-data.frame(x=rnorm(500),z=rep(c("A","B","C","D","E"),100),
y=rnbinom(500,size=1,prob = .5))
mod2<-stan_glmer.nb(y~x+(x|z),data=dat.nb)
model_parameters(mod2, effects="all")

现在出现一条错误消息：

Error in `$<-.data.frame`(`*tmp*`, "parameter", value = c("(Intercept)",  : 
replacement has 3 rows, data has 1

尽管使用parameters版本0.10.1，@BenBolker得到一个空白输出，而不是错误(请参阅注释)。无论哪种方式，此函数似乎都不适用于rstanarm离散分布(请参阅注释)。据我在帮助文档中看到的，没有任何内容表明需要指定负二项式模型。此外，该函数适用于lme4模型：

library(lme4)
mod1<-lmer(y~x+(x|z),data=dat)
model_parameters(mod1, effects="all")
mod2<-glmer.nb(y~x+(x|z),data=dat.nb)
model_parameters(mod2, effects="all")

对于这些模拟数据，存在一些模型收敛问题等，但model_parameters适用于glmer.nb模型，但不适用于stan_glmer.nb模型。知道这里发生了什么吗？

---

我在完全不同的数据集上遇到了同样的问题，但仍然无法弄清楚为什么"替换"在parameters::model_parameters中比"数据"多 2 行(请参阅上面的错误)。另外一行可能是函数不期望的reciprocal_dispersion参数，但不确定为什么函数会出现负二项式 glmm 的错误，这很常见。

请注意，软件包中的tidy_stan函数仍然适用于sjPlot这些模型，但给出了警告：

Warning message:
'tidy_stan' is deprecated.
Use 'parameters::model_parameters()' instead.
See help("Deprecated") 

然而，如上所述，parameters::model_parameters()尚未奏效。

Answer 1

虽然这是一个相当大的挑战，但我终于弄清楚了这个错误(并且有一个简单的修复，如果太长而无法阅读，请转到帖子的末尾)。我通过查找导致错误的指令来池化线程。 从以下开始：

model_parameters(model = mod2, effects="all")

在以下情况下失败：

parameters:::model_parameters.stanreg(model = mod2, effects="all", prior = T)

在以下情况下失败：

params <- parameters:::.extract_parameters_bayesian(mod2, centrality = "median", 
dispersion = F, ci = 0.89, ci_method = "hdi", 
test = "pd", rope_range = "default", rope_ci = 1, 
bf_prior = NULL, diagnostic = "ESS", priors = T, 
effects = "fixed", standardize = NULL)

在以下情况下失败：

parameters <- bayestestR:::describe_posterior.stanreg(mod2, centrality = "median", 
dispersion = F, ci = 0.89, ci_method = "hdi", 
test = "pd", rope_range = "default", rope_ci = 1, 
bf_prior = NULL, diagnostic = "ESS", priors = T)

在以下情况下失败：

priors_data <- bayestestR:::describe_prior.stanreg(mod2)

在以下情况下失败：

priors <- insight:::get_priors.stanreg(mod2)

为了找出它从这里开始失败的阶段，我复制了这个函数的源代码(现在定义为GET_priors)，并放置了一些战略打印：

GET_priors <- function(x) # Modified with tags to see where it fails
{
ps <- rstanarm::prior_summary(mod2)
l <- insight:::.compact_list(lapply(ps[c("prior_intercept", "prior")], 
function(.x) {
if (!is.null(.x)) {
if (is.na(.x$dist)) {
.x$dist <- "uniform"
.x$location <- 0
.x$scale <- 0
.x$adjusted_scale <- 0
}
.x <- do.call(cbind, .x)
as.data.frame(.x)
}
}))

print("STEP1")

cn <- unique(unlist(lapply(l, colnames)))
l <- lapply(l, function(.x) {
missing <- setdiff(cn, colnames(.x))
if (length(missing)) {
.x[missing] <- NA
}
.x
})

print("STEP2")

if(length(l) > 1)
{
prior_info <- do.call(rbind, l)
}
else
{
cn <- colnames(l[[1]])
prior_info <- as.data.frame(l)
colnames(prior_info) <- cn
}

print("STEP3")

flat <- which(prior_info$dist == "uniform")
if (length(flat) > 0) {
prior_info$location[flat] <- NA
prior_info$scale[flat] <- NA
prior_info$adjusted_scale[flat] <- NA
}

print("STEP4")
print(prior_info)
print(insight:::find_parameters(x)$conditional)

prior_info$parameter <- insight:::find_parameters(x)$conditional

print("STEP4.1")

prior_info <- prior_info[, intersect(c("parameter", 
"dist", "location", "scale", "adjusted_scale"), 
colnames(prior_info))]

print("STEP4.2")

colnames(prior_info) <- gsub("dist", "distribution", 
colnames(prior_info))

print("STEP4.3")

colnames(prior_info) <- gsub("df", "DoF", colnames(prior_info))

print("STEP4.4")

priors <- as.data.frame(lapply(prior_info, function(x) {
if (insight:::.is_numeric_character(x)) {
as.numeric(as.character(x))
}
else {
as.character(x)
}
}), stringsAsFactors = FALSE)

print("STEP5")

string <- strsplit(names(priors), "_", fixed = TRUE)
string <- lapply(string, insight:::.capitalize)
names(priors) <- unlist(lapply(string, paste0, collapse = "_"))
priors
}
GET_priors(mod2)
# [1] "STEP1"
# [1] "STEP2"
# [1] "STEP3"
# [1] "STEP4"
#                   dist location scale   adjusted_scale
# prior_intercept normal        0   2.5             <NA>
# prior           normal        0   2.5 2.63656782500616
# [1] "(Intercept)"           "x"                     "reciprocal_dispersion"
#  Error in `$<-.data.frame`(`*tmp*`, "parameter", value = c("(Intercept)",  : 
#  replacement has 3 rows, data has 2 

出于某种奇怪的原因，它试图将 3 行的列添加到包含 2 行的 data.frame 中(因此出现了错误)。但是，失败的模块与先验有关。我们可以通过简单地将 prior 设置为等于 F 来获得结果，避免代码中的所有分支，如下所示：

model_parameters(model = mod2, effects="all", prior = F)
# Fixed effects 
# 
# Parameter             |   Median |            CI |     pd | % in ROPE |  Rhat |  ESS
# ------------------------------------------------------------------------------------
#   (Intercept)           | 8.05e-03 | [-0.11, 0.13] | 54.00% |    81.15% | 1.002 | 1738
# x                     |    -0.12 | [-0.25, 0.00] | 94.67% |    37.18% | 1.000 | 2784
# reciprocal_dispersion |     0.97 | [ 0.75, 1.22] |   100% |        0% | 1.000 | 4463
# 
# # Random effects SD/Cor: z
# 
# Parameter       |    Median |            CI |     pd | % in ROPE |  Rhat |  ESS
# -------------------------------------------------------------------------------
#   (Intercept)     |  3.43e-03 | [ 0.00, 0.03] |   100% |    98.30% | 1.002 | 2077
# x ~ (Intercept) | -9.39e-09 | [-0.01, 0.01] | 50.05% |    99.75% | 1.001 | 2099
# x               |  2.93e-03 | [ 0.00, 0.02] |   100% |    99.08% | 1.001 | 2664
# 
# Using highest density intervals as credible intervals.

事实上，这是一个错误，应该报告(只是错误在依赖项的依赖项中;例如R包"洞察力")。