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

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.

• Please try to this option: "adjust" options are "tukey", "scheffe", "sidak", "bonferroni", "dunnettx", "mvt", and "none". (rcompanion.org/handbook/G_06.html) – J-H Yoon Dec 19 '17 at 2:27
• To me, it's unclear why you think the output from summary(TL_lsmeans1of2\$contrasts) isn't what you want. But, to look at certain contrasts, one solution is to build a matrix which just has the contrasts you are interested in, instead of all pairwise comparisons. Other option is extract the p-values for the comparisons you want and make an adjustment to just these. Pvalue = c(0.0006, 0.0730, 0.0396, 0.0149, 0.1248, 0.4224); Pcorrect = p.adjust(Pvalue, method="fdr"); Pcorrect – Sal Mangiafico Feb 14 '18 at 13:28
• I guess I should mention, too, that you can use the method="none" option in the lsmeans call to get the unadjusted p-values, however many comparisons are being made. – Sal Mangiafico Feb 14 '18 at 15:24

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


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")


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.