Interoperability: {modelbased} and {marginaleffects}

[vincentarelbundock]

Hi all,

As most of you know, I recently released {marginaleffects}, a package with functionality that overlaps and (I think/hope) complements {modelbased} and {ggeffects}:

https://vincentarelbundock.github.io/marginaleffects/

Definitions

There seems to be a ton of variation in how vocabulary is used in this space, so here is how I define the 4 quantities that {marginaleffects} can compute:

• Marginal Effect (Vignette)
  • A partial derivative (slope) of the regression equation with respect to a continuous regressor of interest.
• Adjusted Prediction (Vignette)
  • The outcome predicted by a model for some combination of the regressors' values, such as their means or factor levels (a.k.a. "reference grid").
• Contrast (Vignette)
  • The difference between two adjusted predictions, calculated for meaningfully different regressor values (e.g., College graduates vs. Others).
• Marginal Mean (Vignette)
  • Adjusted predictions of a model, averaged across a "reference grid" of categorical predictors.

To produce all four of these quantities, the {marginaleffects} package requires a model object (mandatory) and a dataset (optional).

What does {modelbased} need?

I admit that I have very little experience with {modelbased}, so I wanted to ask you all: Is there anything in the list above that would be of interest here? Or do you already cover all this territory?

From a cursory review, my sense is that Adjusted Predictions are already very well-covered. There is probably nothing I can contribute here.

Marginal Means also seem well-covered. Are you doing the computation yourself, or wrapping emmeans::emmeans? Are you happy with this, or is there anything missing?

{modelbased} can also compute pairwise contrasts. Is that also by wrapping emmeans? In Issue https://github.com/easystats/modelbased/issues/150, some people expressed an interest in expanding the range of contrasts supported. This could be done by enriching the emmeans wrapper. Currently, {marginaleffects} only supports pairwise contrasts, but the infrastructure would make it relatively easy to allow users to specify specific contrasts. That said, I am unlikely to reach parity with emmeans in the medium run.

Finally, does {modelbased} support "marginal effects" as defined above (i.e., slopes / partial derivatives)? The website suggests that the answer is "yes", but I admit that I am a bit confused by the nomenclature, especially this section of your docs:

Marginal slopes are to numeric predictors what marginal means are to categorical predictors, in the sense that they can eventually be “averaged over” other predictors of the model. The key difference is that, while marginal means return averages of the outcome variable, which allows you to say for instance “the average reaction time in the C1 condition is 1366 ms”, marginal effects return averages of coefficients.

Your concept of "marginal effect" does not sound like the same as mine... In my definition, the partial derivative (slope) of the regression equation is the continuous analogue to a contrast for categorical predictors. Marginal effects measure the change inY that is associated with an infinitesimal change in a continuous X. A contrast measures the change in Y that is associated with a level change in a categorical X. But perhaps this is just nomenclature, and you already support this.

What does the future hold?

I opened this issue to know if there was anything I could bring to the table. My sense is that the slope computation is the most "mature" and useful aspect of my package, so I would guess that this is where most of the value lies. But if you feel that {modelbased} already covers the 4 quantities defined above adequately, then the answer is: "No, there is no added value." That's fine, and would be useful for me to know.

On the other hand, if there are some features that you would like to add but are having trouble with, I would be very curious to hear about it. Some of them might be easy to add.

[strengejacke]

Not sure, but I think that emtrends calculates marginal effects. Maybe you can add an example that compares marginaleffects to modelbased to see if results match? Here's a comparison:

# marginal effect

model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)

# own class/print
modelbased::estimate_slopes(model, trend = "Petal.Length", at = "Species")
#> Estimated Marginal Effects
#> 
#> Species    | Coefficient |   SE |        95% CI | t(144) |      p
#> -----------------------------------------------------------------
#> setosa     |        0.39 | 0.26 | [-0.13, 0.90] |   1.49 | 0.138 
#> versicolor |        0.37 | 0.10 | [ 0.18, 0.56] |   3.89 | < .001
#> virginica  |        0.23 | 0.08 | [ 0.07, 0.40] |   2.86 | 0.005 
#> Marginal effects estimated for Petal.Length

