r - 如何使用 R 中的 arm 或 rstanarm 包模拟感兴趣的数量？

Question

我想知道如何从使用 R 中的 arm 或 rstanarm 包估计的回归模型中模拟感兴趣的数量。我是新手贝叶斯方法和 R 并且使用 Zelig 包已有一段时间了。我问了一个类似的问题before ，但我想知道是否可以使用这些包估计的后验分布来模拟这些数量。

在 Zelig 中，您可以为独立值设置所需的值，并计算结果变量的结果(预期值、概率等)。一个例子:

# Creating a dataset:
set.seed(10)
x <- rnorm(100,20,10)
z <- rnorm(100,10,5)
e <- rnorm(100,0,1)
y <- 2*x+3*z+e
df <- data.frame(x,z,e,y)

# Loading Zelig
require(Zelig)

# Model
m1.zelig <- zelig(y ~ x + z, model="ls", data=df)
summary(m1.zelig)

# Simulating z = 10
s1 <- setx(m1.zelig, z = 10)
simulation <- sim(m1.zelig, x = s1)
summary(simulation)

因此 Zelig 将 x 保持在其均值 (20.56)，并模拟感兴趣的数量 z = 10。在这种情况下，y 约为 71。

同样的模型使用arm:

# Model
require(arm)
m1.arm <- bayesglm(y ~ x + z, data=df)
summary(m1.arm)

并使用 rstanarm:

# Model
require(rstanarm)
m1.stan <- stanlm(y ~ x + z, data=df)
print(m1.stan)

有什么方法可以模拟 z = 10 且 x 等于其均值，并使用这两个包估计的后验分布得到 y 的期望值？非常感谢!

Answer 1

对于 bayesglm，你可以这样做

sims <- arm::sim(m1.arm, n = 1000)
y_sim <- rnorm(n = 1000, mean = sims@coef %*% t(as.matrix(s1)), sd = sims@sigma)
mean(y_sim)

对于(未发布的)rstanarm，它会是类似的

sims <- as.matrix(m1.stan)
y_sim <- rnorm(n = nrow(sims), mean = sims[,1:(ncol(sims)-1)] %*% t(as.matrix(s1)), 
               sd = sims[,ncol(sims)])
mean(y_sim)

一般来说，对于 Stan，您可以将 s1 作为 row_vector 传递，并在 .stan 文件的生成数量 block 中使用它，例如

generated quantities {
  real y_sim;
  y_sim <- normal_rng(s1 * beta, sigma);
}

在这种情况下，y_sim 的后验分布会在您打印 后验摘要时出现。