计算边际效应：为什么ggeffect和ggemmeans给出不同的答案？

Question

计算边际效应：为什么ggeffect和ggemmeans给出不同的答案？

5

例子

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包在计算过程中的内部差异导致的？此外，有没有人知道如何从头开始计算类似示例模型的边际效应的资源？

- Constantin

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- dipetkov · Accepted Answer

你拟合了一个复杂模型：具有随机效应的零膨胀负二项式模型。

观察结果与模型规范无关。通过拟合一个更简单的模型来展示这一点：仅具有固定效应的泊松模型。

library("glmmTMB")
library("ggeffects")

m <- glmmTMB(
  count ~ spp + mined,
  family = poisson,
  data = Salamanders
)

ggemmeans(m, terms = "spp")
#> # Predicted counts of count
#> 
#> spp   | Predicted |       95% CI
#> --------------------------------
#> GP    |      0.73 | [0.59, 0.89]
#> PR    |      0.18 | [0.12, 0.27]
#> DM    |      0.91 | [0.76, 1.10]
#> EC-A  |      0.34 | [0.25, 0.45]
#> EC-L  |      1.35 | [1.15, 1.59]
#> DES-L |      1.43 | [1.22, 1.68]
#> DF    |      0.79 | [0.64, 0.96]

ggeffect(m, terms = "spp")
#> # Predicted counts of count
#> 
#> spp   | Predicted |       95% CI
#> --------------------------------
#> GP    |      0.76 | [0.62, 0.93]
#> PR    |      0.19 | [0.13, 0.28]
#> DM    |      0.96 | [0.79, 1.15]
#> EC-A  |      0.35 | [0.26, 0.47]
#> EC-L  |      1.41 | [1.20, 1.66]
#> DES-L |      1.50 | [1.28, 1.75]
#> DF    |      0.82 | [0.67, 1.00]

文档说明内部，ggemmeans() 调用 emmeans::emmeans()，而 ggeffect() 则调用 effects::Effect()。两者都计算边际效应，但它们在如何对 mined 进行边缘化（即平均）以获得 spp 的效应时做了不同的（默认）选择。 mined 是一个具有两个级别（“是”和“否”）的分类变量。关键是这两个级别不平衡：有比“是”更多的“否”。

xtabs(~ mined + spp, data = Salamanders)
#>      spp
#> mined GP PR DM EC-A EC-L DES-L DF
#>   yes 44 44 44   44   44    44 44
#>   no  48 48 48   48   48    48 48

直觉上，这意味着在“挖掘”之后加权平均值[考虑 (44×yes+48×no)/92]与在“挖掘”之后简单平均值[考虑(yes+no)/2]不同。

让我们通过指定如何边缘化出“挖掘”的方式来检查这种直觉，当我们直接调用emmeans::emmeans()。

# mean (default)
emmeans::emmeans(m, "spp", type = "response", weights = "equal")
#>  spp    rate     SE  df lower.CL upper.CL
#>  GP    0.726 0.0767 636    0.590    0.893
#>  PR    0.181 0.0358 636    0.123    0.267
#>  DM    0.914 0.0879 636    0.757    1.104
#>  EC-A  0.336 0.0497 636    0.251    0.449
#>  EC-L  1.351 0.1120 636    1.148    1.590
#>  DES-L 1.432 0.1163 636    1.221    1.679
#>  DF    0.786 0.0804 636    0.643    0.961
#> 
#> Results are averaged over the levels of: mined 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the log scale

# weighted mean
emmeans::emmeans(m, "spp", type = "response", weights = "proportional")
#>  spp    rate     SE  df lower.CL upper.CL
#>  GP    0.759 0.0794 636    0.618    0.932
#>  PR    0.190 0.0373 636    0.129    0.279
#>  DM    0.955 0.0909 636    0.793    1.152
#>  EC-A  0.351 0.0517 636    0.263    0.469
#>  EC-L  1.412 0.1153 636    1.203    1.658
#>  DES-L 1.496 0.1196 636    1.279    1.751
#>  DF    0.822 0.0832 636    0.674    1.003
#> 
#> Results are averaged over the levels of: mined 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the log scale

第二个选项返回使用ggeffects::ggeffect计算的边际效应。

更新

@Daniel 指出，ggeffects接受weights参数，并将其传递给emmeans。这样，您可以继续使用ggeffects，同时仍然控制如何对预测进行平均以计算边际效应。

请使用以下内容尝试：

ggemmeans(m, terms="spp", weights = "proportional")
ggemmeans(m, terms="spp", weights = "equal")