贝叶斯模型中交互效应的绘制(使用rstanarm)

Question

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

Answer 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)

请注意，列名只是本例中的观察结果。

Answer 2

您也可以使用ggeffects效果-包，特别是边缘效果；或者sjPlot-包用于边缘效果和其他地块类型(对于边际效果，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-argument中的值来绘制其他特定值的边际效果：

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)

​

​

有更多的选项，也有不同的方式定制的情节外观。请参阅相关的帮助文件和软件包小片段。