# emtrends, with modelbased API
modelbased::model_emtrends(model, trend = "Petal.Length", at = "Species")
#>  Species    Petal.Length.trend     SE  df lower.CL upper.CL
#>  setosa                  0.388 0.2602 144  -0.1264    0.902
#>  versicolor              0.374 0.0961 144   0.1843    0.564
#>  virginica               0.234 0.0819 144   0.0725    0.396
#> 
#> Confidence level used: 0.95

# emmeans::emtrends
emmeans::emtrends(model, c("Petal.Length", "Species"), "Petal.Length")
#>  Petal.Length Species    Petal.Length.trend     SE  df lower.CL upper.CL
#>          3.76 setosa                  0.388 0.2602 144  -0.1264    0.902
#>          3.76 versicolor              0.374 0.0961 144   0.1843    0.564
#>          3.76 virginica               0.234 0.0819 144   0.0725    0.396
#> 
#> Confidence level used: 0.95


# marginal mean

model <- lm(Sepal.Width ~ Species + Petal.Length, data = iris)

modelbased::estimate_means(model, at = "Petal.Length=c(1,4,6)")
#> Estimated Marginal Means
#> 
#> Petal.Length | Mean |   SE |       95% CI
#> -----------------------------------------
#> 1.00         | 2.23 | 0.17 | [1.90, 2.57]
#> 4.00         | 3.13 | 0.03 | [3.07, 3.19]
#> 6.00         | 3.73 | 0.14 | [3.45, 4.00]
#> 
#> Marginal means estimated at Petal.Length

modelbased::model_emmeans(model, at = "Petal.Length=c(1,4,6)")
#>  Petal.Length emmean     SE  df lower.CL upper.CL
#>             1   2.23 0.1688 146     1.90     2.57
#>             4   3.13 0.0296 146     3.07     3.19
#>             6   3.73 0.1380 146     3.45     4.00
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95

emmeans::emmeans(model, "Petal.Length", at = list(Petal.Length = c(1, 4, 6)))
#>  Petal.Length emmean     SE  df lower.CL upper.CL
#>             1   2.23 0.1688 146     1.90     2.57
#>             4   3.13 0.0296 146     3.07     3.19
#>             6   3.73 0.1380 146     3.45     4.00
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95



# adjusted predictions

model <- lm(Sepal.Width ~ Species + Petal.Length, data = iris)

vm <- modelbased::visualisation_matrix(iris, at = "Petal.Length=c(1,4,6)")
modelbased::estimate_expectation(model, vm)
#> Model-based Expectation
#> 
#> Petal.Length | Species | Predicted |   SE |       95% CI | Residuals
#> --------------------------------------------------------------------
#> 1.00         |  setosa |      3.29 | 0.05 | [3.19, 3.39] |      0.23
#> 4.00         |  setosa |      4.19 | 0.16 | [3.87, 4.50] |      1.13
#> 6.00         |  setosa |      4.78 | 0.28 | [4.23, 5.33] |      1.72
#> 
#> Variable predicted: Sepal.Width
#> Predictors modulated: Petal.Length=c(1,4,6)
#> Predictors controlled: Sepal.Length

predict(model, newdata = vm, interval = "confidence")
#>        fit      lwr      upr
#> 1 3.290180 3.186116 3.394244
#> 2 4.185113 3.869171 4.501055
#> 3 4.781735 4.232160 5.331310

^{Created on 2021-11-07 by the reprex package (v2.0.1)}

[vincentarelbundock]

Thanks!

The vignettes I linked to above compare adjusted predictions and marginal means explicitly to emmeans, so I know those can be replicated and are already covered by modelbased.

Focusing on marginal effects, it looks like your conjecture is correct: emtrends does indeed calculate marginal effects as I have defined them above. I can replicate your example like this:

library(marginaleffects)
library(modelbased)

model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)

mfx1 <- estimate_slopes(model,
                        trend = "Petal.Length",
                        at = "Species")

mfx2 <- marginaleffects(model,
                        variables = "Petal.Length",
                        newdata = typical(Species = unique(iris$Species)))

mfx1
# Estimated Marginal Effects
# 
# Species    | Coefficient |   SE |        95% CI | t(144) |      p
# -----------------------------------------------------------------
# setosa     |        0.39 | 0.26 | [-0.13, 0.90] |   1.49 | 0.138 
# versicolor |        0.37 | 0.10 | [ 0.18, 0.56] |   3.89 | < .001
# virginica  |        0.23 | 0.08 | [ 0.07, 0.40] |   2.86 | 0.005 
# Marginal effects estimated for Petal.Length

