Search code examples
rmarginal-effects

Testing the difference between marginal effects calculated across factors


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.


Solution

  • 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