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.
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 allpairwise
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
$\endgroup$method="none"
option in the lsmeans call to get the unadjusted p-values, however many comparisons are being made. $\endgroup$