如何从INLA对象正确抽取后验预测分布样本适配bayesplot？

我最近在尝试给INLA对象使用bayesplot包的函数，但卡在了如何正确从后验预测分布里抽取样本。自己做了个初步实现，却发现rstan生成的样本比INLA的变异性大很多，想请教下具体的解决方案。

rstan实现示例（参考bayesplot文档） #

这是用rstanarm完成的标准流程，用来做对比参考：

library(bayesplot)
library(ggplot2)
library(rstanarm)
library(ggpubr)
library(tidyverse)

# 预处理数据：将蟑螂数量转换为百只为单位
roaches$roach100 <- roaches$roach1 / 100
# 搭建泊松模型
fit_poisson <- stan_glm(
  y ~ roach100 + treatment + senior,
  offset = log(exposure2),
  family = poisson(link = "log"),
  data = roaches,
  seed = 1111,
  refresh = 0
)

# 提取观测值
y <- roaches$y
# 从后验预测分布抽取500个样本
yrep_poisson <- posterior_predict(fit_poisson, draws = 500)
# 绘制PPC密度叠加图
p1 <- bayesplot::ppc_dens_overlay(y, yrep_poisson[1:50, ])
p1

适配INLA的初步尝试 #

根据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, yrep, ...)
  }
}

# 测试这个自定义方法
x <- list(y = rnorm(200), yrep = matrix(rnorm(1e5), nrow = 500, ncol = 200))
class(x) <- "foo"
pp_check(x, type = "overlaid")

INLA模型搭建与样本抽取的尝试 #

接下来搭建了INLA的泊松模型，然后尝试抽取样本：

library(INLA)
fit_poisson_inla <- inla(
  y ~ roach100 + treatment + senior,
  offset = log(exposure2),
  data = roaches,
  control.predictor = list(compute = T),
  family = "poisson"
)

# 查看每个观测的后验预测分布（线性预测器的边际分布）
fit_poisson_inla$marginals.linear.predictor
# 以第一个观测的分布为例
fitted.Predictor.1 <- fit_poisson_inla$marginals.linear.predictor[[1]]

因为每个观测的分布只有75个点，我尝试先把线性预测器转换为指数尺度（对应泊松分布的均值），再抽取样本：

# 将线性预测器的边际分布转换为指数尺度（泊松均值）
marginal_dist <- lapply(fit_poisson_inla$marginals.linear.predictor, function(marg) {
  inla.tmarginal(function(x) exp(x), marg) %>% 
    as.data.frame() %>% 
    rename(xx = x)
})

# 抽取500次样本，每次给每个观测抽一个值
yrep_poisson_inla <- as.matrix(
  bind_rows(
    rerun(500, lapply(marginal_dist, function(d) sample(d$xx, 1)) %>% as.data.frame())
  )
)

# 使用自定义的pp_check方法绘图
x <- list(y = y, yrep = yrep_poisson_inla[1:50, ])
class(x) <- "foo"
p2 <- pp_check(x, type = "overlaid")

# 对比rstan和INLA的结果
ggarrange(p1, p2, ncol = 1, nrow = 2, labels = c("rstan", "inla sample"))

核心疑问 #

• 问题1：如何从INLA对象（比如fit_poisson_inla）中正确抽取可传入pp_check的后验预测分布样本矩阵？
• 问题2：当前结果里，yrep_poisson（rstan生成）是离散值，而yrep_poisson_inla是连续值，且rstan样本的变异性远大于INLA。这是我抽样方式有误，还是两种方法本身的差异？如果是复杂模型场景，这种差异会扩大吗？

内容的提问来源于stack exchange，提问作者user63230