Расчет маржинальных эффектов в биномиальном логите с использованием rstanarm

Question

Я пытаюсь получить предельные эффекты, согласно этому посту: 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 не найден, что на самом деле означает, что мне тоже нужно пройти y в newdataчто не должно быть в соответствии с ?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

Answer 1

Для бинарной логит-модели предельный эффект непрерывной переменной является производной вероятности успеха по отношению к этой переменной, которая по правилу цепочки является логистической плотностью (оцениваемой при некоторых значениях предикторов, обычно наблюдаемых значений предикторы), умноженные на коэффициент рассматриваемой переменной. В вашем случае это было бы df <- as.data.frame(glm1) ME <- df$x2 * dlogis(posterior_linpred(glm1)) Поскольку это зависит от наблюдаемых значений предикторов, обычно усредняют по данным AME <- rowMeans(ME) В случае двоичного предиктора вы можете просто вычесть вероятность успеха, когда x1 = 0 от вероятности успеха, когда x1 = 1 с помощью 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)