Solved – Using pairwise comparison on gls object with heterogeneity of variances

I'm comparing total yield from fields under 5 different treatments. As you can see, the variance differed between the treatments (diagnostics from lm() fit):

So I fitted the model with nlme::gls() and a variance structure depending on the treatment:

vs <- varIdent(form= ~1|Treatment)  ls1 <- gls(Yield ~ Treatment, weights=vs, data=prod, method="ML") 

When doing an anova on this model and one without 'Treatment', the effect of treatment on yield proved to be significant (p<0.05). So I wanted to see which treatments differed from each other and found the package emmeans (used to be lsmeans) which can work with gls objects.

When I run emmeans(ls1) I get the desired output and at the bottom it says "P value adjustment: tukey method for comparing a family of 5 estimates"

I learned that one of the assumptions of Tukey's test is homoscedasticity, which is not the case here (that's why I used gls() to begin with).

So I cannot use the results of emmeans()?
What other options do I have?

This analyses is for a paper, can I just describe the results I've plotted without any post hoc tests and thus without reporting p-values?
I've read several times that post hoc tests and reporting p-values are overrated, but then how do I report? "Figure x clearly show that treatment A differs from treatment B"? Does that suffice?

EDIT:
I checked the assumptions again with the gls() modelfit and these are the standardized residuals vs the fitted values now:

I assume I can use Tukey's test then?

model.matrix.gls <- function(object, ...) {     model.matrix(terms(object), data = getData(object), ...) } model.frame.gls <- function(object, ...) {     model.frame(formula(object), data = getData(object), ...) } terms.gls <- function(object, ...) {     terms(model.frame(object), ...) }  phAnova<-glht(ls1, linfct=mcp(Treatment="Tukey")) summary(phAnova) 

Similar Posts:

• Solved – How to compare two levels of one factor
• Solved – Interpreting R lm() output with factor as variable
• Solved – Easy post-hoc tests when meta-analyzing with the `metafor` package in r
• Solved – R: Post-hoc test on lmer. emmeans and multcomp packages
• Solved – multcomp() vs emmeans() for multiple comparisons