3

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

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 的期望值？非常感谢！

4

1 回答 1

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 作为 a 传递row_vector并在 .stan 文件的生成数量块中使用它，例如

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

在这种情况下，当您进行后验总结y_sim时会出现后验分布。print

于 2014-06-25T18:30:12.700 回答