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.

Sign up for GitHub

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

Jump to bottom

estimate_means(backend = "marginaleffects") not calculating marginal means #273

Open
strengejacke opened this issue Dec 13, 2024 · 12 comments · May be fixed by #275
Open

estimate_means(backend = "marginaleffects") not calculating marginal means #273

strengejacke opened this issue Dec 13, 2024 · 12 comments · May be fixed by #275
Assignees
@strengejacke

Comments

@strengejacke
Copy link
Member

See:

data(Salamanders, package = "glmmTMB")
m <- glm(count ~ mined + spp, family = poisson(), data = Salamanders)
modelbased::estimate_means(m, "mined")
#> Estimated Marginal Means
#> 
#> mined | rate |   SE |       95% CI
#> ----------------------------------
#> yes   | 0.24 | 0.03 | [0.20, 0.30]
#> no    | 1.86 | 0.08 | [1.70, 2.03]
#> 
#> Marginal means estimated at mined
modelbased::estimate_means(m, "mined", backend = "marginaleffects")
#> Estimated Marginal Means
#> 
#> mined | spp | Mean |   SE |       95% CI
#> ----------------------------------------
#> yes   |  GP | 0.26 | 0.04 | [0.19, 0.33]
#> no    |  GP | 2.01 | 0.19 | [1.63, 2.39]
#> 
#> Marginal means estimated at mined
modelbased::estimate_expectation(m, insight::get_datagrid(m, "mined"))
#> Model-based Expectation
#> 
#> mined | Predicted |   SE |       95% CI
#> ---------------------------------------
#> yes   |      0.26 | 0.04 | [0.20, 0.34]
#> no    |      2.01 | 0.19 | [1.66, 2.43]
#> 
#> Variable predicted: count
#> Predictors modulated: mined
#> Predictors controlled: spp (GP)

Created on 2024-12-13 with reprex v2.1.1

I think we need avg_predictions() here?
@DominiqueMakowski
Copy link
Member

modelbased/R/get_marginalmeans.R

Lines 19 to 24 in 46fbc0d

