How can I program correction for multiple comparisons using lsmeans in R while not correcting for more than necessary

Question

I posted this question on stack overflow here enter link description here but did not get useful feedback. So, I'll post it again here, in the hopes that this is a better venue for it. I'm wondering if there is a way to correct for multiple comparisons using lsmeans that will allow me to correct for only a certain number of comparisons?

The best model to explain my data is a complete model, and I'm looking at interactions between three factors (year, river and sex) (there are 4 rivers and 4 different years). I could do the following (TL is total length, and is a response variable)

TL_interact = lm(TL ~ YearRiverSex, data = Walleye_alldata_nu) 

followed by the lsmeans

turtle = lsmeans(possum, specs = pairwise~PeriodRiverSex, adjust = "FDR")

But this corrects for 496 comparisons. The actual number of comparisons that I am interested in is 56. These can be achieved by using the slightly altered code

TL_lsmeans1of2 = lsmeans(TL_interact, specs = pairwise~Year|Sex*River, adjust = "FDR")

then using summary I can manipulate this data as a dataframe

TL_lsmeans1of2_df = summary(TL_lsmeans1of2$contrasts)

a sample of the output from this is

Sex = F, River = Chalifour:
 contrast            estimate        SE   df t.ratio p.value
 Y2002/03 - Y2015  74.1746032 18.994645 1497   3.905  0.0006
 Y2002/03 - Y2016  33.4750958 16.963837 1497   1.973  0.0730
 Y2002/03 - Y2017  45.7222222 19.604254 1497   2.332  0.0396
 Y2015 - Y2016    -40.6995074 14.468508 1497  -2.813  0.0149
 Y2015 - Y2017    -28.4523810 17.489789 1497  -1.627  0.1248
 Y2016 - Y2017     12.2471264 15.260011 1497   0.803  0.4224

Sex = M, River = Chalifour:
 contrast            estimate        SE   df t.ratio p.value
 Y2002/03 - Y2015  49.4788034  6.054656 1497   8.172  <.0001
 Y2002/03 - Y2016  36.4539394  5.893992 1497   6.185  <.0001
 Y2002/03 - Y2017  55.7266667  6.352504 1497   8.772  <.0001
 Y2015 - Y2016    -13.0248640  5.645109 1497  -2.307  0.0254
 Y2015 - Y2017      6.2478632  6.122289 1497   1.021  0.3077
 Y2016 - Y2017     19.2727273  5.963447 1497   3.232  0.0019

and then, selecting certain rows of the output for

TL_lsmeans2of2 = lsmeans(TL_interact, specs = pairwise~Year*Sex|River, adjust = "FDR")

I won't talk about this second lsmeans output b/c the same solution to my question will be applied there. However, with the lines that I select from this second lsmeans statement, together with the results of the first lsmeans statement, I have 56 comparisons.

I tried to adjust the p-value using the t.ratio and the degrees of freedom using the pt command using

TL_lsmeans1of2_df$p.val.adjusted = pt(TL_lsmeans1of2_df$t.ratio, 
TL_lsmeans1of2_df$df)

But the results are just wacky

Sex = F, River = Chalifour:
 contrast            estimate        SE   df t.ratio p.value p.val.adjusted
 Y2002/03 - Y2015  74.1746032 18.994645 1497   3.905  0.0006    0.999950797
 Y2002/03 - Y2016  33.4750958 16.963837 1497   1.973  0.0730    0.975678658
 Y2002/03 - Y2017  45.7222222 19.604254 1497   2.332  0.0396    0.990090497
 Y2015 - Y2016    -40.6995074 14.468508 1497  -2.813  0.0149    0.002486354
 Y2015 - Y2017    -28.4523810 17.489789 1497  -1.627  0.1248    0.051995069
 Y2016 - Y2017     12.2471264 15.260011 1497   0.803  0.4224    0.788822739

How could I get a p-value using an appropriate correction for 56 comparisons? And, if you think it is better for me to code my contrasts differently (I saw this post using 0s and 1s [https://rcompanion.org/rcompanion/h_01.html], I would also gladly take any suggestions on how to approach this.

---

I think that I might not have been completely clear. I have a different solution though. I will do standard t-tests for each contrast, and then extract those p-values and FDR correct them.

Answer 1

I think what's most useful to you is that the [] operator is defined for lsmeans() results (also emmeans::emmeans()). So you can do:

TL_lsm = lsmeans(TL_interact, ~ Year*Sex|River)
TL_lsm_pairs = pairs(TL_lsm)
wanted = c(1,2,3,5,8,13,21)   ## whichever ones you want
summary(TL_lsm_pairs[wanted], adjust = "fdr")

There are also other contrast families besides all pairwise ones. For example,

contrast(TL_lsm, "consec")

compares only consecutive levels. There is also "trt.vs.ctrl1" and "trt.vs.ctrlk" for comparing each level with the first or last level. Perhaps one of the available contrast families is already what you want.

Finally, in emmeans, there is a provision to do simple comparisons (comparisons of each factor at each combination of levels of the other two factors). That may not help in this example, but it's often useful for getting simple effects:

library(emmeans)
TL_emm = emmeans(TL_interact, ~Year*Sex*River)
pairs(TL_emm, simple = "each")

Answer 2

You don't want to use individual t-tests. That will ignore the complexity of the model. There's no point in building the model if you are just going to ignore it.

Instead, use lsmeans (now emmeans). Just use the method="none" option, and then correct the p-values by any method you want.

As you mention, another approach is to define the contrasts you are interested in, and then test only those contrasts.