为什么`parameters::model_parameters`会为负二项式`rstanarm`模型抛出错误答案

【问题标题】：Why `parameters::model_parameters` throws an error for negative binomial `rstanarm` model为什么`parameters::model_parameters`会为负二项式`rstanarm`模型抛出错误
【发布时间】：2020-10-13 19:29:45
【问题描述】：

我正在尝试为stan_glmer.nb (rstanarm) 输出创建一个表，但是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 得到一个空白输出，而不是错误（参见 cmets）。无论哪种方式，这个函数似乎都不适用于rstanarm 离散分布（参见 cmets）。据我在帮助文档中看到的，没有任何迹象表明需要指定负二项式模型。此外，该功能适用于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() 还不起作用。

【问题讨论】：

性能 v 0.6.1 我得到 Stan/NB 模型的空白/空输出，而不是错误...
这个：dat<-data.frame(x=rnorm(500),z=rep(c("A","B","C","D","E"),100), y=.2+x*.7) 对我不起作用.. y 表达式中的“x”从何而来？
@BenBolker 我不知道我是怎么不小心写了performance 而不是parameters，但它看起来不像是performance 被调用了——很抱歉造成混乱。为清楚起见，对问题进行了编辑。
@MrSmithGoesToWashington 抱歉，我已将 x=rnorm(500) 拉到数据框参数之外，现在应该可以使用了。当我创建问题时，我必须已经运行了这条线并且在我的工作环境中有 x。
也许需要考虑： neg bin 是离散的 .. 尝试使用另一个离散分布我得到同样的错误 ..

标签： r lme4 mixed-models sjplot rstanarm

【解决方案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（因此出现错误）。但是，失败的模块与先验有关。我们可以通过简单地设置先验等于 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”）。

【讨论】：