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

Question

我试图根据这个 post 获得边际效应：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)])因为根据上面的post（据我了解），为了计算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

Answer 1

对于二元对数模型，连续变量的边际效应是该变量成功概率的导数，根据链式法则，它是逻辑密度（在预测变量的某些值下评估，通常是预测变量的观察值）乘以相关变量的系数。在你的情况下，那将是 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)

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

Calculating marginal effects in binomial logit using rstanarm

r

stan

rstan

rstanarm