Test homogeneity in lmer models

A reviewer is asking me to test for homoscedasticity with "appropriate tests", as opposed to visual inspection of residual plots. I haven't found any such test (I guess he wants a p-value). Does such a test exist and if not, is there something I can cite to justify not presenting a p-value?

Yes you could use for example Levene's test using the leveneTest() function from the car package. Here's an example with the Machines dataset from the nlme package:

library(lme4)
library(nlme); data(Machines)
library(car)

mod <- lmer(score ~ Machine + (1|Worker), data=Machines)

> leveneTest(residuals(mod) ~ Machines$Machine) Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F) group 2 0.0811 0.9222 51  Since the result is not significant, the assumption of equal variances (homoscedasticity) is met. Also check ?leveneTest for more options. Let's compare this with a boxplot of the residuals: boxplot(residuals(mod) ~ Machines$Machine)


Since the reviewer seems to want a "formal test", it will probably be difficult to convince him accepting your visual inspection, despite, in my opinion, this would be the way to go. Maybe someone else has an actual reference why checking those assumptions visually is superior compared to "formal tests".

Edit to address comment by @D_Williams below:

A good and strongly cited paper by Zuur et al. 2010 may help to convince your reviewer regarding visual inspection of residuals to test for homogeneity of variances. Also here's is the link to the book Mixed Effects Models and Extensions in Ecology with R.

• Yes, a great citation for just this is Zuur et al. generalized linear mixed models in ecology. It is a book. You might add that variances are always unequal to some degree. The question is the point at which this become problematic. I would suggest to use brms to fit a Bayesian model in which you can fit a so called unequal variance model. Then, compare this via LOO and posterior predictive checks to a model that assumes equal variances. I generally always fit the unequal variance model, because it makes the most sense. Jan 11, 2017 at 3:35
• @D_Williams Thanks for your comment! Yes, now that you mentioned Zuur et al, I also remembered they had a paper on statistical data analysis: onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/… I will add this paper to my answer. It's highly cited and hence has some weight to it. Jan 11, 2017 at 3:46
• @D_Williams I am not familiar with Bayesian models but this is on my "to-do" list. Thanks for the comment! Jan 11, 2017 at 3:47
• No problem. Brms uses same arguments as lme4. for unequal variances, suppose we have an interaction and we want to model the variances for each level of the interaction: y ~ x1 * x2, sigma ~ x1 * x2, where sigma...is for modeling the variances, or standars deviations. Very straightforward and highly recommended. Jan 11, 2017 at 3:53
• Thank you both. I hope the citation will satisfy her/him. Most likely not, but fingers crossed. If you have any other relevant citations, please add them nonetheless. Jan 11, 2017 at 13:47