nnet and emmeans odds ratio discrepancy
Closed this issue · 8 comments
Model parameters of the odds ratio taken from the model differ from what emmeans shows
To reproduce
df <- tibble(mcType = factor(rep(c("licit", "illicit"),
times = c(534,1207))),
cond = factor(c(rep(c("pain","mhsu","allOther","sleep"),
times = c(280,141,82,31)),
rep(c("pain","mhsu","allOther","sleep"),
times = c(491,360,208,148))),
levels = c("pain","sleep","mhsu","allOther")))
mm <- multinom(cond ~ mcType,
data = df))
library(modelbased)
model_parameters(mm, exponentiate = TRUE)
This differs from estimates using emmeans:
multi_an <- emmeans(mm, ~ mcType*licit)
d=pairs(regrid(multi_an, transform = "logit"), type = "response", reverse=TRUE)
Expected behavior
I would expect the odds ratios calculated from nnet to be very similar to those produced by emmeans but they are not. Am I missing something critical here?
Sure. This is just for the allother vs. pain comparisons
model_parameters(mm, exponentiate = TRUE)
# Response level: allother
Parameter | Odds Ratio | SE | 95% CI | t(1735) | p
--------------------------------------------------------------------
(Intercept) | 0.42 | 0.04 | [0.36, 0.50] | -10.38 | < .001
mcType [licit] | 0.69 | 0.10 | [0.51, 0.93] | -2.46 | 0.014
multi_an <- emmeans(mm, ~ mcType*licit)
d=pairs(regrid(multi_an, transform = "logit"), type = "response", reverse=TRUE)
contrast odds.ratio SE df null
licit pain / illicit pain 1.6075 0.1681 6 1
illicit sleep / illicit pain 0.2038 0.0244 6 1
illicit sleep / licit pain 0.1268 0.0156 6 1
licit sleep / illicit pain 0.0899 0.0174 6 1
licit sleep / licit pain 0.0559 0.0125 6 1
licit sleep / illicit sleep 0.4410 0.0903 6 1
illicit mhsu / illicit pain 0.6198 0.0661 6 1
illicit mhsu / licit pain 0.3856 0.0413 6 1
illicit mhsu / illicit sleep 3.0413 0.3643 6 1
illicit mhsu / licit sleep 6.8961 1.3479 6 1
licit mhsu / illicit pain 0.5232 0.0598 6 1
licit mhsu / licit pain 0.3255 0.0543 6 1
licit mhsu / illicit sleep 2.5674 0.3380 6 1
licit mhsu / licit sleep 5.8214 1.2923 6 1
licit mhsu / illicit mhsu 0.8442 0.0984 6 1
illicit allOther / illicit pain 0.3036 0.0341 6 1
illicit allOther / licit pain 0.1889 0.0218 6 1
illicit allOther / illicit sleep 1.4898 0.1872 6 1
illicit allOther / licit sleep 3.3780 0.6761 6 1
illicit allOther / illicit mhsu 0.4898 0.0550 6 1
illicit allOther / licit mhsu 0.5803 0.0721 6 1
licit allOther / illicit pain 0.2645 0.0353 6 1
licit allOther / licit pain 0.1646 0.0291 6 1
licit allOther / illicit sleep 1.2980 0.1930 6 1
licit allOther / licit sleep 2.9432 0.6798 6 1
licit allOther / illicit mhsu 0.4268 0.0578 6 1
licit allOther / licit mhsu 0.5056 0.0877 6 1
licit allOther / illicit allOther 0.8713 0.1239 6 1
t.ratio p.value
4.538 0.0407
-13.293 0.0001
-16.747 <.0001
-12.412 0.0002
-12.883 0.0002
-3.997 0.0700
-4.486 0.0428
-8.900 0.0013
9.285 0.0011
9.879 0.0007
-5.666 0.0143
-6.727 0.0060
7.161 0.0043
7.935 0.0025
-1.453 0.8081
-10.612 0.0005
-14.443 0.0001
3.172 0.1669
6.082 0.0101
-6.353 0.0081
-4.379 0.0476
-9.956 0.0007
-10.213 0.0006
1.754 0.6642
4.674 0.0356
-6.283 0.0085
-3.934 0.0747
-0.969 0.9639
Let's take the intercept for that comparison which is illicit allOther / illicit pain
Odd's ratio from. model_parameters
(Intercept) | 0.42
In emmeans
illicit allOther / illicit pain 0.3036
But why should those be the same, and how do you know it's those two. I did not find a function model_parameters(), it must not be in that modelbased package. But this is a model for a two-level factor for a 4-level multinomial response. I get
> coef(mm)
(Intercept) mcTypelicit
sleep -1.1992431 -1.0014884
mhsu -0.3103369 -0.3756443
allOther -0.8589398 -0.3691759
Each of these are on a scale of offsets from the pain condition being zero, and they are on the log scale, not the logit scale.
Here are the estimated probabilities:
> EMM = emmeans(mm, ~mcType|cond)
> EMM
cond = pain:
mcType prob SE df lower.CL upper.CL
illicit 0.4068 0.01414 6 0.3722 0.4414
licit 0.5243 0.02161 6 0.4715 0.5772
cond = sleep:
mcType prob SE df lower.CL upper.CL
illicit 0.1226 0.00944 6 0.0995 0.1457
licit 0.0581 0.01012 6 0.0333 0.0828
cond = mhsu:
mcType prob SE df lower.CL upper.CL
illicit 0.2983 0.01317 6 0.2660 0.3305
licit 0.2641 0.01908 6 0.2174 0.3107
cond = allOther:
mcType prob SE df lower.CL upper.CL
illicit 0.1723 0.01087 6 0.1457 0.1989
licit 0.1535 0.01560 6 0.1154 0.1917
Confidence level used: 0.95
So if you're looking to relate emmeans() results to model coefficients, look at ratios (not odds rations) of probabilities compared with pain. For example, with the 1st row of coefficients...
> c(parm = exp(-1.1992431), prob.ratio = .1226/.4068)
parm prob.ratio
0.3014223 0.3013766
c(parm = exp(-1.1992431 - 1.0014884), prob.ratio = .0581/.5243)
parm prob.ratio
0.1107221 0.1108144
There is no bug here.
Here's the whole shebang:. remember the coefficients:
> coef(mm)
(Intercept) mcTypelicit
sleep -1.1992431 -1.0014884
mhsu -0.3103369 -0.3756443
allOther -0.8589398 -0.3691759
First, the intercepts, in the "illicit" entries:
> coefs = contrast(regrid(EMM, "log"), "trt.vs.ctrl1", by = "mcType")
> update(coefs, by = "contrast")
contrast = sleep - pain:
mcType estimate SE df t.ratio p.value
illicit -1.199 0.0938 6 -12.789 <.0001
licit -2.201 0.1893 6 -11.627 <.0001
contrast = mhsu - pain:
mcType estimate SE df t.ratio p.value
illicit -0.310 0.0694 6 -4.473 0.0079
licit -0.686 0.1033 6 -6.643 0.0011
contrast = allOther - pain:
mcType estimate SE df t.ratio p.value
illicit -0.859 0.0827 6 -10.382 0.0001
licit -1.228 0.1256 6 -9.781 0.0001
Results are given on the log (not the response) scale.
P value adjustment: dunnettx method for 2 tests
Then the mcTypelicit coefficients:
> contrast(coefs, "revpairwise", by = "contrast")
contrast = sleep - pain:
contrast1 estimate SE df t.ratio p.value
licit - illicit -1.001 0.211 6 -4.741 0.0032
contrast = mhsu - pain:
contrast1 estimate SE df t.ratio p.value
licit - illicit -0.376 0.124 6 -3.019 0.0234
contrast = allOther - pain:
contrast1 estimate SE df t.ratio p.value
licit - illicit -0.369 0.150 6 -2.455 0.0494
Results are given on the log (not the response) scale.
Thanks for taking the time to do all this. This all makes sense to me. I will close the thread.
It's pretty easy to think the coefficients in these models are on the logit scale. And in a way they are, if you kind of look at them sideways.