lstrends reporting partial results

[joelem]

Hello, thank you for such an amazing package.

I came across an issue that might be an actual issue or maybe me misunderstanding the underlying computations.

I have a multinom() model that predicts probability of being in 10 groups from continuous variable "x". When I try to estimate the slopes per group, lstrends only reports the first 5 groups (groups 6 to 10 are missing). However, when I estimate a marginal mean/probability per group using emmeans from the same model, then it correctly reports all 10 groups. I tried restarting R and not activating other packages in case there were some interferences, but the issue persists. I am on emmeans version 1.4.6

Code & link to example data below:

library(emmeans)
library(nnet)
testdata = read.csv('https://www.dropbox.com/s/n8thorvflolqi5v/testdata.csv?dl=1', header = T)
test = multinom(group ~ x, dat =testdata)
emtrends(test, ~group, var = 'x')
emmeans(test, ~x | group)

[rvlenth]

Yikes. Thanks for the report. Not only is there the wrong number of groups, but the results are really wacko. An approximation to the right answer can be obtained by comparing EMMs at x values that differ by 1:

> confint(pairs(emmeans(test, ~x|group, at = list(x=59:60)), rev = TRUE), by = NULL)
 contrast group  estimate       SE df  lower.CL upper.CL
 60 - 59  1     -0.000570 0.000641 18 -0.001916 0.000776
 60 - 59  2     -0.000895 0.000595 18 -0.002145 0.000354
 60 - 59  3      0.001104 0.000688 18 -0.000341 0.002550
 60 - 59  4      0.000569 0.000476 18 -0.000430 0.001569
 60 - 59  5     -0.000988 0.000663 18 -0.002381 0.000405
 60 - 59  6     -0.000100 0.000601 18 -0.001364 0.001163
 60 - 59  7      0.000516 0.000479 18 -0.000491 0.001523
 60 - 59  8     -0.000513 0.000490 18 -0.001543 0.000518
 60 - 59  9      0.000183 0.000344 18 -0.000540 0.000906
 60 - 59  10     0.000694 0.000242 18  0.000186 0.001203

Confidence level used: 0.95 

I think this must have something to do with the internal bookkeeping. A couple iof versions ago, I changed to code for emtrends() quite a bit to make it more efficient when estimating polynomial trends. It involves a hook in the ref_grid() function, and I think I overlooked the possibility of a multivariate response and what it does with those. Note that your example really messes up if I try to estimate the quadratic trend:

> emtrends(test, ~group, var = 'x')
 group x.trend    SE df lower.CL upper.CL
 1      -0.242 0.182 18   -0.625    0.141
 2      -0.564 0.156 18   -0.892   -0.237
 3      -0.966 0.173 18   -1.330   -0.603
 4      -0.397 0.115 18   -0.639   -0.154
 5      -1.524 0.144 18   -1.826   -1.222

Results are averaged over the levels of: x 
Confidence level used: 0.95 

> emtrends(test, ~group, var = 'x', max.degree = 2)
 Error in emmeans(object = new("emmGrid", model.info = list(call = multinom(formula = group ~  : 
  Irregular reference grid: Marginal means cannot be determined 

My guess is if you had an odd number of levels in the multivariate response, you'd get an error too.

Thanks for the report. I'll post here again when I have some answers.

[rvlenth]

It turns out that I was getting all the right calculations, but the multivariate response disturbed the pattern I was expecting. So I had to replace the code that determined which results to pick out
for each value of h in computing values of EMM(x+h) - EMM(x).

> emtrends(test, ~group, var = 'x')
 group   x.trend       SE df  lower.CL upper.CL
 1     -0.000571 0.000640 18 -0.001916 0.000775
 2     -0.000895 0.000594 18 -0.002144 0.000353
 3      0.001104 0.000689 18 -0.000343 0.002551
 4      0.000569 0.000476 18 -0.000431 0.001570
 5     -0.000988 0.000663 18 -0.002380 0.000404
 6     -0.000101 0.000601 18 -0.001365 0.001163
 7      0.000516 0.000480 18 -0.000492 0.001524
 8     -0.000513 0.000490 18 -0.001542 0.000517
 9      0.000183 0.000344 18 -0.000540 0.000906
 10     0.000696 0.000243 18  0.000185 0.001207

Confidence level used: 0.95 

I am pleased to see that these happen to be close to the values I computed roughly earlier.

I will push the revised code up soon.

[joelem]

Awesome, thank you so much for the quick resolution & again for the package!

[rvlenth]

You can do:

remotes::install_github("rvlenth/emmeans", dependencies = TRUE, build_opts = "")

to install it from github. If using Windows, you need to have installed the latest version of Rtools as well.

I think this is resolved, so will close this issue. But you can comment further or even reopen it if something comes up. Thanks for your report.

[rvlenth]

I'm going to go ahead and close this issue as I think we have it fixed. The current CRAN version of emmeans includes this update.

[rvlenth]

Just alerting you to an embarrassing error in emmeans's support of multinomial models (class "multinom"). In those objects store the regression coefficients as the transpose to how they are stored in most multivariate models. My code accounts for that, but it did not account for the fact that this alters the order of rows and columns in the variance-covariance matrix of these coefficients. Consequently, emmeans versions up through 1.5.3 computed incorrect SEs of estimates produced by emmeans() or related functions. (The estimates themselves are correct, just not the standard errors.) The code is corrected in version 1.5.4, just submitted to CRAN.

My apologies for any difficulties this has created in your past analyses or publications.

Sincerely,

Russ Lenth

[rvlenth]

Note: Please see Issue #253 relating to multinomial models. A new version 1.5.4 has been submitted to CRAN that corrects this error.