Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[paper] writing #332

Open
wants to merge 20 commits into
base: main
Choose a base branch
from
Open

[paper] writing #332

wants to merge 20 commits into from

Conversation

DominiqueMakowski
Copy link
Member

Let's get this bad boi out.

@mattansb @bwiernik what would you write about the two backends, emmeans and marginaleffects, their main differences etc.

Do they both rely on the delta method?

@DominiqueMakowski
Copy link
Member Author

@strengejacke I think I'll need your help for the marginal section, the different types of marginalization etc. it's quite blurry in my head to be honest ^^

@DominiqueMakowski DominiqueMakowski linked an issue Jan 13, 2025 that may be closed by this pull request
@DominiqueMakowski
Copy link
Member Author

@easystats/core-team remaining bits:

  • add details about types of marginalization
  • Drawbacks/benefits of both backends
  • Polishing

@mattansb
Copy link
Member

emmeans relies on the delta method for non-linear transformations when those are done as part of the estimation.

@DominiqueMakowski where do you want my input, and what should it include? I must admit I have 0 familiarity with modelbased - I use marginaleffects, emmeans and ggeffects directly for plotting and estimation, so I have no idea what would be relevant here.

@DominiqueMakowski
Copy link
Member Author

I must admit I have 0 familiarity

Your torment over you can make the switch now 😛

Mostly if you could take a look at the technical details section that talks about the backends

@strengejacke
Copy link
Member

strengejacke commented Jan 14, 2025

I must admit I have 0 familiarity with modelbased

modelbased currently covers classical predict() (non-focal held constant at representative values) or estimates marginal means (non-focal averaged/weighted), the latter using emmeans resp. marginaleffects via estimate_means(). estimate_means() should return identical results for both emmeans and marginaleffects backends, i.e. currently we support marginaleffects by mimicking the emmeans-approach of marginalzation.

@DominiqueMakowski DominiqueMakowski mentioned this pull request Jan 14, 2025

The modelbased package simplifies the extraction of these effects, providing a clear interface to understand how different predictors interact with outcomes in various scenarios.

[details about types of marginalization]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@strengejacke I think here we could mention our 2 main types of marginalization (essentially copypasta our nice docs on that)

Copy link
Member

@strengejacke strengejacke Jan 28, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can start working/helping on this paper by mid of February, or earlier if it's just minor stuff.

paper/paper.md Outdated
The term "Least-Squares Means" was somewhat misleading as it suggested a method specific to least-squares estimation, hence its renaming to `emmeans` in 2016 to clarify its scope for a wider range of models including generalized linear models, mixed models, and Bayesian models.
- `marginaleffects` (REF) was more recently introduced and also employs the delta method to approximate variance estimation. It is compatible with a wider range of models and allows for more flexibility in the specification of the marginal effects to be estimated.

