emmeans() for {ordbetareg} models only reports output on the log odds scale
mz555 opened this issue · 6 comments
I fit a model using the {ordbetareg} package.
ordbeta.m1 <- ordbetareg(
bf(Count ~ condition * Profession + (1 | part_id)),
data = MAGL.data.long,
true_bounds = c(0, 3))
When I try to compute pairwise contrasts for the interaction term of interest, emmeans() will only provide responses on the log odds scale, regardless of what the type = argument is. This was confirmed by the {ggeffects} package dev who suggested I make this ticket.
emm.RQ1 <- emmeans(ordbeta.m1, ~condition | Profession, type = "response")) |>
contrast("trt.vs.ctrl", ref = 1)
Output:
Profession = writ:
condition emmean lower.HPD upper.HPD
1 -0.3116 -0.3683 -0.2564
2 -0.3609 -0.4209 -0.3041
3 -0.0280 -0.0880 0.0298
4 -0.2011 -0.2547 -0.1491
5 -0.1917 -0.2413 -0.1359
Point estimate displayed: median
HPD interval probability: 0.95
And changed to something else, type = "scale":
Profession = writ:
condition emmean lower.HPD upper.HPD
1 -0.3116 -0.3683 -0.2564
2 -0.3609 -0.4209 -0.3041
3 -0.0280 -0.0880 0.0298
4 -0.2011 -0.2547 -0.1491
5 -0.1917 -0.2413 -0.1359
Point estimate displayed: median
HPD interval probability: 0.95
Describe the bug
I would like for the output to either be on the DV scale (counts here) or the probability scale (0-1), but this only is possible in the marginaleffects package for now.
Thank you for the quick response. I have also contacted the {ordbetareg} dev. For context, their package is just a wrapper and some special families for the {brms} package, which i believe is natively supported by {emmeans} and does work well with beta family models.
Just to add: the same issue can be seen for glmmTMB(family=ordbeta()). So it might be emmeans not detecting the response / link inverse correctly? Just a guess
Please provide a reproducible example. I don't have the MAGL.data.long data
Here's one. But here it looks like glmmTMB gets it right? ordbetareg just wraps brms.
library(ordbetareg)
#> Loading required package: brms
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.21.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#>
#> ar
library(glmmTMB)
#>
#> Attaching package: 'glmmTMB'
#> The following object is masked from 'package:brms':
#>
#> lognormal
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
data("pew")
model_data <- select(pew, therm,
education = "F_EDUCCAT2_FINAL",
region = "F_CREGION_FINAL",
income = "F_INCOME_FINAL"
)
model_data$therm2 <- datawizard::normalize(model_data$therm)
m <- ordbetareg(
formula = therm ~ education + income + (1 | region),
data = model_data,
cores = 2, chains = 2, backend = "rstan"
)
#> Normalizing using the observed bounds of 0 - 100. If these are incorrect, please pass the bounds to use to the true_bounds parameter.
#> Warning: Rows containing NAs were excluded from the model.
#> Compiling Stan program...
#> Start sampling
#> Warning: There were 2 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
m2 <- glmmTMB(
formula = therm2 ~ education + income + (1 | region),
data = model_data,
family = ordbeta()
)
emmeans::emmeans(m, "income", regrid = "mu")
#> income emmean lower.HPD upper.HPD
#> Less than $10,000 0.3653 0.1158 0.653
#> 10 to under $20,000 0.2769 0.0382 0.502
#> 20 to under $30,000 0.4522 0.2470 0.702
#> 30 to under $40,000 0.2851 0.0811 0.537
#> 40 to under $50,000 0.2311 0.0132 0.466
#> 50 to under $75,000 0.2145 0.0150 0.430
#> 75 to under $100,000 0.1021 -0.0839 0.310
#> 100 to under $150,000 [OR] 0.0881 -0.1216 0.288
#> $150,000 or more 0.1106 -0.0943 0.329
#> (VOL) Don't know/Refused -0.1048 -0.3589 0.142
#>
#> Results are averaged over the levels of: education
#> Point estimate displayed: median
#> HPD interval probability: 0.95
emmeans::emmeans(m, "income", regrid = "response")
#> income emmean lower.HPD upper.HPD
#> Less than $10,000 0.3653 0.1158 0.653
#> 10 to under $20,000 0.2769 0.0382 0.502
#> 20 to under $30,000 0.4522 0.2470 0.702
#> 30 to under $40,000 0.2851 0.0811 0.537
#> 40 to under $50,000 0.2311 0.0132 0.466
#> 50 to under $75,000 0.2145 0.0150 0.430
#> 75 to under $100,000 0.1021 -0.0839 0.310
#> 100 to under $150,000 [OR] 0.0881 -0.1216 0.288
#> $150,000 or more 0.1106 -0.0943 0.329
#> (VOL) Don't know/Refused -0.1048 -0.3589 0.142
#>
#> Results are averaged over the levels of: education
#> Point estimate displayed: median
#> HPD interval probability: 0.95
emmeans::emmeans(m2, "income", regrid = "mu")
#> income response SE df asymp.LCL asymp.UCL
#> Less than $10,000 0.590 0.0279 Inf 0.535 0.644
#> 10 to under $20,000 0.566 0.0236 Inf 0.520 0.613
#> 20 to under $30,000 0.609 0.0231 Inf 0.564 0.654
#> 30 to under $40,000 0.569 0.0229 Inf 0.525 0.614
#> 40 to under $50,000 0.556 0.0233 Inf 0.510 0.602
#> 50 to under $75,000 0.552 0.0197 Inf 0.514 0.591
#> 75 to under $100,000 0.526 0.0197 Inf 0.487 0.564
#> 100 to under $150,000 [OR] 0.521 0.0201 Inf 0.482 0.561
#> $150,000 or more 0.527 0.0209 Inf 0.486 0.568
#> (VOL) Don't know/Refused 0.473 0.0273 Inf 0.419 0.526
#>
#> Results are averaged over the levels of: education
#> Confidence level used: 0.95
emmeans::emmeans(m2, "income", regrid = "response")
#> income response SE df asymp.LCL asymp.UCL
#> Less than $10,000 0.590 0.0279 Inf 0.535 0.644
#> 10 to under $20,000 0.566 0.0236 Inf 0.520 0.613
#> 20 to under $30,000 0.609 0.0231 Inf 0.564 0.654
#> 30 to under $40,000 0.569 0.0229 Inf 0.525 0.614
#> 40 to under $50,000 0.556 0.0233 Inf 0.510 0.602
#> 50 to under $75,000 0.552 0.0197 Inf 0.514 0.591
#> 75 to under $100,000 0.526 0.0197 Inf 0.487 0.564
#> 100 to under $150,000 [OR] 0.521 0.0201 Inf 0.482 0.561
#> $150,000 or more 0.527 0.0209 Inf 0.486 0.568
#> (VOL) Don't know/Refused 0.473 0.0273 Inf 0.419 0.526
#>
#> Results are averaged over the levels of: education
#> Confidence level used: 0.95Created on 2024-04-29 with reprex v2.1.0
This issue has been languishing here, but I honestly don't know what to do with it, so I am closing it. emmeans does not support the ordbetareg package, and if that package's code is just a wrapper for brms code, that's still not an emmeans issue because brms provides its own emmeans support. There are several options for that brms code. I suggest you look at ? emm_basis.brmsfit.