Dyl*_*mes 5 r lme4 mixed-models sjplot rstanarm
我正在尝试为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
尽管在parameters0.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模型。知道这里发生了什么吗？
我在一个完全不同的数据集上遇到了同样的问题，但仍然无法弄清楚为什么“替换”比“数据”多 2 行parameters::model_parameters（参见上面的错误）。另一行可能是该reciprocal_dispersion函数不期望的参数，但不确定为什么该函数会有负二项式 glmms 的错误，这很常见。
请注意，包中的tidy_stan函数sjPlot仍然适用于这些模型，但会发出警告：
Warning message:
'tidy_stan' is deprecated.
Use 'parameters::model_parameters()' instead.
See help("Deprecated")
然而，parameters::model_parameters()如上所述，还不起作用。
小智 4
尽管这是一个很大的挑战，但我终于找到了这个错误（并且有一个简单的修复方法，如果太长而无法阅读，请转到帖子的末尾）。我通过查找导致错误的指令来汇集线程。从...开始：
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 包“insight”）。