Update on Post-hoc tests following pGLS

[MaximePolicarpo]

Dear Liam,

I recently performed a pGLS, to assess if the copy number in one gene family (taste receptors genes) was associated to the habitat of birds.
I used the caper package and this command :

pgls_habitat <-
  pgls(T2R_number ~ Habitat, 
       data = caper_T2R_eco, 
       lambda = "ML",
       bounds=list(lambda=c(0.5,1)))

I found a significant p-value, and , because the habitat has multiple categories, I now want to make a post-hoc test to see which categories are different. This is when I came across this page of your blog : http://blog.phytools.org/2022/12/post-hoc-tests-for-generalized.html

First, there are apparently some changes in either the nlme or multcomp package, which don't allow to reproduce your example anymore. I came across a solution (that seem to yield similar results to yours). Intead of running glht(primate.ancova,linfct=mcp(Activity_pattern="Tukey")), which returns an error, I run

emmeans_result <- emmeans(primate.ancova, ~ Activity_pattern)
pairwise_comparisons <- contrast(emmeans_result, method = "pairwise", adjust = "tukey")

Now, lets go back to my T2R problem. After finding that the pGLS result was significant, here is what I did as a post-hoc test (same as your blog page):

corBM<-corBrownian(phy=my_bird_tree,form=~species)

my.ancova<-
  gls(T2R_number~
      Habitat,
      data=T2R_eco,
      correlation=corBM)
emmeans_result <- emmeans(my.ancova, ~ AVONET.Habitat)
pairwise_comparisons <- contrast(emmeans_result, method = "pairwise", adjust = "tukey")

I now have few significant pairwise comparisons. I had two main questions :

• Do you still think that such a post-hoc test is valid after your investigations ?
• If the pGLS is also significant, with a lower AIC, using two categorical predictors (Habitat + lifestyle), is it valid to run the same post-hoc method, but adding the same two categorical predictors in the gls, and doing the contrast for each predictor independently ? Like this:

my.ancova<-
  gls(T2R_number~
      Habitat + Lifestyle,
      data=T2R_eco,
      correlation=corBM)

emmeans_result_1 <- emmeans(my.ancova, ~ AVONET.Habitat)
pairwise_comparisons_1 <- contrast(emmeans_result_1, method = "pairwise", adjust = "tukey")

emmeans_result_2 <- emmeans(my.ancova, ~ Lifestyle)
pairwise_comparisons_2 <- contrast(emmeans_result_2, method = "pairwise", adjust = "tukey")

Thank you in advance for any help !

All the best,

Maxime Policarpo