Testing the difference between marginal effects calculated across factors

Question

I'm trying to test the difference between two marginal effects. I can get R to calculate the effects, but I can't find any resource explaining how to test their difference.

I've looked in the margins documentations and other marginal effects packages but have not been able to find something that tests the difference.

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

I want to produce a test that decides whether or not the marginal effects in each row of marg are significantly different; i.e., that the slopes in the marginal effects plots are different. This appears to be true because the confidence intervals do not overlap -- indicating that the effect of displacement is different for am=0 vs am=1.

We discuss in the comments below that we can test contrasts using emmeans, but that is a test of the average response across am=0 and 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 

Here the p-value is large indicating that average response when am=0 is not significantly different than when am=1.

Is it reasonable to do this (like testing the difference of two means)?

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

This p-value appears to agree with the analysis of non overlapping confidence intervals.

Answer 1

For the question of "Are the slopes statistically different, indicating that the effect of displacement is different for am=0 vs am=1?" question, you can get the p-value associated with the comparison directly from the summary of the lm() fit.

> 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

Notice that the p-value for the as.factor(am)1:disp term is 0.01044, which matches the output from pairs(emt) in Russ Lenth's answer.

(posting as an answer because insufficient reputation to post as a comment, yet)

Answer 2

If I understand your question, it can be answered using 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

Answer 3

Not sure, but probably you're looking at contrasts or pairwise comparisons of marginal effects? You can do this using the emmeans package:

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

Or simply use 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

^{Created on 2019-06-07 by the reprex package (v0.3.0)}