Open vincentarelbundock opened 3 years ago
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)
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
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.
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:
marginaleffects
1 is trivial; 2 is probably not enough to justify much work on anyone's part.
Marginal effects are not yet implemented
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
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
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 ;-)
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)
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")
The plot above appears to be mislabeled. It should be: "Marginal effect of Petal.Length on Sepal.Width"
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:
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!
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.
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.
Many potential evolutions could be considered to avoid too many duplication and to pool maintaining efforts
we'll see how it goes :)
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.
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 🌻
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:
emtrends
column if/when I figure out how to deal with categorical outcome models.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:
marginaleffects
.This is a great (and impressive!) comparison! Can you remind me the input we need here?
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").predict()
, such as re.form
for specifying random effects in a lme4
or brms
model.This gives you 4 things under one unified interface:
marginaleffects
: Slopes computed using numDeriv
with standard errors obtained via the delta method. Similar to emmeans::emtrends
.emmeans::contrast
.predictions
: Predicted outcomes for user-specified values of the regressors. This is essentially a wrapper around insight::get_predicted
.marginalmeans
: Similar to emmeans::emmeans
.The next steps of development will probably involve:
insight::get_predicted
's predict
argument: expectation
, link
, prediction
, etc.marginalmeans
can only compute standard errors on the link scale.
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)
Here is a simple formulation of the delta method (modified from
msm::deltamethod()
):
Very cool!
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:
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 theemmeans
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 withemmeans
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:
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.