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):
enter image description here

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?

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

I assume I can use Tukey's test then?

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:

Similar Posts:

Rate this post

Leave a Comment