0

我使用R2jags. 我喜欢这种jags语法，但我发现生成的输出R2jags不容易使用。我最近读到了这个rstanarm包。它具有许多有用的功能，并得到tidybayes和bayesplot软件包的良好支持，便于模型诊断和可视化。但是，我不喜欢用rstanarm. 理想情况下，我想两全其美，即编写模型R2jags并将输出转换为Stanreg对象以使用rstanarm函数。

那可能吗？如果是这样，怎么做？

4

1 回答 1

1

我认为那么问题不一定是它是否可能 - 我怀疑它可能是。问题实际上是你准备花多少时间去做这件事。您所要做的就是尝试在结构中复制由 生成的对象rstanarm，以使其R2jags输出可能。这将使一些后处理任务可能会起作用。

如果我可能如此大胆，我怀疑更好地利用您的时间是将R2jags对象变成可以与您想要使用的后处理功能一起使用的东西。例如，只需对 JAGS 输出进行少量修改，即可使所有mcmc_*()绘图功能bayesplot正常工作。这是一个例子。下面是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)

现在，expect 的mcmc_*()绘图函数bayesplot是 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输出转换为适合您想要使用的后处理功能的格式。

编辑- 增加了pp_check()from 的可能性bayesplot。

在这种情况下，后部绘制y在theta参数中。所以，我们制作一个有元素的对象，y并yrep使它成为类foo

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

然后我们可以为 class 的对象编写一个pp_check方法foo。这直接来自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")

在此处输入图像描述

于 2020-09-09T13:58:38.593 回答