rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
364 stars 32 forks source link

Failures to disable Tukey adjustment when it should #275

Closed rvlenth closed 3 years ago

rvlenth commented 3 years ago

Here are two situations where Tukey adjustments are allowed, but shouldn't be.

Case A

library(emmeans)

warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)
EMM = emmeans(warp.lm, ~ tension | wool)

test(pairs(EMM), by = NULL)
##  contrast wool estimate   SE df t.ratio p.value
##  L - M    A      20.556 5.16 48  3.986  0.0007 
##  L - H    A      20.000 5.16 48  3.878  0.0009 
##  M - H    A      -0.556 5.16 48 -0.108  0.9936 
##  L - M    B      -0.556 5.16 48 -0.108  0.9936 
##  L - H    B       9.444 5.16 48  1.831  0.1704 
##  M - H    B      10.000 5.16 48  1.939  0.1389 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Here, we no longer have a family of pairwise comparisons (we have two such families combined into one). So the Tukey adjustment shouldn't be used. Worse, it defaults to Tukey. Even worse, it remembers the adjustment for one group, and uses that adjustment for all of them combined -- making this adjustment extremely liberal.

Case B

Continuing with the same objects...

ALL = rbind(pairs(EMM, simple = "each"))

test(ALL, adjust = "tukey")
##  wool tension contrast estimate   SE df t.ratio p.value
##  A    .       L - M      20.556 5.16 48  3.986  0.0018 
##  A    .       L - H      20.000 5.16 48  3.878  0.0026 
##  A    .       M - H      -0.556 5.16 48 -0.108  0.9999 
##  B    .       L - M      -0.556 5.16 48 -0.108  0.9999 
##  B    .       L - H       9.444 5.16 48  1.831  0.3467 
##  B    .       M - H      10.000 5.16 48  1.939  0.2922 
##  .    L       A - B      16.333 5.16 48  3.167  0.0196 
##  .    M       A - B      -4.778 5.16 48 -0.926  0.8683 
##  .    H       A - B       5.778 5.16 48  1.120  0.7727 
## 
## P value adjustment: tukey method for comparing a family of 4.772 estimates

test(ALL, adjust = "mvt")
##  wool tension contrast estimate   SE df t.ratio p.value
##  A    .       L - M      20.556 5.16 48  3.986  0.0020 
##  A    .       L - H      20.000 5.16 48  3.878  0.0026 
##  A    .       M - H      -0.556 5.16 48 -0.108  1.0000 
##  B    .       L - M      -0.556 5.16 48 -0.108  1.0000 
##  B    .       L - H       9.444 5.16 48  1.831  0.3793 
##  B    .       M - H      10.000 5.16 48  1.939  0.3192 
##  .    L       A - B      16.333 5.16 48  3.167  0.0206 
##  .    M       A - B      -4.778 5.16 48 -0.926  0.9119 
##  .    H       A - B       5.778 5.16 48  1.120  0.8258 
## 
## P value adjustment: mvt method for 9 tests

Created on 2021-05-03 by the reprex package (v1.0.0) In this situation, at least the default isn't Tukey (it is Bonferroni, not shown). But it allows the Tukey adjustment, which it shouldn't. And note that it is a bit liberal: the P values are smaller than for "mvt" which is the exact method.

These errors happen because of remembered attributes leftover from the pairs() calls, and failure to forget them when they don't apply. Case B is pretty easily fixed because we should just always have rbind() disable pairwise. Case A is trickier because we need to continue to allow Tukey when the by groups are what they are supposed to be.

rvlenth commented 3 years ago

OK, I think we've got it. Basic approach was to remember the by variables used whenever we construct paired comparisons (misc$estType = "pairs"), and make sure those haven't changed. In the rbind case, we changed misc$estType to "rbind" so it is never going to allow Tukey adjustments.

library(emmeans)

warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)
EMM = emmeans(warp.lm, ~ tension | wool)

