我正在尝试在bayesplot对象上通过inla包实现函数，并且有点不确定如何从后验预测分布中得出。我想我差不多了，但是rstan的抽奖比inla的抽奖更具可变性。
在rstan中，使用bayesplot vignette的简化示例，我可以：
library(bayesplot)
library(ggplot2)
library(rstanarm)
library(ggpubr)
library(tidyverse)
#rstan model set up
roaches$roach100 <- roaches$roach1 / 100 # pre-treatment number of roaches (in 100s)
fit_poisson <- stan_glm(y ~ roach100 + treatment + senior,offset = log(exposure2),family = poisson(link = "log"),data = roaches,seed = 1111,refresh = 0)
#In order to use the PPC functions from the bayesplot package we need a vector y of outcome values:
y <- roaches$y
#and a matrix yrep of draws from the posterior predictive distribution,yrep_poisson <- posterior_predict(fit_poisson,draws = 500)
#then plot:
p1 <- bayesplot::ppc_dens_overlay(y,yrep_poisson[1:50,])
p1
我想将该图复制到inla对象上。根据{{1}}小插图，您可以执行此操作，因为它们提供了定义简单的bayesplot方法的代码，该方法创建了例如类的拟合模型对象。 pp_check：
foo要使用pp_check.foo <- function(object,type = c("multiple","overlaid"),...) {
type <- match.arg(type)
y <- object[["y"]]
yrep <- object[["yrep"]]
stopifnot(nrow(yrep) >= 50)
samp <- sample(nrow(yrep),size = ifelse(type == "overlaid",50,5))
yrep <- yrep[samp,]
if (type == "overlaid") {
ppc_dens_overlay(y,yrep,...)
} else {
ppc_hist(y,...)
}
}
，我们只需列出包含pp_check.foo和y组件的列表，并将其赋予foo类：
yrepx <- list(y = rnorm(200),yrep = matrix(rnorm(1e5),nrow = 500,ncol = 200))
class(x) <- "foo"
#create plot above:
pp_check(x,type = "overlaid")
#create same model but in inla:
library(inla)
fit_poisson_inla <- inla(y ~ roach100 + treatment + senior,control.predictor = list(compute = T),family = "poisson")
返回每个inla_object_name$marginals.fitted.values的后验预测分布：
y我认为从中反复采样将满足我的需要，但是实际上创建每个分布时，每个观察值只能使用75个值（fit_poisson_inla$marginals.fitted.values
#so to get distribution for first oberservation:
fitted.Predictor.1 <- fit_poisson_inla$marginals.fitted.values[[1]]
），而实际上我希望从所有范围的值中采样。我认为我们可以通过使用线性预测变量使用dim(fitted.Predictor.1)来做到这一点（第4.3 here节）：
inla.tmarginal 
我的问题是，如何从这个fitted_dist <- fit_poisson_inla$marginals.linear.predictor
#should i have used "inla.rmarginal(n,marginal)"?
marginal_dist <- lapply(fitted_dist,function(y) inla.tmarginal(function(x) {exp(x)},y)) %>% map(~ as.data.frame(.) %>% rename(.,xx = x))
#resample 500 times
yrep_poisson_inla <- as.matrix(bind_rows(rerun(500,lapply(marginal_dist,function(x) sample(x$xx,1)) %>% as.data.frame())))
#convert to class foo for pp_check
x <- list(y = y,yrep = yrep_poisson_inla[1:50,])
class(x) <- "foo"
p2 <- pp_check(x,type = "overlaid")
#plot
ggarrange(p1,p2,ncol = 1,nrow = 2,labels = c("rstan","inla sample"))
（inla）对象的后验预测分布中正确得出平局矩阵并传递到fit_poisson_inla？ pp_check产生离散值，而yrep_poisson产生连续值。 yrep_poisson_inla抽奖中的变化远远大于rstan（第二个图）。我做的是否正确，这只是一些抽样问题，还是不同方法的产物？在更复杂的示例中，差异可能很大。
谢谢