Odds ratios to interpret a multinomial regression model

[ChrisBilder]

Hello and thank you for creating this package. I found a small problem with the calculation of odds ratios to interpret explanatory variables in a multinomial regression model. I have verified that version 1.8.8 of emmeans does the calculation correctly, but version 2.0.2 does not. The problem is simply the newer version does not apply the exponential function to the difference of two log odds.

Below is code that estimates a multinomial regression model to the classic pneumoconiosis data (see the pneumo data frame of the VGAM package and exercise 23 of Bilder and Loughin, 2024). A modified format version of the data set is at chrisbilder.com/stat875/graded/project4.csv. The purpose is to use time (years in coal mine) to estimate resp (normal, mild, or severe pneumoconiosis). Version 1.8.8 estimates the odds ratios for time correctly, such as 0.658 for the mild vs. normal comparison with a 5-year reduction in time. Version 2.0.2 returns the difference in the log odds for mild vs. normal. While one can simply use exp() to obtain the odds ratio, this was not needed with earlier versions of the package.

Thank you,

Chris

project4 <- read.csv(file = 'https://chrisbilder.com/stat875/graded/project4.csv')
head(project4, n = 3)

library(nnet)

mod.fit.nom <- multinom(formula = resp ~ time, data = project4, weights = count)
summary(mod.fit.nom)

library(emmeans)
calc.est <- emmeans(object = mod.fit.nom, specs = ~ resp + time, mode = "latent",  at = list(time = c(15, 20)))
test.info <- contrast(object = calc.est, interaction = list("pairwise", "pairwise"))
confint(object = test.info, df = Inf, type = "response", level = 0.95)

[rvlenth]

I've been able to verify this. Evidently, somehow the log transformation got dropped from the object.
The short-term fix is to put it back in:

> test.info <- update(test.info, tran = "log")
> confint(object = test.info, df = Inf, type = "response", level = 0.95)
 resp_pairwise   time_pairwise response     SE  df asymp.LCL asymp.UCL
 mild - normal   15 - 20          0.658 0.0503 Inf     0.567     0.765
 mild - severe   15 - 20          1.137 0.1120 Inf     0.937     1.380
 normal - severe 15 - 20          1.727 0.1420 Inf     1.470     2.030

Degrees-of-freedom method: user-specified 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

Thanks for reporting this. It looks easy to fix; but more important, we'll try to understand why this happened.

[rvlenth]

Here's what happened: With version 1.10.4, we introduced an mvregrid() function that handles multivariate transformations. At that time the transformation was changed from "log" to "clr" (centered log ratio, in the compositions package). So now we see, with your objects:

> summary(calc.est, type = "response")
 resp   time response    SE df lower.CL upper.CL
 mild     15   -0.885 0.235  4    -1.54   -0.233
 normal   15    2.153 0.170  4     1.68    2.624
 severe   15   -1.268 0.259  4    -1.99   -0.550
 mild     20   -0.789 0.193  4    -1.32   -0.254
 normal   20    1.832 0.137  4     1.45    2.213
 severe   20   -1.043 0.211  4    -1.63   -0.456

Unknown transformation "clr": no transformation done 
Confidence level used: 0.95 
Warning message:
In summary.emmGrid(calc.est, type = "response") :
  It looks like you have a multivariate transformation -- not handled by type = 'response'.
Perhaps instead you wanted to use 'mvregrid()'?

You have stumbled upon what I now see is an incomplete implementation of that clr multivariate transformation. One aspect is that "clr" as not recognized as a log-family transformation in contrast(), and that's why we don't back-transform it correctly.

I've gone ahead and modified the contrast() and multinom support code to solve the latter problem; so now it recognizes alr, clr, and ilr transformations as log-ratio transformations:

> test.info <- contrast(object = calc.est, interaction = list("pairwise", "pairwise"))
> confint(object = test.info, df = Inf, type = "response", level = 0.95)
 resp_pairwise   time_pairwise odds.ratio     SE  df asymp.LCL asymp.UCL
 mild / normal   15 / 20            0.658 0.0503 Inf     0.567     0.765
 mild / severe   15 / 20            1.137 0.1120 Inf     0.937     1.380
 normal / severe 15 / 20            1.727 0.1420 Inf     1.470     2.030

Degrees-of-freedom method: user-specified 
Confidence level used: 0.95 
Intervals are back-transformed from the log odds ratio scale 

I'm considering whether in the example above, I should just do the multivariate reverse transformation, instead opf just warning about it.

In any case, you'll see a notice when I push an update (probably to the devel branch) that references this issue.

[rvlenth]

@ChrisBilder -- you should be able to install the modified code now via

remotes::install_github("rvlenth/emmeans@devel", dependencies = TRUE, build_vignettes = TRUE)

--- and it should work as shown at the end of the preceding comment.

These changes will eventually be incorporated in the next CRAN update.