test(pairs(EMM), by = NULL)
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
##  contrast wool estimate   SE df t.ratio p.value
##  L - M    A      20.556 5.16 48  3.986  0.0014 
##  L - H    A      20.000 5.16 48  3.878  0.0019 
##  M - H    A      -0.556 5.16 48 -0.108  1.0000 
##  L - M    B      -0.556 5.16 48 -0.108  1.0000 
##  L - H    B       9.444 5.16 48  1.831  0.3665 
##  M - H    B      10.000 5.16 48  1.939  0.3030 
## 
## P value adjustment: sidak method for 6 tests

ALL = rbind(pairs(EMM, simple = "each"))
test(ALL, adjust = "tukey")
## Note: adjust = "tukey" was changed to "sidak"
## because "tukey" is only appropriate for one set of pairwise comparisons
##  wool tension contrast estimate   SE df t.ratio p.value
##  A    .       L - M      20.556 5.16 48  3.986  0.0021 
##  A    .       L - H      20.000 5.16 48  3.878  0.0029 
##  A    .       M - H      -0.556 5.16 48 -0.108  1.0000 
##  B    .       L - M      -0.556 5.16 48 -0.108  1.0000 
##  B    .       L - H       9.444 5.16 48  1.831  0.4958 
##  B    .       M - H      10.000 5.16 48  1.939  0.4181 
##  .    L       A - B      16.333 5.16 48  3.167  0.0238 
##  .    M       A - B      -4.778 5.16 48 -0.926  0.9817 
##  .    H       A - B       5.778 5.16 48  1.120  0.9398 
## 
## P value adjustment: sidak method for 9 tests

Created on 2021-05-03 by the reprex package (v1.0.0)

I also changed famSize returned by rbind to the number of estimates. Apparently I once allowed -- as in Case B above -- for a Tukey adjustment in rbind() results. Or maybe it is an artifact of the olden days when I also used famSize in the Scheffe adjustment, but now I compute the dimension of the estimates.

rvlenth commented 3 years ago

I think this is now resolved, so am closing this issue.

SchmidtPaul commented 2 years ago

While I think I understand your reasoning above and also read #91 as well as your comment here, I am irritated why I get neither Tukey-p-values, nor Tukey-CIs in the case below. Isn't it just three simple groups and thus a set of pairwise comparisons?

library(emmeans)

mod <- lm(weight ~ group, data = PlantGrowth)
emm <- emmeans(mod, ~ group, adjust = "tukey")
con <- contrast(emm, adjust = "tukey")

emm
#> Note: adjust = "tukey" was changed to "sidak"
#> because "tukey" is only appropriate for one set of pairwise comparisons
#>  group emmean    SE df lower.CL upper.CL
#>  ctrl    5.03 0.197 27     4.53     5.53
#>  trt1    4.66 0.197 27     4.16     5.16
#>  trt2    5.53 0.197 27     5.02     6.03
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 3 estimates
con
#> Note: adjust = "tukey" was changed to "sidak"
#> because "tukey" is only appropriate for one set of pairwise comparisons
#>  contrast    estimate    SE df t.ratio p.value
#>  ctrl effect   -0.041 0.161 27  -0.255  0.9921
#>  trt1 effect   -0.412 0.161 27  -2.560  0.0484
#>  trt2 effect    0.453 0.161 27   2.814  0.0268
#> 
#> P value adjustment: sidak method for 3 tests

Created on 2022-07-14 by the reprex package (v2.0.1)

rvlenth commented 2 years ago

You did effect contrasts (the default), not pairwise comparisons. Use

con <- contrast(emm, "pairwise")

for which the tukey adjustment is the default.

Sent from my iPad

On Jul 14, 2022, at 3:55 AM, Paul Schmidt @.***> wrote:

con <- contrast(emm, adjust = "tukey")

SchmidtPaul commented 2 years ago

Ah yes, you are right. Sorry and thanks. (corrected reprex below)

But I am still not sure why I then cannot get the corresponding Tukey-CLs for the emmeans.

library(emmeans)

mod <- lm(weight ~ group, data = PlantGrowth)
emm <- emmeans(mod, ~ group, adjust = "tukey")
con <- contrast(emm, method = "pairwise", adjust = "tukey")

