将R2jags对象转换为Stanreg (rstanarm)对象

Question

我用R2jags做了一个模型。我喜欢jags语法，但是我发现R2jags产生的输出并不容易使用。我最近读到了关于rstanarm包的内容。它有许多有用的功能，并且得到了tidybayes和bayesplot包的良好支持，以便轻松地进行模型诊断和可视化。但是，我不喜欢用rstanarm编写模型的语法。理想情况下，我希望在这两个领域中取得最佳效果，即用R2jags编写模型，并将输出转换为Stanreg对象以使用rstanarm函数。

这有可能吗？如果是这样的话，是怎么做的？

Answer 1

我认为问题不一定是它是否可能--我怀疑它可能是可能的。真正的问题是你准备花多少时间来做这件事。您所要做的就是尝试在结构上复制由rstanarm创建的对象，在一定程度上，这对于R2jags输出是可能的。这将使得一些后处理任务可能会起作用。

恕我直言，我怀疑更好地利用您的时间是将R2jags对象转换为您想要使用的后处理函数所使用的对象。例如，只需对JAGS输出进行很小的修改，就可以使用bayesplot中的所有mcmc_*()绘图函数。下面是一个例子。下面是jags()函数帮助中的示例模型。

# An example model file is given in:
model.file <- system.file(package="R2jags", "model", "schools.txt")

# data
J <- 8.0
y <- c(28.4,7.9,-2.8,6.8,-0.6,0.6,18.0,12.2)
sd <- c(14.9,10.2,16.3,11.0,9.4,11.4,10.4,17.6)


jags.data <- list("y","sd","J")
jags.params <- c("mu","sigma","theta")
jags.inits <- function(){
  list("mu"=rnorm(1),"sigma"=runif(1),"theta"=rnorm(J))
}

jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
                n.iter=5000, model.file=model.file, n.chains = 2)

现在，bayesplot的mcmc_*()绘图函数期望的是MCMC绘制的矩阵列表，其中的列名提供了参数的名称。默认情况下，jags()将所有这些元素放入一个矩阵中。在上面的例子中，总共有5000次迭代，其中2500次作为burnin (保留2500次采样)，在这种情况下，n.thin被设置为2 (jags()有一种识别细化参数的算法)，但在任何情况下，jagsfit$BUGSoutput$n.keep元素都会标识保持多少次迭代。在本例中，它是1250。所以你可以用它从输出中生成一个包含两个矩阵的列表。

jflist <- list(jagsfit$BUGSoutput$sims.matrix[1:jagsfit$BUGSoutput$n.keep, ], 
               jagsfit$BUGSoutput$sims.matrix[(jagsfit$BUGSoutput$n.keep+1):(2*jagsfit$BUGSoutput$n.keep), ])

现在，您只需调用一些绘图函数：

mcmc_trace(jflist, regex_pars="theta")

​

​

或

mcmc_areas(jflist, regex_pars="theta")

​

​

因此，与其尝试复制rstanarm生成的所有输出，不如尝试将jags输出转换为适合您想要使用的后处理函数的格式，这样可能会更好地利用您的时间。

编辑-为bayesplot中的pp_check()添加了可能性。

在这种情况下，y的后部提取位于theta参数中。因此，我们创建了一个具有y和yrep元素的对象，并将其设置为foo类

x <- list(y = y, yrep = jagsfit$BUGSoutput$sims.list$theta)
class(x) <- "foo"

然后，我们可以为类foo的对象编写一个pp_check方法。这直接来自于bayesplot::pp_check()的帮助文件。

pp_check.foo <- function(object, ..., type = c("multiple", "overlaid")) {
  y <- object[["y"]]
  yrep <- object[["yrep"]]
  switch(match.arg(type),
         multiple = ppc_hist(y, yrep[1:min(8, nrow(yrep)),, drop = FALSE]),
         overlaid = ppc_dens_overlay(y, yrep[1:min(8, nrow(yrep)),, drop = FALSE]))
}

然后，只需调用该函数：

pp_check(x, type="overlaid")

​

​