emmeans() for {ordbetareg} models only reports output on the log odds scale

[mz555]

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.

[rvlenth]

Traveling. It will be a while but offhand, I think the problem is just that these objects are not supported in emmeans. The developer of ordbetareg could write the needed methods. there is a vignette on that.

[mz555]

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.

[strengejacke]

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

[rvlenth]

Please provide a reproducible example. I don't have the MAGL.data.long data

[strengejacke]

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.95

^{Created on 2024-04-29 with reprex v2.1.0}

[rvlenth]

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.