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!
Best Answer
Both of the
lsmeans
statements you show generate lists oflsmobj
s, 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:
This creates a list of 4
lsmobj
s (originally two lists of two)The user wants to combine the three comparisons in
lsmlist[[2]]
with the six inlsmlist[[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.First, bind together the linear functions for the two sets of comparisons:
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 namedcontrast
with 9 levels for the 9 contrasts (something fancier could be done here).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:Now, we can look at the resulting summary:
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 (thegrid
part). It is also a problem that different users expect different defaults.