r - 绘制贝叶斯模型中的交互作用(使用 rstanarm)

标签 r bayesian stan rstan

我试图在 rstanarm() 的贝叶斯线性模型中展示一个变量的影响如何随另一个变量的值而变化。我能够拟合模型并从后验中抽取以查看每个参数的估计值,但尚不清楚如何给出一个变量在相互作用中的影响作为其他变化和相关不确定性的某种图(即边际效应图)。以下是我的尝试:

library(rstanarm)

# Set Seed
set.seed(1)

# Generate fake data
w1 <- rbeta(n = 50, shape1 = 2, shape2 = 1.5)
w2 <- rbeta(n = 50, shape1 = 3, shape2 = 2.5)

dat <- data.frame(y = log(w1 / (1-w1)),
                  x = log(w2 / (1-w2)),
                  z = seq(1:50))

# Fit linear regression without an intercept:
m1 <- rstanarm::stan_glm(y ~ 0 + x*z, 
                         data = dat,
                         family = gaussian(),
                         algorithm = "sampling",
                         chains = 4,
                         seed = 123,
                         )


# Create data sets with low values and high values of one of the predictors
dat_lowx <- dat
dat_lowx$x <- 0

dat_highx <- dat
dat_highx$x <- 5

out_low <- rstanarm::posterior_predict(object = m1,
                                   newdata = dat_lowx)

out_high <- rstanarm::posterior_predict(object = m1,
                                        newdata = dat_highx)

# Calculate differences in posterior predictions
mfx <- out_high - out_low

# Somehow get the coefficients for the other predictor?

最佳答案

在这种(线性、高斯、恒等链接、无截距)情况下,

mu = beta_x * x + beta_z * z + beta_xz * x * z 
   = (beta_x + beta_xz * z) * x 
   = (beta_z + beta_xz * x) * z

因此,要绘制 xz 的边际效应,您只需要适当的范围和系数的后验分布,您可以通过以下方式获得

post <- as.data.frame(m1)

然后

dmu_dx <- post[ , 1] + post[ , 3] %*% t(sort(dat$z))
dmu_dz <- post[ , 2] + post[ , 3] %*% t(sort(dat$x))

然后您可以使用类似下面的方法来估计数据中每个观察的单个边际效应,它计算了 x 对每个观察的 mu 的影响在您的数据中以及 z 对每次观察的 mu 的影响。

colnames(dmu_dx) <- round(sort(dat$x), digits = 1)
colnames(dmu_dz) <- dat$z
bayesplot::mcmc_intervals(dmu_dz)
bayesplot::mcmc_intervals(dmu_dx)

请注意,在这种情况下,列名只是观察值。

关于r - 绘制贝叶斯模型中的交互作用(使用 rstanarm),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/47113400/

相关文章:

linux - 使用 R 在 Linux 中的 Stan 程序中找不到 prep_call_sampler

r - R 中的数值积分 : bad integrand behaviour error

java - 在 RHEL 6 上安装 R R-java-devel 需要 java-devel

r - on.exit 示例,添加参数

r - 如何从 bsts R 包中提取包含概率

python - 如何降低计算成本

c# - 使用 Accord.Net 的编码对象编码第二个数据集

c++ - 从 C++ 程序调用 Stan 例程

stan - 如何从 stanfit 对象中提取估计值

r - gbm::interact.gbm 与 dismo::gbm.interactions