search

rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
358 stars 31 forks source link

nnet and emmeans odds ratio discrepancy #404

Closed jgeller112 closed 1 year ago

jgeller112 commented 1 year ago

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?

rvlenth commented 1 year ago

Please show me the results of this code.

Russ

jgeller112 commented 1 year ago

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
rvlenth commented 1 year ago

So what is discrepant? Which of the 28 odds ratios am I supposed to compare with the two model parameters?

jgeller112 commented 1 year ago

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
rvlenth commented 1 year ago

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.

rvlenth commented 1 year ago

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.
jgeller112 commented 1 year ago

Thanks for taking the time to do all this. This all makes sense to me. I will close the thread.

rvlenth commented 1 year ago

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.