为什么 `parameters::model_parameters` 会为负二项式 `rstanarm` 模型抛出错误
Why `parameters::model_parameters` throws an error for negative binomial `rstanarm` model
我正在尝试为 stan_glmer.nb (rstanarm) 输出创建一个 table，但是包 parameters 中的 model_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 参数，但不确定为什么函数会有负二项式 glmms 的错误，这很常见。
请注意，sjPlot 包中的 tidy_stan 函数仍然适用于这些模型，但会给出警告：
Warning message:
'tidy_stan' is deprecated.
Use 'parameters::model_parameters()' instead.
See help("Deprecated")
然而，如上所述，parameters::model_parameters() 还没有生效。
虽然这是一个很大的挑战，但我终于弄清楚了这个错误（并且有一个简单的修复方法，如果阅读时间太长，请转到 post 的末尾）。我通过查找导致错误的指令来合并线程。
开始于：
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（因此出现错误）。但是，失败的模块与先验有关。我们可以通过简单地将先验设置为 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-package“insight”）。
我正在尝试为 stan_glmer.nb (rstanarm) 输出创建一个 table，但是包 parameters 中的 model_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 参数，但不确定为什么函数会有负二项式 glmms 的错误，这很常见。
请注意，sjPlot 包中的 tidy_stan 函数仍然适用于这些模型，但会给出警告：
Warning message:
'tidy_stan' is deprecated.
Use 'parameters::model_parameters()' instead.
See help("Deprecated")
然而，如上所述，parameters::model_parameters() 还没有生效。
虽然这是一个很大的挑战，但我终于弄清楚了这个错误（并且有一个简单的修复方法，如果阅读时间太长，请转到 post 的末尾）。我通过查找导致错误的指令来合并线程。 开始于：
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（因此出现错误）。但是，失败的模块与先验有关。我们可以通过简单地将先验设置为 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-package“insight”）。