测试跨因素计算的边际效应之间的差异
Testing the difference between marginal effects calculated across factors
我正在尝试检验两个边际效应之间的差异。我可以让 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 是不同的。
我们在下面的评论中讨论了我们可以使用 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 值似乎与非重叠置信区间的分析一致。
不确定，但您可能正在查看边际效应的对比或成对比较？你可以使用 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-06-07 创建
如果我理解你的问题，可以使用 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
对于“斜率在统计上是否不同，表明 am=0 与 am=1 时位移的影响不同？”的问题。问题，您可以直接从 lm() 拟合的摘要中获得与比较相关的 p-value。
> 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-value 是 0.01044，这与 Russ Lenth 的回答中 pairs(emt) 的输出匹配。
（post 作为答案，因为声誉不足 post 作为评论，但）
我正在尝试检验两个边际效应之间的差异。我可以让 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 是不同的。
我们在下面的评论中讨论了我们可以使用 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 值似乎与非重叠置信区间的分析一致。
不确定，但您可能正在查看边际效应的对比或成对比较？你可以使用 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-06-07 创建如果我理解你的问题，可以使用 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
对于“斜率在统计上是否不同，表明 am=0 与 am=1 时位移的影响不同？”的问题。问题，您可以直接从 lm() 拟合的摘要中获得与比较相关的 p-value。
> 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-value 是 0.01044，这与 Russ Lenth 的回答中 pairs(emt) 的输出匹配。
（post 作为答案，因为声誉不足 post 作为评论，但）