mfx2[, c("Species", "dydx", "std.error")]
#      Species      dydx  std.error
# 1     setosa 0.3878739 0.26016778
# 2 versicolor 0.3743068 0.09614966
# 3  virginica 0.2343482 0.08186667

But how would I go about estimating the marginal effect of wt on am when hp=100 and mpg=20, on the response or link scales?

mod <- glm(am ~ wt + hp + mpg, family = binomial, data = mtcars)

marginaleffects(mod, variables = "wt", newdata = typical(hp = 100, mpg = 20))
#   rowid     type term        dydx std.error      wt  hp mpg   predicted
# 1     1 response   wt -0.04098692 0.1340003 3.21725 100  20 0.005930419

marginaleffects(mod, variables = "wt", newdata = typical(hp = 100, mpg = 20), type = "link")
#   rowid type term      dydx std.error      wt  hp mpg predicted
# 1     1 link   wt -6.954924  3.352973 3.21725 100  20 -5.121712

And how do I get the “average marginal effects”, defined as:

Calculate the marginal effect of each variable for every combination of regressors actually observed in the dataset, and take the average of those marginal effects:

These "average marginal effects" are the default output in Stata and in the margins package. We get it like this:

marginaleffects(mod) |> summary()
# Average marginal effects 
#       type Term    Effect Std. Error z value   Pr(>|z|)     2.5 %  97.5 %
# 1 response   hp  0.003452   0.003342  1.0329    0.30165 -0.003098  0.0100
# 2 response  mpg  0.050576   0.065718  0.7696    0.44154 -0.078229  0.1794
# 3 response   wt -0.286129   0.046250 -6.1866 6.1489e-10 -0.376777 -0.1955
# 
# Model type:  glm 
# Prediction type:  response

[strengejacke]

marginal effects on the link-scale:

library(modelbased)
mod <- glm(am ~ wt + hp + mpg, family = binomial, data = mtcars)

estimate_slopes(mod,
                trend = "wt",
                at = list(hp = 100, mpg = 20))
#> Estimated Marginal Effects
#> 
#> hp     |   mpg | Coefficient |   SE |          95% CI |     z |  df |     p
#> ---------------------------------------------------------------------------
#> 100.00 | 20.00 |       -6.95 | 3.35 | [-13.53, -0.38] | -2.07 | Inf | 0.038
#> Marginal effects estimated for wt

@DominiqueMakowski @bwiernik - any ideas on average marginal effects? I think this is not yet implemented.

[vincentarelbundock]

That's super useful, thanks!

Maybe there's no value for you after all. Frankly, it is beginning to sound like I should have done more market research before jumping into this project ;)

The main things that might (or might not) still distinguish marginaleffects from emmeans are:

1. Average marginal effects
2. Model coverage potentially larger in marginaleffects

1 is trivial; 2 is probably not enough to justify much work on anyone's part.

[bwiernik]

Marginal effects are not yet implemented

[strengejacke]

Looks like marginal effects are indeed implemented (see https://github.com/easystats/modelbased/issues/158#issuecomment-962686661), at least on the link scale, but not average marginal effects, shown in the last example here: https://github.com/easystats/modelbased/issues/158#issuecomment-962676923

[DominiqueMakowski]

Actually modelbased computes something named "average marginal effects":

  model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)
  slopes <- modelbased::estimate_slopes(model, at = "Petal.Length", length = 50)
#> No numeric variable was specified for slope estimation. Selecting `trend = "Petal.Length"`.
  summary(slopes)
#> Average Marginal Effects
#> 
#> Start |  End | Petal.Length | Coefficient |   SE |         95% CI | t(142.33) |     p
#> -------------------------------------------------------------------------------------
#> 1.00  | 1.96 |         1.48 |        0.12 | 0.30 | [-0.47,  0.72] |      0.29 | 0.420
#> 2.08  | 3.05 |         2.57 |       -0.78 | 0.19 | [-1.15, -0.41] |     -4.25 | 0.001
#> 3.17  | 3.65 |         3.41 |       -0.10 | 0.26 | [-0.61,  0.42] |     -0.29 | 0.345
#> 3.77  | 4.25 |         4.01 |        0.54 | 0.20 | [ 0.15,  0.93] |      2.71 | 0.011
#> 4.37  | 6.90 |         5.64 |        0.07 | 0.23 | [-0.39,  0.53] |      0.58 | 0.430
#> Marginal effects estimated for Petal.Length
  plot(slopes)

