Contrasts on pairs of means and slopes

[samuelcharberet]

Hello,

Thank you very much for the package.

I currently work with mgcv::gam objects, which have formulae of the form (numeric ~ s(numeric, by = factor) + factor).

I would like to compare the slopes and mean for each factor in pairs.

Is this doable in marginaleffects ?

[vincentarelbundock]

Yes it almost certainly is.

Which tutorials have you read? What have you tried?

If you have not read any yet, I would start with these:

• https://marginaleffects.com/vignettes/predictions.html
• https://marginaleffects.com/vignettes/comparisons.html
• https://marginaleffects.com/vignettes/slopes.html

[samuelcharberet]

I have tried emmeans and emtrends, but I don't know how it can handle gams, still I did as follows to get statistical tests of comparison of mean and slopes between pairs of the levels of the factor

  emm <- emmeans::emmeans(gam, ~factor)
  emt <- emmeans::emtrends(gam, "factor", var = "continuous_predictor")

  con_emm = summary(emmeans::contrast(emm , method = "pairwise"))
  con_emt = summary(emmeans::contrast(emt , method = "pairwise"))

[vincentarelbundock]

OK, but this is the marginaleffects help forum, so not the right place to ask questions about emmeans.

I'm pretty sure marginaleffects can do what you need, but it would be nice if you could invest a bit of time to learn about package first.

The answer is probably already in one of the tutorials (we spent a lot of time writing). If it is not, reading them will empower you to ask better, more precise questions.

[samuelcharberet]

Thanks for your answer. Of course, I have already looked at the three vignettes you linked and listened to YouTube presentations of marginaleffects, but I did not find an answer to my question. I also described my problem more precisely in CrossValidated but received no answer.

I am not interested in doing predictions or compute slopes from my model. I would like to statistically compare the slopes and the means for every level of the factor included in the model. The emmean example serves as a explanation of the problem.

[vincentarelbundock]

You may be interested in this vignette, where I show several examples of emmeans vs. marginaleffects:

https://marginaleffects.com/vignettes/alternative_software.html

I think it may help to consider an actual minimal replicable example, instead of generic descriptions.

Here’s a GAM model, some emmeans code, and an exact replication with marginaleffects.

Is that what you were looking for?

library(mgcv)
library(emmeans)
library(marginaleffects)

dat <- transform(mtcars, am = factor(am))
mod <- gam(mpg ~ s(hp, by = am) + am, data = dat)

emm <- emmeans::emmeans(mod, ~am)
emt <- emmeans::emtrends(mod, "am", var = "hp")
summary(emmeans::contrast(emm, method = "pairwise"))
#  contrast  estimate  SE   df t.ratio p.value
#  am0 - am1    -3.19 1.5 26.9  -2.119  0.0434

summary(emmeans::contrast(emt, method = "pairwise"))
#  contrast  estimate     SE   df t.ratio p.value
#  am0 - am1   0.0126 0.0267 26.9   0.473  0.6397

predictions(mod, 
    newdata = datagrid(am = unique),
    hypothesis = "pairwise")
# 
#  Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
#      3.19        1.5 2.12   0.0341 4.9  0.24   6.13
# 
# Term: Row 1 - Row 2
# Type:  response 
# Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

slopes(mod, 
    variables = "hp",
    newdata = datagrid(am = unique),
    hypothesis = "pairwise")
# 
#  Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
#   -0.0127     0.0267 -0.477    0.634 0.7 -0.0651 0.0396
# 
# Term: Row 1 - Row 2
# Type:  response 
# Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high