[What's the difference / benefits / drawbacks of using one or ther other?]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mattansb any other interesting facts to mention here?

Copy link
Member

@strengejacke strengejacke Jan 28, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reference for marginaleffects would be https://www.jstatsoft.org/article/view/v111i09

Arel-Bundock, V., Greifer, N., & Heiss, A. (2024). How to Interpret Statistical Models Using marginaleffects for R and Python. Journal of Statistical Software, 111(9), 1–32. https://doi.org/10.18637/jss.v111.i09

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh there's a lot that I can add.

I think the biggest difference would be:

  • emmeans used a reference grid (that has by default a "complete" combination of all predictors) to generate expected values, by default conditioning on the mean of numeric values and marginalizing over categorical/binary variables, using a linear function of the model's coefficients (and v-cov matrix to get SEs) to give "predictions at the mean" (predictions for an average observation).
  • marginaleffects uses unit level predictions to generate two counterfactual values - the difference of which is then taken (with SEs computed using the delta method), and averaged across all units. By default, the original model frame is used.

Of course, emmeans can also the delta method and can build complex reference grids (that aren't actually "grid" like), and marginaleffects can also generate linear perditions at the mean.

Using the delta method is more computationally costly than using a linear combination (though marginaleffects is very efficient). Using linear combinations with orthogonal "grids" also often means that results from emmeans directly correspond to a models coefficients (which is a benefit for those who are used to looking at regressions tables to understand their models - this can be shown with an example).

...

Would you like me to add all of this in? Just some of this?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would you like me to add all of this in? Just some of this?

Anything you think is relevant, but I think it'll be good to be quite thorough and detailed here as the inner workings of marginal stuff are not often clearly explained so having details is good!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright - I'll give it some time tomorrow.

@mattansb
Copy link
Member

I've added the expanded "technical details" section.

Some general comments:

Statement of need

I think the "statement of need" needs some restructuring, because the actual need it addresses is that emmeans and marginaleffects are sometimes hard to use - "the modelbased package is designed to be user-friendly, with a focus on simplicity and flexibility". That is, it doesn't provide any functionality in and of itself, but it provides a (very!) convenient wrapper with good defaults.

Predictions

"For generalized linear models (GLMs), while the distinction between prediction and confidence intervals do not apply" - why not? I teach prediction intervals for non-gaussian models.
See also:
https://fromthebottomoftheheap.net/2017/05/01/glm-prediction-intervals-i/
https://fromthebottomoftheheap.net/2017/05/01/glm-prediction-intervals-ii/

Marginal effects

The different "marginals" need to be stated better here - the "marginal" in "marginal means" is not the same "marginal" that is in "marginal effects".
In fact, I'm not sure, from reading the paper, if estimate_contrasts() gives marginal differences or average differences? Etc..

@strengejacke
Copy link
Member

The different "marginals" need to be stated better here - the "marginal" in "marginal means" is not the same "marginal" that is in "marginal effects".
In fact, I'm not sure, from reading the paper, if estimate_contrasts() gives marginal differences or average differences? Etc..

Since I updated my mixed models slides for teaching after working on modelbased, maybe we can use some of the "definitions" I use on my slides (hopefully, these are accurate):

After fitting a model, it is useful generate model-based estimates. Usually, there are three quantities of interest:

  • marginal means: to better understand the relationship of predictors with the outcome. These are the average expected values, or adjusted predictions, of the response variable for different combinations of predictor values.
  • marginal effects: to evaluate the average strength of an effect (slope, average coefficient).
  • marginal contrasts: to look for (statistically significant) differences between groups, which helps us, for instance, analyzing “social inequalities”.

The key difference:

  • marginal means (or: adjusted predictions) return averages of the outcome variable, which allows you to say for instance “the average health score for a person at the age of sixty is 80 points”.
  • marginal effects return averages of coefficients, which allows you to say, for instance, “the average effect of age (or ageing) on the health score is an decrease 5 points”.
  • marginal contrasts return the average difference between two groups, typically indicated by different factor levels, which allows you to say “the average difference in health scores for lower and higher status groups at the same age is 15 points”.

Briefly:

  • marginal means show the association of the dependent variable and one or more independent variables;
  • marginal effects show you the average strength of a coefficient;
  • marginal contrasts show you average differences in the outcome between groups

@strengejacke
Copy link
Member

marginal differences or average differences?

What would be the difference between these differences?

@mattansb
Copy link
Member

mattansb commented Jan 29, 2025

library(marginaleffects)

mtcars$vs <- factor(mtcars$vs)
mtcars$gear <- factor(mtcars$gear)

mod <- glm(am ~ vs * hp * gear, 
           family = binomial(),
           data = mtcars)

# unit level differences
comparisons(mod, variables = "vs",
            type = "response")

# average differences (mean unit level differences)
avg_comparisons(mod, variables = "vs",
                type = "response")
#> 
#>  Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
#>    -0.166      0.195 -0.851    0.395 1.3 -0.547  0.216
#> 
#> Term: vs
#> Type:  response 
#> Comparison: 1 - 0
#> 




# marginal means over a regular grid
avg_predictions(mod, variables = "vs",
                # a regular grid
                newdata = datagrid(hp = mean, gear = levels),
                type = "response")

# marginal differences (differences between marginal means)
avg_predictions(mod, variables = "vs",
                # a regular grid
                newdata = datagrid(hp = mean, gear = levels),
                type = "response",
                hypothesis = ~ reference)
#> 
#>  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
#>       1 - 0   -0.302      0.239 -1.26    0.207 2.3 -0.772  0.167
#> 
#> Type:  response
#>  

@mattansb
Copy link
Member

I think we should avoid using the ambiguous term "marginal" as much as possible with regard to effects/slopes/differences (probably still okay with regard to "marginal means" since that can only be the "average" kind, not the "conditional" kind).

@strengejacke
Copy link
Member

(probably still okay with regard to "marginal means" since that can only be the "average" kind, not the "conditional" kind).

Unless you deal with mixed models, where you can have "marginal" and "conditional" effects/predictions/means (or whatever you'd like to call them).

@strengejacke
Copy link
Member

Ok, I see, I think your distinction refers to the marginalize argument. I personally think, and maybe we should state that clearly, that we use "marginal" in terms of "marginalizing out" / "integrating out" / "averaging over/across". Then, marginal means/effects/contrasts makes sense. Whether the means/effects/contrasts refer to a specific person, an average person or the "broader population" (counterfactuals) depends on how we "marginalizing out" / "integrating out" / "averaging over/across".

@strengejacke
Copy link
Member

Here are the modelbased equivalents:

library(modelbased)
library(marginaleffects)

mtcars$vs <- factor(mtcars$vs)
mtcars$gear <- factor(mtcars$gear)

mod <- glm(am ~ vs * hp * gear, family = binomial(), data = mtcars)

# average differences (mean unit level differences)
avg_comparisons(mod, variables = "vs")
estimate_contrasts(mod, "vs", marginalize = "population")

# marginal means over a regular grid
avg_predictions(mod, variables = "vs", newdata = datagrid(vs = levels, hp = mean, gear = levels))
estimate_means(mod, "vs", predict = "response")

# note the corrected prediction type
avg_predictions(mod, variables = "vs",
                newdata = datagrid(vs = levels, hp = mean, gear = levels),
                type = "invlink(link)")
estimate_means(mod, "vs")

# marginal differences (differences between marginal means)
avg_predictions(mod, variables = "vs",
                newdata = datagrid(hp = mean, gear = levels),
                type = "response",
                hypothesis = ~ reference)
estimate_contrasts(mod, "vs", comparison = ~ reference)

@mattansb
Copy link
Member

That's fine - but since there are these two meanings, I think it is best we stear clear of this term whenever possible. I would even try and find a better argument name for marginalize and explain in more detail what this argument does and why it's defaults are better / different than emmeans / marginaleffects'.

@DominiqueMakowski
Copy link
Member Author

a better argument name for marginalize

any suggestions?

@strengejacke
Copy link
Member

why it's defaults are better / different than emmeans / marginaleffects'.

I wouldn't say the defaults are better - since there are no "defaults" in marginaleffects, and the current default in modelbased is similar to emmeans results (despite using marginaleffects as backend by default).

find a better argument name for marginalize

How would you call the process of averaging non-focal terms?

Vincent proposes two definitions, "marginal" is either a minimal change of a derivative, or an average of some unit-level estimates (integral):
image

Since marginalize describes how to average over non-focal terms, I find the name quite ok, but I'm open for other suggestions.

@mattansb
Copy link
Member

since there are no "defaults" in marginaleffects

The default in marginaleffects is to use the model frame and produce unit level predictions/effects for each observation and potentially average across some/all groups of observations to get average / conditional effects.

It sounds like the default in modelbased is to make regular grids (like emmeans does by default). So that can also be a talking point in the paper.

Doesn't margenalize control the "grid type", similar to the string options of the newdata= argument in marginaleffects? Can we use that instead? grid_type= or grid=?

@strengejacke
Copy link
Member

Doesn't margenalize control the "grid type"

Not fully. That's why Dom and I think it's clearer to avoid a rather "technical" approach to what's going on, and instead a more "practical", or: "what's the question I'd like to answer?" approach.

marginalize has two options (producing different data grids) that either produce results like predict() (predictions for a specific subject, non-focals hold constant at mean and reference levels) or emmeans() (predictions for an average subject, non-focals hold constant at mean or averaged over factor levels), and one option that does not produce a grid, but returns "counterfactual" predictions (predictions for a broader / more representative population, extrapolation to a hypothetical population).

In all three cases you get an estimated average value for your response, however, for different "persons"/"subjects"/"groups" - because in all three cases, non-focal terms are treated differently (according to how they're "averaged" - or "marginalized").

That's the idea behind this argument name, and to avoid a "technical" thinking of marginal means, and instead looking at predictions in terms of "what is the sample/entity/person we're making inferences about"?

## Predictions

At a fundamental level, `modelbased` and similar packages leverage model *predictions*.
These predictions can be of different types, depending on the model and the question at hand.
For instance, for linear regressions, predictions can be associated with **confidence intervals** (`predict="expectation"`) or **prediction intervals** (`predict="prediction"`).
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not only linear regression? Remove that subordinary sentence?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes

It provides a unified interface to extract marginal means, marginal effects, and model predictions from a wide range of statistical models.
It is built on top of the `emmeans` (REF) and `marginaleffects` (REF) packages, which are the two most popular packages for extracting these quantities of interest from statistical models.
In line with the `easystats`' *raison d'être*, the `modelbased` package is designed to be user-friendly, with a focus on simplicity and flexibility.
The two probably most popular R packages for extracting these quantities of interest from statistical models are `emmeans` (REF) and `marginaleffects` (REF). These packages are enormously feature rich and (almost) cover all imaginable needs for post-hoc analysis of statistical models. However, these packages are not always easy to use, especially for users who are not familiar with the underlying statistical concepts. The `modelbased` package aims to unlock this currently underused potential by providing a unified interface to extract marginal means, marginal effects, contrasts, comparisons and model predictions from a wide range of statistical models. It is built on top of the two aforementioned `emmeans` and `marginaleffects` packages. In line with the `easystats`' *raison d'être*, the `modelbased` package is designed to be user-friendly, with a focus on simplicity and flexibility.


# Key concepts
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We might add a short sentence clarifying the conceptual differences between predictions and marginal means?


Mention types of predictions (expectations, predictions, response vs. the link scale)
## Marginal means and effects
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think somewhere we can introduce the concept of focal and non-focal terms. For marginalization, we can refer to the wording of "non-focal terms", which is probably better than "non-focal predictors" or "other predictors", because sometimes, we have interaction terms.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed

@mattansb
Copy link
Member

Hmmm to me the word "marginalize" is very technical, and marginalize = "average" is very confusing. First, both options are potentially computing an average, depending on the non-focal predictor. Second, both words (can) mean "to average".

Since the options are "average over regular grid" or "average over the model frame" and the word "marginalize" has two meanings, I would have:

  • average = "regular" / average = "grid" (what is currently marginalize = "average")
  • average = "population" (what is currently marginalize = "population")

(Or to be even more clear, I would have average_over= as the argument name.)

@strengejacke
Copy link
Member

Hmmm to me the word "marginalize" is very technical, and marginalize = "average" is very confusing. First, both options are potentially computing an average, depending on the non-focal predictor. Second, both words (can) mean "to average".

Since the options are "average over regular grid" or "average over the model frame" and the word "marginalize" has two meanings, I would have:

  • average = "regular" / average = "grid" (what is currently marginalize = "average")
  • average = "population" (what is currently marginalize = "population")

(Or to be even more clear, I would have average_over= as the argument name.)

Sounds good, we should think about this - @DominiqueMakowski WDYT?

@DominiqueMakowski
Copy link
Member Author

There's two decisions to be made I think:

  • Change name of argument: I'm somewhat positive about average_over even though it's a bit "out of scheme" and odd-looking. I prefer average_over much than average though. average_over / marginalize > average

  • Change "average" option: I don't like "grid" because it's not indicative of what it's used for, as Daniel mentioned, I like the idea of focusing on the use-case rather than the technical thing that happens. We toyed in the past with "specific", but maybe we could do like "average_individual" / "average_observation" even though that could be too specific and inaccurate?

@strengejacke
Copy link
Member

strengejacke commented Jan 30, 2025

I also prefer average_over. To clarify, what would be the new options / "translation" from current options? I think some of them would still fit, e.g.

marginalize = "average" -> average_over = "individual"
marginalize = "population" -> average_over = "population"

(and marginalize = "specific" -> average_over = "none")

@mattansb
Copy link
Member

mattansb commented Jan 30, 2025

I'm confused by marginalize = "average" having anything to do with an "individual" - what individual? Only when none of the non-focal predictors are factors will the give an estimate of a single theoretical individual - else you are averaging over a grid of marginal means.

Edit: Oh, you mean like "an average person"? (something that is weird to think about when averaging over a grid of marginal means)

how about:

marginalize = "average" -> average_over = "at mean"
marginalize = "population" -> average_over = "population"


What does marginalize = "specific" do? I haven't seen that documented.

@strengejacke
Copy link
Member

strengejacke commented Jan 30, 2025

What does marginalize = "specific" do? I haven't seen that documented.

That's like predict(newdata = <grid>). marginalize="average" is like emmeans default, and marginalize="population" is counterfactual (avg_predictions(variables=...)).

See also docs for ?estimate_means:

Screenshot_20250130-204101

@strengejacke
Copy link
Member

something that is weird to think about when averaging over a grid of marginal means

Why? How would you interpret the different types of averaging? I.e. to which subjects do your results refer to? (It's not helpful interpreting results in technical terms, I think, when a general audience should understand it)

@mattansb
Copy link
Member

Oh, now I get it - this isn't about averaging really, averaging is sometimes used, but it's about what units are the estimates generated (and then possibly averaged across).

So how about:

marginalize = "specific" -> estimate_over = "individual" (for a single unit)
marginalize = "average" -> estimate_over = "average" (for an average unit / average across marginal means)
marginalize = "population" -> estimate_over = "population" (for average across all units in a sample)

Or estimate_for= / over=

@strengejacke
Copy link
Member

strengejacke commented Jan 30, 2025

So how about:

Yes, that sounds appropriate, too! 🙂

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Contrasting slopes: differences between backends
3 participants