^{Created on 2021-11-08 by the reprex package (v2.0.1)}

It's for derivatives for now, it computes the average effect over segments of significant vs. non-significant. In other words, it breaks the effect into parts and gives the average (linear) effect in that segment. It's very useful for reporting complex relationships such as the ones from GAMs

I feel like what marginaleffects means by "average marginal effect" can be obtained via emmeans, and hence via modelbased ? pinging emmeans expert @mattansb

[strengejacke]

The average marginal effect should just be one number, the "average overall one and only effect to rule them all", or something like this. Looking at splines, it seems obvious that an "average marginal effect" not always makes sense, but this seems to be true for some other things done in econometrics as well ;-)

[strengejacke]

library(marginaleffects)
model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)
marginaleffects(model) |> summary()
#> Average marginal effects 
#>       type         Term Effect Std. Error z value Pr(>|z|)   2.5 % 97.5 %
#> 1 response Petal.Length 0.2065    0.08426    2.45 0.014267 0.04133 0.3716
#> 
#> Model type:  gam 
#> Prediction type:  response

^{Created on 2021-11-08 by the reprex package (v2.0.1)}

[vincentarelbundock]

The average marginal effect should just be one number, the “average overall one and only effect to rule them all”, or something like this. Looking at splines, it seems obvious that an “average marginal effect” not always makes sense, but this seems to be true for some other things done in econometrics as well ;-)

I laughed. But, to be fair, I don’t think many economists (or political scientists) would argue that the average marginal effects makes sense in all context. It is mainly used in models (a) without interactions or splines, (b) where heterogeneous treatment effects are not the primary interest, but (c) functional form makes it so that the marginal effect depends on the values of regressors. In those cases, it doesn’t make sense to report the marginal effect for a fictional individual with characteristics set at the mean (the emmeans default), because that individual often does not (and could not) exist in real life. So we compute the marginal effects for every actually observed individuals, and then report the mean of that. Seems reasonable.

I would call the graph shown by @DominiqueMakowski a “Conditional Marginal Effects” plot. The y-axis is the partial derivative of Sepal.Width with respect to Petal.Length, and the x-axis is Petal.Length. That partial derivative of the GAM regression equation is different for different values of Petal.Length, so we get squiggles. The equivalent in marginaleffects would be:

library(marginaleffects)
model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris) 

plot_cme(model, effect = "Petal.Length", condition = "Petal.Length")

[vincentarelbundock]

The plot above appears to be mislabeled. It should be: "Marginal effect of Petal.Length on Sepal.Width"

[vincentarelbundock]

Don't think we need to get bogged down in a debate over Average Marginal Effects. This is really a trivial issue.

The only meaningful difference is that emmeans and marginaleffects cover different model types:

• emmeans models
• marginaleffects models

My infrastructures also makes it very easy to add support for new models, so the list will increase in the near future.

Unless range of supported models is a killer feature for you, I don't think it makes sense to move ahead with any interoperability or integration plans. I will thus close this issue. But feel free to add comments if you want; I'll keep my eye on it.

Thanks all for the helpful comments!

[strengejacke]

This is probably more a long-term goal, but given the possibly more light weight character of your package, I would say we should keep this issue open and investigate further about your suggested interoperability.

[vincentarelbundock]

This is probably more a long-term goal, but given the possibly more light weight character of your package, I would say we should keep this issue open and investigate further about your suggested interoperability.

Sounds good. I can also ping you if/when the promise of more supported models gets realized.

[DominiqueMakowski]

Many potential evolutions could be considered to avoid too many duplication and to pool maintaining efforts

• sharing some common core low-level set of functions (via insight?)
• modelbased swaping emmeans for marginaleffects for its workhorse
• marginaleffects putting some internals in modelbased and using it as its workhourse since it already relies on insight
• ...

we'll see how it goes :)

[vincentarelbundock]

Many potential evolutions could be considered to avoid too many duplication and to pool maintaining efforts

