2

我正在等待bife（二进制固定效果）包得到支持，ggeffects但与此同时，我想知道如何复制 ggeffects 的功能，并且（非常第二）想了解我是否应该使用平均预测效应 (APE) 或预测概率。

在正常情况下，我会运行 aglm然后使用ggeffects.

testglm <- iris %>%
  mutate(LengthDummy = if_else(Sepal.Length > 5, "Long", "Short") %>% as_factor(.)) %>%
  glm(LengthDummy ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = .,  family = binomial)
ggeffects::ggeffect(testglm, terms = "Petal.Length")

你应该得到这张表：

# Predicted probabilities of LengthDummy
Petal.Length | Predicted |       95% CI
---------------------------------------
        1.00 |      1.00 | [0.41, 1.00]
        1.50 |      1.00 | [0.24, 1.00]
        2.50 |      0.58 | [0.04, 0.98]
        3.50 |      0.01 | [0.00, 0.09]
        4.00 |      0.00 | [0.00, 0.03]
        4.50 |      0.00 | [0.00, 0.01]
        5.50 |      0.00 | [0.00, 0.00]
        7.00 |      0.00 | [0.00, 0.00]

这个输出，其中 Petal.Length 是数据集中 (8) 个值的范围，并且所有其他模型预测变量都保持在它们的平均值，生成并且可以通过管道传输到plot()( ggeffects::ggeffect(testglm, terms = "Petal.Length")%>% plot()) 或到ggplot2()。

ggeffects::ggeffect(testglm, terms = "Petal.Length") %>% 
ggplot(., aes(x, predicted))+ 
  geom_line()

然而，这一次，我的模型（显然不是 iris 项目，这是一个示例数据集）需要 FE。所以，我的问题是，我如何像glmggeffects一样“从 logit 模型的预测变量中计算所有可能级别和值的预测值”？bife

library(bife)
fetest <- iris%>% mutate(LengthDummy= if_else(Sepal.Length>5, "Long", "Short") %>% as_factor(.)) %>%  
                 bife(LengthDummy~ Sepal.Width+ Petal.Length+Petal.Width | Species, data = ., "logit")

我可以获得平均预测效果，它可以使用APEsIris <- summary(get_APEs(fetest)). 这当然会产生 Petal.Length 的单一平均值，而不是 8。我可以用 se 和 CI 计算数据帧，但它不是等价的（单一估计）


    object <- APEsIris
    estimate <- object[["delta"]]
    se <- sqrt(diag(object[["vcov"]]))
    z <- estimate/se
    p <- 2 * pnorm(-abs(z))
    conf.low <- estimate-1.96*se
    conf.high <- estimate+1.96*se
    term <- names(estimate)
    APEsdf <- data.frame(term, estimate, se, conf.low, conf.high, z, p)
    colnames(cm) <- c("term", "estimate", 
                      "std.error", "conf.low", "conf.high", "statistic", "p.value")
    #colnames(cm) <- c("term", "estimate", "Std. error", "Lower bound", "Upper bound", "statistic", "Pr(> |z|)")
    rownames(cm) <- NULL
    APEsdf <- as.data.frame(APEsdf)

那么预测呢？

predictedbifes<-predict.bife(fetest)在数据集中产生 150 个预测，即观察数。即使我将它们绘制成图表，我也不明白我将如何操纵预测以随着 Petal.Length 的观察值而变化，并将所有其他预测变量保持在平均值。

简而言之，我怎样才能重现ggeffects计算以产生一个边际效应表，我可以用它来绘制这个固定效应 logit 模型中的不同变量？

4

0 回答 0