在贝叶斯模型中绘制交互效果（使用 rstanarm）答案

【问题标题】：Plotting interaction effects in Bayesian models (using rstanarm)在贝叶斯模型中绘制交互效果（使用 rstanarm）
【发布时间】：2017-11-04 17:18:56
【问题描述】：

我试图在 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?

【问题讨论】：

标签： r bayesian stan rstan

【解决方案1】：

在这种（线性，高斯，恒等链接，无截距）情况下，

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

所以，要绘制x 或z 的边际效应，您只需要一个适当的范围和系数的后验分布，您可以通过

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)

请注意，在这种情况下，列名只是观察结果。

【讨论】：

抱歉，我不清楚我在这里看到了什么。看起来对 dmu_dx 的分配正在计算后验的每次绘制的边际效应。这很好，但这是一个平均边际效应。如果我想指定“将 z 从 a 增加到 b 对 y 的影响的影响”，那么解决方案不是更像 dmu_dx1
对dmu_dx（或dmu_dz）的赋值为dat中的每个点创建了一个边际效应，并将相同的东西赋值给colnames(dmu_dx)将这个值放在@的垂直轴上987654336@。这对于真正连续预测变量的边际效应（瞬时变化率）有效。如果您想要虚拟变量从 0 到 1 或更一般的 a 和 b 变化的后验分布，您可以按照您的建议进行操作。
我跟着这个，它做了一个这样的情节teblunthuis.cc/outgoing/Goodrich_Marginal_Effects_Plot.png我有太多的数据。应该绘制一个样本。

【解决方案2】：

您也可以使用ggeffects-package，尤其是对于边际效应；或 sjPlot-package 用于边际效应和其他绘图类型（对于边际效应，sjPlot 只是简单地包装了 ggeffects 中的函数）。

要绘制交互的边际效应，请使用 sjPlot::plot_model() 和 type = "int"。使用mdrt.values 定义为连续调节变量绘制哪些值，并使用ppd 让预测基于线性预测变量的后验分布或后验预测分布。

library(sjPlot)
plot_model(m1, type = "int", terms = c("x", "z"), mdrt.values = "meansd")

plot_model(m1, type = "int", terms = c("x", "z"), mdrt.values = "meansd", ppd = TRUE)

或绘制其他特定值的边际效应，使用type = "pred" 并在terms-参数中指定值：

plot_model(m1, type = "pred", terms = c("x", "z [10, 20, 30, 40]"))
# same as:
library(ggeffects)
dat <- ggpredict(m1, terms = c("x", "z [10, 20, 30, 40]"))
plot(dat)

有更多选项，以及自定义情节外观的不同方式。请参阅相关的帮助文件和包小插曲。

【讨论】：