• sharing some common core low-level set of functions (via insight?)
• modelbased swaping emmeans for marginaleffects for its workhorse
• marginaleffects putting some internals in modelbased and using it as its workhourse since it already relies on insight
• ...

we'll see how it goes :)

Yes! All plausible options.

In the medium run, the vast majority of my efforts will be dedicated to developing model-specific methods and computation tools. More a substitute to emmeans than modelbased, so not much of a risk for duplication on that front.

But I agree that there are many possibilities, and that it will make sense to check back in in a few weeks/months to see how things develop.

[mattansb]

I have no idea what's going on so I'll just say: Yes, you can get true (average) marginal effects in emmeans with some trickery. I've discussed this with @vincentarelbundock on twitter (of all places!)

But I def like the idea of this collab here 🌻

[vincentarelbundock]

Sorry to revive a dead thread, but I wrote this above:

Sounds good. I can also ping you if/when the promise of more supported models gets realized.

I'm gearing up for a new release in the coming weeks and did a new accounting of supported models. Currently, marginaleffects can compute marginal effects (slopes) for 46 model types, including 22 model types that emmeans::emtrends does not currently support.

Caveats:

1. I might add two models in the emtrends column if/when I figure out how to deal with categorical outcome models.
2. marginaleffects is still not at feature parity in other dimensions, so I'm still not saying this is a slam dunk.

Here's a handy table which summarizes which models are supported by different packages, and which tables have been confirmed to yield numerically identical results under different packages:

• ✓: a green check means that the results of at least one model are equal to a reasonable tolerance.
• x: a red cross means that the results are not identical; extra caution is warranted.
• U: a grey U means that computing slopes for a model type is unsupported by alternative packages, but supported by marginaleffects.
• An empty cell means means that no comparison has been made yet.

[bwiernik]

This is a great (and impressive!) comparison! Can you remind me the input we need here?

[vincentarelbundock]

Can you remind me the input we need here?

Inputs:

• model object such as the output of MASS::polr()
• type string which is accept as an argument for the corresponding method predict.polr()
• newdata data frame which is the "grid" over which predictions are made or slopes calculated. This can be created easily with marginaleffects::datagrid or modelbased::visualization_matrix. If newdata is omitted, slopes are computed for every row of the original data (this allows users to easily compute "average marginal effects").
• Optional arguments to be pushed forward to predict(), such as re.form for specifying random effects in a lme4 or brms model.

This gives you 4 things under one unified interface:

1. marginaleffects: Slopes computed using numDeriv with standard errors obtained via the delta method. Similar to emmeans::emtrends.
2. Reverse pairwise contrasts for logical, factor, and character variables. Similar to emmeans::contrast.
3. predictions: Predicted outcomes for user-specified values of the regressors. This is essentially a wrapper around insight::get_predicted.
4. marginalmeans: Similar to emmeans::emmeans.

The next steps of development will probably involve:

1. Compatibility with insight::get_predicted's predict argument: expectation, link, prediction, etc.
2. More contrast types.
3. marginalmeans can only compute standard errors on the link scale.

[mattansb]

3. marginalmeans can only compute standard errors on the link scale.

Here is a simple formulation of the delta method (modified from msm::deltamethod()):

dmethod <- function (trans, .coef, .vcov) {
  
  .vcov <- as.matrix(.vcov)
  
  syms <- paste("x", seq_len(.coef), sep = "")
  g <- lapply(paste0(trans, "(", syms, ")"), reformulate)
  
  for (i in seq_len(.coef)) {
    assign(syms[i], .coef[i])
  }
  
  gdashmu <- t(sapply(g, function(form) {
    as.numeric(attr(eval(deriv(form, syms)), "gradient"))
  }))
  
  gdashmu %*% .vcov %*% t(gdashmu)
}



m <- glm(cyl ~ 1, mtcars, family = poisson())

sqrt(dmethod("exp", coef(m), vcov(m))) # SE on response scale
#>           [,1]
#> [1,] 0.4397264

emmeans::emmeans(m, ~ 1, trans = "response")
#>  1       rate   SE  df asymp.LCL asymp.UCL
#>  overall 6.19 0.44 Inf      5.33      7.05
#> 
#> Confidence level used: 0.95

^{Created on 2022-01-03 by the reprex package (v2.0.1)}

[vincentarelbundock]

Here is a simple formulation of the delta method (modified from msm::deltamethod()):

Very cool!