if (marginal) {
means <- marginaleffects::predictions(model,
newdata = insight::get_data(model),
by = at_specs$varname,
conf_level = ci
)

We don't use avg_predictions at all it seems, should we change that line to use it? I still don't have a good grasp on the API on marginaleffects so let's do what you think is correct to try to be consistent between emmeans, marginaleffects, ggeffects and modelbased

@strengejacke
Copy link
Member Author

I suggest we try to mimic emmeans. Not sure about mixed models, I think it works be cool to marginalize over random effects here, which is slightly different from emmeans behavior, I think

@strengejacke
Copy link
Member Author

I'm also still not 100% sure what marginaleffects does in which situation, but I suggest following:

marginaleffects::avg_predictions(m, by = "focalterm", type = "response")
data(Salamanders, package = "glmmTMB")

# marginal means - similar to emmeans
m <- glm(count ~ mined + spp, family = poisson(), data = Salamanders)
marginaleffects::avg_predictions(m, by = "mined", type = "invlink(link)")
#> 
#>  mined Estimate Pr(>|z|)     S 2.5 % 97.5 %
#>    yes    0.243   <0.001 127.9 0.196   0.30
#>    no     1.859   <0.001 141.7 1.702   2.03
#> 
#> Type:  invlink(link) 
#> Columns: mined, estimate, p.value, s.value, conf.low, conf.high
modelbased::estimate_means(m, "mined")
#> Estimated Marginal Means
#> 
#> mined | rate |   SE |       95% CI
#> ----------------------------------
#> yes   | 0.24 | 0.03 | [0.20, 0.30]
#> no    | 1.86 | 0.08 | [1.70, 2.03]
#> 
#> Marginal means estimated at mined

# marginal means - averaged over random effects (not similar to emmeans)
m <- glmmTMB::glmmTMB(count ~ mined + spp + (1 | site), family = poisson(), data = Salamanders)
marginaleffects::avg_predictions(m, by = "mined", type = "response")
#> 
#>  mined Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#>    yes    0.292     0.0651 4.48   <0.001 17.1 0.164  0.419
#>    no     2.264     0.3862 5.86   <0.001 27.7 1.507  3.021
#> 
#> Type:  response 
#> Columns: mined, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
modelbased::estimate_means(m, "mined")
#> Estimated Marginal Means
#> 
#> mined | rate |   SE |       95% CI
#> ----------------------------------
#> yes   | 0.18 | 0.04 | [0.12, 0.28]
#> no    | 1.75 | 0.30 | [1.25, 2.46]
#> 
#> Marginal means estimated at mined

# counterfactual - simular to simulate response
ggeffects::predict_response(m, "mined", type = "simulate")
#> # Predicted counts of count
#> 
#> mined | Predicted |     95% CI
#> ------------------------------
#> yes   |      0.26 | 0.00, 1.69
#> no    |      2.54 | 0.00, 8.11
marginaleffects::avg_predictions(m, variables = "mined", type = "response")
#> 
#>  mined Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#>    yes    0.262     0.0585 4.48   <0.001 17.1 0.148  0.377
#>    no     2.524     0.4306 5.86   <0.001 27.7 1.680  3.368
#> 
#> Type:  response 
#> Columns: mined, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

# for linear mixed models, emmeans and marginaleffects same
data(sleepstudy, package = "lme4")
m <- lme4::lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
marginaleffects::avg_predictions(m, by = "Days", type = "response")
#> 
#>  Days Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
#>     0      251       6.82 36.8   <0.001 984.4   238    265
#>     1      262       6.79 38.6   <0.001   Inf   249    275
#>     2      272       7.09 38.4   <0.001   Inf   258    286
#>     3      283       7.71 36.7   <0.001 977.2   268    298
#>     4      293       8.56 34.3   <0.001 853.0   277    310
#>     5      304       9.58 31.7   <0.001 730.3   285    323
#>     6      314      10.73 29.3   <0.001 623.5   293    335
#>     7      325      11.97 27.1   <0.001 535.6   301    348
#>     8      335      13.28 25.2   <0.001 464.6   309    361
#>     9      346      14.63 23.6   <0.001 407.5   317    374
#> 
#> Type:  response 
#> Columns: Days, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
modelbased::estimate_means(m, "Days")
#> Estimated Marginal Means
#> 
#> Days |   Mean |    SE |           95% CI
#> ----------------------------------------
#> 0.00 | 251.41 |  6.82 | [237.01, 265.80]
#> 1.00 | 261.87 |  6.79 | [247.55, 276.19]
#> 2.00 | 272.34 |  7.09 | [257.37, 287.31]
#> 3.00 | 282.81 |  7.71 | [266.55, 299.06]
#> 4.00 | 293.27 |  8.56 | [275.22, 311.32]
#> 5.00 | 303.74 |  9.58 | [283.53, 323.96]
#> 6.00 | 314.21 | 10.73 | [291.57, 336.85]
#> 7.00 | 324.68 | 11.97 | [299.42, 349.94]
#> 8.00 | 335.14 | 13.28 | [307.13, 363.16]
#> 9.00 | 345.61 | 14.63 | [314.75, 376.47]
#> 
#> Marginal means estimated at Days

Created on 2024-12-13 with reprex v2.1.1

@strengejacke
Copy link
Member Author

Not sure why default type "response" returns different results than type = "link" with inverse-link... This here equals to emmeans. I guess the example above really "averaged" across random effects.

data(Salamanders, package = "glmmTMB")

# marginal means - averaged over random effects (not similar to emmeans)
m <- glmmTMB::glmmTMB(count ~ mined + spp + (1 | site), family = poisson(), data = Salamanders)
out <- marginaleffects::avg_predictions(m, by = "mined", type = "link")
modelbased::estimate_means(m, "mined")
#> Estimated Marginal Means
#> 
#> mined | rate |   SE |       95% CI
#> ----------------------------------
#> yes   | 0.18 | 0.04 | [0.12, 0.28]
#> no    | 1.75 | 0.30 | [1.25, 2.46]
#> 
#> Marginal means estimated at mined
insight::link_inverse(m)(out$estimate)
#> [1] 0.188521 1.767239

Created on 2024-12-13 with reprex v2.1.1

@strengejacke
Copy link
Member Author

@mattansb ?

@strengejacke
Copy link
Member Author

@DominiqueMakowski As far as I understood, basically, emmeans and predict return conditional predictions, while marginaleffects returns marginal predictions: https://www.andrewheiss.com/blog/2022/11/29/conditional-marginal-marginaleffects/

it would still be estimate_means(), but marginalized over random effects.

@strengejacke
Copy link
Member Author

And we can then just add hypothesis = "..." for the contrast/comparison function (estimate_contrasts()), and avg_slopes() for the trend of continuous focal terms.

@strengejacke strengejacke self-assigned this Dec 13, 2024
strengejacke added a commit that referenced this issue Dec 13, 2024
@strengejacke
esrtimate_means(backend = "marginaleffects") not calculating margin…
9e47d1b 
…al means

Fixes #273
@strengejacke strengejacke mentioned this issue Dec 13, 2024
Closed
@strengejacke strengejacke changed the title esrtimate_means(backend = "marginaleffects") not calculating marginal means estimate_means(backend = "marginaleffects") not calculating marginal means Dec 13, 2024
strengejacke added a commit that referenced this issue Dec 13, 2024
@strengejacke
estimate_means(backend = "marginaleffects") not calculating margina…
2d202a6 
…l means

Fixes #273
@strengejacke strengejacke linked a pull request Dec 13, 2024 that will close this issue
Open
@mattansb
Copy link
Member

This is Jensen's inequality - the mean of a transformed estimates is not that same as the transformed mean of estimates.

@strengejacke
Copy link
Member Author

I thought bias-correction is only an issue for mixed models?

@mattansb
Copy link
Member

It happens in fixed effect GLMs too

@strengejacke
Copy link
Member Author

Ok, so this is only an issue if we do not use the "link" and "linkinverse", but "response" directly?
image

@mattansb
Copy link
Member

If those are the only 3 options, then yes (:

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@strengejacke strengejacke

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

estimate_means(backend = "marginaleffects") not calculating marginal means easystats/modelbased
3 participants
@DominiqueMakowski @strengejacke @mattansb