rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
341 stars 30 forks source link

mvcontrast returns list, not emm_list #410

Closed WillemSleegers closed 1 year ago

WillemSleegers commented 1 year ago

Describe the bug

mvcontrast() sometimes returns a list and not an emm_list.

To reproduce

moats_lm <- lm(yield ~ Variety + Block, data = MOats)
moats_emm <- emmeans(moats_lm, ~ Variety | rep.meas)
moats_mvcontrast <- mvcontrast(moats_emm, "consec", show.ests = TRUE)
class(moats_mvcontrast)

Expected behavior

class(moats_mvcontrast) should return an object with class emm_list.

Additional context

I'm the developer of tidystats, which takes statistical objects, extracts the statistics, and re-organizes them to store and report them. I'm working on adding support for emmeans and the object not returning an object of class emm_list makes things a bit tricky.

rvlenth commented 1 year ago

Hmmm, that doesn't sound good. I'll try to look at this soon.

rvlenth commented 1 year ago

Looking at this more closely, this is what we have:

> class(moats_mvcontrast)
[1] "list"

> class(moats_mvcontrast[[1]])
[1] "emmGrid"
attr(,"package")
[1] "emmeans"

> class(moats_mvcontrast[[2]])
[1] "summary_emm" "data.frame" 

So this cannot be of class emm_list. The call with show.ests = TRUE is in fact asking for two different sets of output, and in a programmatic setting, I think it is best to do them separately anyway:

> mvc <- mvcontrast(moats_emm, "consec")
> mvc
 contrast                 T.square df1 df2 F.ratio p.value
 Marvellous - Golden Rain    3.082   4   7   0.539  0.9174
 Victory - Marvellous        9.181   4   7   1.607  0.4726

P value adjustment: sidak 

> mvc_est <- summary(contrast(moats_emm, "consec"), by = "contrast")
> mvc_est
contrast = Marvellous - Golden Rain:
 rep.meas estimate    SE df t.ratio p.value
 0            6.67  7.84 10   0.851  0.8169
 0.2         10.00  9.34 10   1.071  0.6859
 0.4          2.50 12.30 10   0.203  0.9985
 0.6          2.00 10.33 10   0.194  0.9988

contrast = Victory - Marvellous:
 rep.meas estimate    SE df t.ratio p.value
 0          -15.17  7.84 10  -1.936  0.2354
 0.2        -18.83  9.34 10  -2.017  0.2083
 0.4         -6.33 12.30 10  -0.515  0.9580
 0.6         -8.33 10.33 10  -0.807  0.8403

Results are averaged over the levels of: Block 
P value adjustment: mvt method for 4 tests 

As shown, mvc_est is of class summary_emm. If you want an emmGrid instead, use

> mvc_emm <- update(contrast(moats_emm, "consec"), by = "contrast")

FWIW, the documentation does say, under Value:

An object of class summary_emm containing the multivariate test results; or a list of the estimates and the tests if show.ests is TRUE.

So I will claim that this is working as documented

WillemSleegers commented 1 year ago

Ah clear! My apologies.

rvlenth commented 1 year ago

I think this is cleared-up, so am closing