lsmeans (R): Adjust for multiple comparisons with interaction terms

Question

I have a lsmeans problem in R. I want to do a post-hoc analysis of an interaction, similar to examples provided in the lsmeans documentation.

I am puzzled by the fact that the p-values are the same whether I use

warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)

(A):

lsmeans(warp.lm, list(pairwise ~ wool|tension, pairwise ~ tension|wool))

or (B):

lsmeans(warp.lm, list(pairwise ~ wool|tension))

Should the p-values not be higher in the case of (A) than (B)? (A) makes 9 post-hoc comparisons, (B) only 3.

Can I be sure that (A) is adjusted for all 9 comparisons?

Thank you!

Answer 1

Provision in upcoming version of lsmeans

The next update of lsmeans (2.20 or later) will include an rbind method for ref.grid and lsmobj objects. It makes it easy to combine two or ore objects into one family, and defaults to the "mvt" adjustment method. Here is the present example:

> w.t <- pairs(lsmeans(warp.lm, ~ wool | tension))
> t.w <- pairs(lsmeans(warp.lm, ~ tension | wool))

> rbind(w.t, t.w)
 tension wool contrast   estimate       SE df t.ratio p.value
 L       .    A - B    16.3333333 5.157299 48   3.167  0.0203
 M       .    A - B    -4.7777778 5.157299 48  -0.926  0.9119
 H       .    A - B     5.7777778 5.157299 48   1.120  0.8258
 .       A    L - M    20.5555556 5.157299 48   3.986  0.0018
 .       A    L - H    20.0000000 5.157299 48   3.878  0.0027
 .       A    M - H    -0.5555556 5.157299 48  -0.108  1.0000
 .       B    L - M    -0.5555556 5.157299 48  -0.108  1.0000
 .       B    L - H     9.4444444 5.157299 48   1.831  0.3793
 .       B    M - H    10.0000000 5.157299 48   1.939  0.3191

P value adjustment: mvt method for 9 tests 

If you do specify Tukey instead of mvt, it comes fairly close to the same result:

> test(rbind(w.t, t.w), adjust = "tukey")
 tension wool contrast   estimate       SE df t.ratio p.value
 L       .    A - B    16.3333333 5.157299 48   3.167  0.0196
 M       .    A - B    -4.7777778 5.157299 48  -0.926  0.8683
 H       .    A - B     5.7777778 5.157299 48   1.120  0.7727
 .       A    L - M    20.5555556 5.157299 48   3.986  0.0018
 .       A    L - H    20.0000000 5.157299 48   3.878  0.0026
 .       A    M - H    -0.5555556 5.157299 48  -0.108  0.9999
 .       B    L - M    -0.5555556 5.157299 48  -0.108  0.9999
 .       B    L - H     9.4444444 5.157299 48   1.831  0.3467
 .       B    M - H    10.0000000 5.157299 48   1.939  0.2922

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

The curiously fractional family size results from ${4.772 \choose 2} = 9$.

Answer 2

Both of the lsmeans statements you show generate lists of lsmobjs, and each element of those lists is handled separately. If you want to incorporate an overall adjustment for two or more lists combined, it is technical and it takes a bit of work.

First, save the list:

lsmlist = lsmeans(warp.lm, list(pairwise ~ wool|tension, 
                  pairwise ~ tension|wool))

This creates a list of 4 lsmobjs (originally two lists of two)

> names(lsmlist)
[1] "lsmeans of wool | tension"                          
[2] "pairwise differences of contrast, tension | tension"
[3] "lsmeans of tension | wool"                          
[4] "pairwise differences of contrast, wool | wool"

The user wants to combine the three comparisons in lsmlist[[2]] with the six in lsmlist[[4]] and have an overall multiplicity adjustment for those 9 comparisons.

To start, create a new lsmobj from one of the results, and fix it up.

mydiffs = lsmlist[[4]]

First, bind together the linear functions for the two sets of comparisons:

mydiffs@linfct = rbind(lsmlist[[2]]@linfct, lsmlist[[4]]@linfct)

We also need to define the grid slot, which defines the factors associated with each linear function. To make it simple, I just define a factor named contrast with 9 levels for the 9 contrasts (something fancier could be done here).

mydiffs@grid = data.frame(contrast = 1:9)

Finally, fix up the auxiliary info that does the bookkeeping. We now use our new contrast factor as the only variable, with no "by" variables. For the combined family of 9 contrasts, the Tukey adjustment makes no sense so we use the multivariate $t$ ("mvt") method:

mydiffs = update(mydiffs, pri.vars = "contrast", by.vars = NULL,
                 adjust="mvt")

Now, we can look at the resulting summary:

> mydiffs

 contrast   estimate       SE df t.ratio p.value
        1 16.3333333 5.157299 48   3.167  0.0205
        2 -4.7777778 5.157299 48  -0.926  0.9119
        3  5.7777778 5.157299 48   1.120  0.8258
        4 20.5555556 5.157299 48   3.986  0.0019
        5 20.0000000 5.157299 48   3.878  0.0026
        6 -0.5555556 5.157299 48  -0.108  1.0000
        7 -0.5555556 5.157299 48  -0.108  1.0000
        8  9.4444444 5.157299 48   1.831  0.3796
        9 10.0000000 5.157299 48   1.939  0.3190

P value adjustment: mvt method for 9 tests

I could consider adding a feature for combining pieces of lsm.list objects, but it is not straightforward -- primarily because of complications in obtaining meaningful labels for the results (the grid part). It is also a problem that different users expect different defaults.