Con*_*tin 5 r mixed-models marginal-effects glmmtmb
library(glmmTMB)
library(ggeffects)
## Zero-inflated negative binomial model
(m <- glmmTMB(count ~ spp + mined + (1|site),
ziformula=~spp + mined,
family=nbinom2,
data=Salamanders,
na.action = "na.fail"))
summary(m)
ggemmeans(m, terms="spp")
spp | Predicted | 95% CI
--------------------------------
GP | 1.11 | [0.66, 1.86]
PR | 0.42 | [0.11, 1.59]
DM | 1.32 | [0.81, 2.13]
EC-A | 0.75 | [0.37, 1.53]
EC-L | 1.81 | [1.09, 3.00]
DES-L | 2.00 | [1.25, 3.21]
DF | 0.99 | [0.61, 1.62]
ggeffects::ggeffect(m, terms="spp")
spp | Predicted | 95% CI
--------------------------------
GP | 1.14 | [0.69, 1.90]
PR | 0.44 | [0.12, 1.63]
DM | 1.36 | [0.85, 2.18]
EC-A | 0.78 | [0.39, 1.57]
EC-L | 1.87 | [1.13, 3.07]
DES-L | 2.06 | [1.30, 3.28]
DF | 1.02 | [0.63, 1.65]
为什么 ggeffect 和 ggemmeans 给出的边际效应结果不同？这是否只是包 emmeans 和effects 计算它们的内部方式？另外，有谁知道一些关于如何从头开始计算示例中模型的边际效应的资源吗？
您拟合了一个复杂模型：具有随机效应的零膨胀负二项式模型。
\n您观察到的情况与型号规格关系不大。让我们通过拟合一个更简单的模型来展示这一点：仅具有固定效应的泊松模型。
\nlibrary("glmmTMB")\nlibrary("ggeffects")\n\nm <- glmmTMB(\n count ~ spp + mined,\n family = poisson,\n data = Salamanders\n)\n\nggemmeans(m, terms = "spp")\n#> # Predicted counts of count\n#> \n#> spp | Predicted | 95% CI\n#> --------------------------------\n#> GP | 0.73 | [0.59, 0.89]\n#> PR | 0.18 | [0.12, 0.27]\n#> DM | 0.91 | [0.76, 1.10]\n#> EC-A | 0.34 | [0.25, 0.45]\n#> EC-L | 1.35 | [1.15, 1.59]\n#> DES-L | 1.43 | [1.22, 1.68]\n#> DF | 0.79 | [0.64, 0.96]\n\nggeffect(m, terms = "spp")\n#> # Predicted counts of count\n#> \n#> spp | Predicted | 95% CI\n#> --------------------------------\n#> GP | 0.76 | [0.62, 0.93]\n#> PR | 0.19 | [0.13, 0.28]\n#> DM | 0.96 | [0.79, 1.15]\n#> EC-A | 0.35 | [0.26, 0.47]\n#> EC-L | 1.41 | [1.20, 1.66]\n#> DES-L | 1.50 | [1.28, 1.75]\n#> DF | 0.82 | [0.67, 1.00]\n文档解释了内部ggemmeans()调用emmeans::emmeans()whileggeffect()调用effects::Effect().
emmeans和effects都计算边际效应，但它们做出了不同的（默认）选择，如何边缘化（即平均）挖掘以获得spp的效果。
\n挖掘是一个具有两个级别的分类变量：“是”和“否”。关键是这两个层面不平衡：“不”的数量略多于“是”的数量。
\nxtabs(~ mined + spp, data = Salamanders)\n#> spp\n#> mined GP PR DM EC-A EC-L DES-L DF\n#> yes 44 44 44 44 44 44 44\n#> no 48 48 48 48 48 48 48\n直观上，这意味着过度挖掘的加权平均值[想想(44 \xc3\x97 yes + 48 \xc3\x97 no) / 92 ] 与过度挖掘的简单平均值不同[想想(yes + no) / 2 ]。
\n让我们通过指定直接调用时如何边缘化挖掘emmeans::emmeans()来检查直觉。
# mean (default)\nemmeans::emmeans(m, "spp", type = "response", weights = "equal")\n#> spp rate SE df lower.CL upper.CL\n#> GP 0.726 0.0767 636 0.590 0.893\n#> PR 0.181 0.0358 636 0.123 0.267\n#> DM 0.914 0.0879 636 0.757 1.104\n#> EC-A 0.336 0.0497 636 0.251 0.449\n#> EC-L 1.351 0.1120 636 1.148 1.590\n#> DES-L 1.432 0.1163 636 1.221 1.679\n#> DF 0.786 0.0804 636 0.643 0.961\n#> \n#> Results are averaged over the levels of: mined \n#> Confidence level used: 0.95 \n#> Intervals are back-transformed from the log scale\n\n# weighted mean\nemmeans::emmeans(m, "spp", type = "response", weights = "proportional")\n#> spp rate SE df lower.CL upper.CL\n#> GP 0.759 0.0794 636 0.618 0.932\n#> PR 0.190 0.0373 636 0.129 0.279\n#> DM 0.955 0.0909 636 0.793 1.152\n#> EC-A 0.351 0.0517 636 0.263 0.469\n#> EC-L 1.412 0.1153 636 1.203 1.658\n#> DES-L 1.496 0.1196 636 1.279 1.751\n#> DF 0.822 0.0832 636 0.674 1.003\n#> \n#> Results are averaged over the levels of: mined \n#> Confidence level used: 0.95 \n#> Intervals are back-transformed from the log scale\n第二个选项返回用 计算的边际效应ggeffects::ggeffect。
更新
\n@Daniel 指出ggeffects接受weights参数并将其传递给emmeans. 通过这种方式，您可以继续使用ggeffects并仍然控制如何平均预测来计算边际效应。
自己尝试一下：
\nggemmeans(m, terms="spp", weights = "proportional")\nggemmeans(m, terms="spp", weights = "equal")\n