emm
#> Note: adjust = "tukey" was changed to "sidak"
#> because "tukey" is only appropriate for one set of pairwise comparisons
#>  group emmean    SE df lower.CL upper.CL
#>  ctrl    5.03 0.197 27     4.53     5.53
#>  trt1    4.66 0.197 27     4.16     5.16
#>  trt2    5.53 0.197 27     5.02     6.03
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 3 estimates
con
#>  contrast    estimate    SE df t.ratio p.value
#>  ctrl - trt1    0.371 0.279 27   1.331  0.3909
#>  ctrl - trt2   -0.494 0.279 27  -1.772  0.1980
#>  trt1 - trt2   -0.865 0.279 27  -3.103  0.0120
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates

Created on 2022-07-14 by the reprex package (v2.0.1)

rvlenth commented 2 years ago

because the family of means is not a family of pairwise comparisons

Sent from my iPad

On Jul 14, 2022, at 8:17 AM, Paul Schmidt @.***> wrote:



Ah yes, you are right. Sorry and thanks. (corrected reprex below)

But I am still not sure why I then cannot get the corresponding Tukey-CLs for the emmeans.

library(emmeans)

mod <- lm(weight ~ group, data = PlantGrowth) emm <- emmeans(mod, ~ group, adjust = "tukey") con <- contrast(emm, method = "pairwise", adjust = "tukey")

emm

> Note: adjust = "tukey" was changed to "sidak"

> because "tukey" is only appropriate for one set of pairwise comparisons

> group emmean SE df lower.CL upper.CL

> ctrl 5.03 0.197 27 4.53 5.53

> trt1 4.66 0.197 27 4.16 5.16

> trt2 5.53 0.197 27 5.02 6.03

>

> Confidence level used: 0.95

> Conf-level adjustment: sidak method for 3 estimates

con

> contrast estimate SE df t.ratio p.value

> ctrl - trt1 0.371 0.279 27 1.331 0.3909

> ctrl - trt2 -0.494 0.279 27 -1.772 0.1980

> trt1 - trt2 -0.865 0.279 27 -3.103 0.0120

>

> P value adjustment: tukey method for comparing a family of 3 estimates

Created on 2022-07-14 by the reprex packagehttps://reprex.tidyverse.org (v2.0.1)

— Reply to this email directly, view it on GitHubhttps://github.com/rvlenth/emmeans/issues/275#issuecomment-1184437973, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGMJPL5LVFY52DPXZZ3ENR3VUAHPRANCNFSM44BOT6PQ. You are receiving this because you modified the open/close state.Message ID: @.***>

rvlenth commented 2 years ago

You could do confint(con) to get Tukey-adjusted confidence intervals for the pairwise differences. Is that what you were wanting?

Sent from my iPad

On Jul 14, 2022, at 8:17 AM, Paul Schmidt @.***> wrote:



Ah yes, you are right. Sorry and thanks. (corrected reprex below)

But I am still not sure why I then cannot get the corresponding Tukey-CLs for the emmeans.

library(emmeans)

mod <- lm(weight ~ group, data = PlantGrowth) emm <- emmeans(mod, ~ group, adjust = "tukey") con <- contrast(emm, method = "pairwise", adjust = "tukey")

emm

> Note: adjust = "tukey" was changed to "sidak"

> because "tukey" is only appropriate for one set of pairwise comparisons

> group emmean SE df lower.CL upper.CL

> ctrl 5.03 0.197 27 4.53 5.53

> trt1 4.66 0.197 27 4.16 5.16

> trt2 5.53 0.197 27 5.02 6.03

>

> Confidence level used: 0.95

> Conf-level adjustment: sidak method for 3 estimates

con

> contrast estimate SE df t.ratio p.value

> ctrl - trt1 0.371 0.279 27 1.331 0.3909

> ctrl - trt2 -0.494 0.279 27 -1.772 0.1980

> trt1 - trt2 -0.865 0.279 27 -3.103 0.0120

>

> P value adjustment: tukey method for comparing a family of 3 estimates

Created on 2022-07-14 by the reprex packagehttps://reprex.tidyverse.org (v2.0.1)

— Reply to this email directly, view it on GitHubhttps://github.com/rvlenth/emmeans/issues/275#issuecomment-1184437973, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGMJPL5LVFY52DPXZZ3ENR3VUAHPRANCNFSM44BOT6PQ. You are receiving this because you modified the open/close state.Message ID: @.***>