Using pairwise comparison on gls object with heterogeneity of variances?

Question

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?

Answer 1

You may use the glht() function in the multcomp package if you run the following code before:

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)

Look at this link, down the page: http://rstudio-pubs-static.s3.amazonaws.com/13472_0daab9a778f24d3dbf38d808952455ce.html.