使用 rstanarm 计算二项式 logit 中的边际效应答案

【问题标题】：Calculating marginal effects in binomial logit using rstanarm使用 rstanarm 计算二项式 logit 中的边际效应
【发布时间】：2017-07-11 14:30:51
【问题描述】：

根据这篇文章，我正在尝试获得边际效应：http://andrewgelman.com/2016/01/14/rstanarm-and-more/

td <- readRDS("some data")

CHAINS <- 1
CORES <- 1
SEED <- 42
ITERATIONS <- 2000
MAX_TREEDEPTH <- 9

md <- td[,.(y,x1,x2)] # selection the columns i need. y is binary


glm1 <- stan_glm(y~x1+x2,
                 data = md,
                 family = binomial(link="logit"),
                 prior = NULL,
                 prior_intercept = NULL,
                 chains = CHAINS,
                 cores = CORES,
                 seed = SEED,
                 iter = ITERATIONS,
                 control=list(max_treedepth=MAX_TREEDEPTH)
)

# launch_shinystan(glm1) 


tmp <- posterior_predict(glm1,newdata=md[,.(x1,x2)])

问题

运行此代码后，我收到以下错误：我收到一个错误，y not found，这实际上意味着我还需要在newdata 中传递y，根据?posterior_predict 不应该是这种情况

推理

我需要tmp <- posterior_predict(glm1,newdata=md[,.(x1,x2)])，因为根据上面的帖子（据我所知），为了计算 x1 的边际效应（如果我假设 x1 是二进制的）将是

temp <- md
temp[,x1:=0]
temp[,x2:=mean(x2)]
number_0 <- posterior_predict(glm1,newdata=temp)

temp <- md
temp[,x1:=1]
temp[,x2:=mean(x2)]
number_1 <- posterior_predict(glm1,newdata=temp)

marginal_effect_x1 <- number_1 - number_0

【问题讨论】：

虽然与您的问题无关，但仅使用 1 条链并不是一个好主意。而且您不必将max_treedepth 从其默认值（rstanarm 中的 15 与 rstan 中的 10）减少；将其设置为较高的值不会在未达到时造成任何伤害。

标签： r stan rstan rstanarm

【解决方案1】：

对于二元 logit 模型，连续变量的边际效应是关于该变量的成功概率的导数，根据链式法则，它是逻辑密度（在预测变量的某些值下评估，通常是预测变量的观测值）乘以相关变量的系数。在你的情况下，那将是 df <- as.data.frame(glm1) ME <- df$x2 * dlogis(posterior_linpred(glm1)) 由于这取决于预测变量的观察值，因此通常对数据进行平均 AME <- rowMeans(ME) 在二进制预测器的情况下，您可以通过 x1 = 1 时的成功概率减去 x1 = 0 时的成功概率 nd <- md nd$x1 <- 0 p0 <- posterior_linpred(glm1, newdata = nd, transform = TRUE) nd$x1 <- 1 p1 <- posterior_linpred(glm1, newdata = nd, transform = TRUE) ME <- p1 - p0 AME <- rowMeans(ME)

【讨论】：

当我运行二进制代码时，我得到一个向量而不是一个标量。但边际效应只是一个标量，对吧？
AME 是一个向量，因为它是一个后验分布。从概念上讲，平均边际效应是一个标量，但你不知道它是什么标量，所以你有一个关于它是什么的信念分布。这与贝叶斯估计框架中的任何其他数量相同。
另外，predict 函数旨在优化后使用，而不是 MCMC 或 ADVI。 predict 函数使用系数的边际中位数进行预测，它忽略了系数估计值与其周围不确定性之间的相关性。因此，它可能与 posterior_linpred 和 posterior_predict 完全不同。
在 logit 情况下，posterior_predict 产生 0 和 1。对于大量模拟，0 和 1 的平均值将与posterior_linpred 的平均值相同，这会产生transform = TRUE 时的概率。但是对于有限的抽奖次数，使用posterior_linpred 的噪音要小一些
您仍然会得到略有不同的数字。从概念上讲，您正在寻找这些边际效应来区分线性预测变量。