测试跨因子计算的边际效应之间的差异答案

【问题标题】：Testing the difference between marginal effects calculated across factors测试跨因子计算的边际效应之间的差异
【发布时间】：2019-06-07 16:09:32
【问题描述】：

我正在尝试测试两个边际效应之间的差异。我可以让 R 来计算效果，但我找不到任何资源来解释如何测试它们的差异。

我查看了边距文档和其他边际效应包，但找不到可以测试差异的东西。

data("mtcars")

mod<-lm(mpg~as.factor(am)*disp,data=mtcars)

(marg<-margins(model = mod,at = list(am = c("0","1"))))
 at(am)     disp    am1
      0 -0.02758 0.4518
      1 -0.05904 0.4518
summary(marg)
 factor     am     AME     SE       z      p   lower   upper
    am1 1.0000  0.4518 1.3915  0.3247 0.7454 -2.2755  3.1791
    am1 2.0000  0.4518 1.3915  0.3247 0.7454 -2.2755  3.1791
   disp 1.0000 -0.0276 0.0062 -4.4354 0.0000 -0.0398 -0.0154
   disp 2.0000 -0.0590 0.0096 -6.1353 0.0000 -0.0779 -0.0402

我想生成一个测试来确定每行 marg 中的边际效应是否显着不同；即，边际效应图中的斜率不同。这似乎是正确的，因为置信区间不重叠 - 表明位移的影响对于 am=0 和 am=1 是不同的。

我们在下面的 cmets 中讨论了我们可以使用 emmeans 测试对比，但这是对 am=0 和 am=1 的平均响应的测试。

emm<-emmeans(mod,~ as.factor(am)*disp)
emm
 am disp emmean    SE df lower.CL upper.CL
  0  231   18.8 0.763 28     17.2     20.4
  1  231   19.2 1.164 28     16.9     21.6
cont<-contrast(emm,list(`(0-1)`=c(1,-1)))
cont
 contrast estimate   SE df t.ratio p.value
 (0-1)      -0.452 1.39 28 -0.325  0.7479

这里的 p 值很大，表明 am=0 时的平均响应与 am=1 时的平均响应没有显着差异。

这样做是否合理（例如测试两种方法的差异）？

smarg<-summary(marg)
(z=as.numeric((smarg$AME[3]-smarg$AME[4])/sqrt(smarg$SE[3]^2+smarg$SE[4]^2)))
[1] 2.745
2*pnorm(-abs(z))
[1] 0.006044

这个 p 值似乎与非重叠置信区间的分析一致。

【问题讨论】：

请不要使用重叠的置信区间来比较均值。这是不合适的。均值的标准误与均值的差异不同，而且它们不匹配！

标签： r marginal-effects

【解决方案1】：

如果我理解你的问题，可以使用emtrends回答：

library(emmeans)

emt = emtrends(mod, "am", var = "disp")

emt  # display the estimated slopes

## am disp.trend      SE df lower.CL upper.CL
##  0    -0.0276 0.00622 28  -0.0403  -0.0148
##  1    -0.0590 0.00962 28  -0.0787  -0.0393
##
## Confidence level used: 0.95 

pairs(emt) # test the difference of slopes

## contrast estimate     SE df t.ratio p.value
## 0 - 1      0.0315 0.0115 28 2.745   0.0104

【讨论】：

【解决方案2】：

对于“斜率在统计上是否不同，表明位移对 am=0 和 am=1 的影响不同”的问题？问题，您可以直接从 lm() 拟合的摘要中获取与比较相关的 p 值。


> summary(mod)

Call:
lm(formula = mpg ~ as.factor(am) * disp, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.6056 -2.1022 -0.8681  2.2894  5.2315 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         25.157064   1.925053  13.068 1.94e-13 ***
as.factor(am)1       7.709073   2.502677   3.080  0.00460 ** 
disp                -0.027584   0.006219  -4.435  0.00013 ***
as.factor(am)1:disp -0.031455   0.011457  -2.745  0.01044 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.907 on 28 degrees of freedom
Multiple R-squared:  0.7899,    Adjusted R-squared:  0.7674 
F-statistic: 35.09 on 3 and 28 DF,  p-value: 1.27e-09

请注意，as.factor(am)1:disp 术语的 p 值为 0.01044，这与 Russ Lenth 的答案中 pairs(emt) 的输出相匹配。

（作为答案发布，因为声誉不足以作为评论发布）

【讨论】：

【解决方案3】：

不确定，但您可能正在研究边际效应的对比或成对比较？您可以使用 emmeans 包来做到这一点：

library(margins)
library(emmeans)
library(magrittr)

data("mtcars")
mod <- lm(mpg ~ as.factor(am) * disp, data = mtcars)

marg <- margins(model = mod, at = list(am = c("0", "1")))
marg
#> Average marginal effects at specified values
#> lm(formula = mpg ~ as.factor(am) * disp, data = mtcars)
#>  at(am)     disp    am1
#>       0 -0.02758 0.4518
#>       1 -0.05904 0.4518

emmeans(mod, c("am", "disp")) %>% 
  contrast(method = "pairwise")
#>  contrast                    estimate   SE df t.ratio p.value
#>  0,230.721875 - 1,230.721875   -0.452 1.39 28 -0.325  0.7479

emmeans(mod, c("am", "disp")) %>% 
  contrast()
#>  contrast            estimate    SE df t.ratio p.value
#>  0,230.721875 effect   -0.226 0.696 28 -0.325  0.7479 
#>  1,230.721875 effect    0.226 0.696 28  0.325  0.7479 
#> 
#> P value adjustment: fdr method for 2 tests

或者干脆使用summary():

library(margins)

data("mtcars")
mod <- lm(mpg ~ as.factor(am) * disp, data = mtcars)

marg <- margins(model = mod, at = list(am = c("0", "1")))
marg
#> Average marginal effects at specified values
#> lm(formula = mpg ~ as.factor(am) * disp, data = mtcars)
#>  at(am)     disp    am1
#>       0 -0.02758 0.4518
#>       1 -0.05904 0.4518

summary(marg)
#>  factor     am     AME     SE       z      p   lower   upper
#>     am1 1.0000  0.4518 1.3915  0.3247 0.7454 -2.2755  3.1791
#>     am1 2.0000  0.4518 1.3915  0.3247 0.7454 -2.2755  3.1791
#>    disp 1.0000 -0.0276 0.0062 -4.4354 0.0000 -0.0398 -0.0154
#>    disp 2.0000 -0.0590 0.0096 -6.1353 0.0000 -0.0779 -0.0402

^{由reprex package (v0.3.0) 于 2019 年 6 月 7 日创建}

【讨论】：

我到了 summary(marg) 部分，但我想表明 disp(at am=0) - disp (at am=1) 与 0 显着不同。在这里，查看置信区间，这似乎是真的，但我不知道如何要求 R 进行测试。
是的，这就是 emmeans 所做的（参见我上面的示例）。另请参阅this vignette。
这不是回答不同的问题吗？例如，emmeans 不是告诉我们，对于 am=1（相对于 am=0），平均响应要大 0.32 吗？这并没有告诉我们边际效应图的斜率不同，还是这样？
我在上面（在原始帖子中）添加了更多代码和输出以阐明我的意思。