6

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?

mariachi
  • 213
  • 1
  • 3
  • 7

1 Answers1

7

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)

enter image description here

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.

Stefan
  • 4,977
  • 1
  • 18
  • 38
  • 1
    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. – D_Williams Jan 11 '17 at 03: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: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/full I will add this paper to my answer. It's highly cited and hence has some weight to it. – Stefan Jan 11 '17 at 03:46
  • @D_Williams I am not familiar with Bayesian models but this is on my "to-do" list. Thanks for the comment! – Stefan Jan 11 '17 at 03:47
  • 2
    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. – D_Williams Jan 11 '17 at 03: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. – mariachi Jan 11 '17 at 13:47
  • @Stefan I've been unable to run the leveneTest function. R keeps saying not all arguments are the same length. My variable of interest (treatment) has two levels and both have the same number of rows. I build a model with just value ~ treatment + (1|ID) and used the same command as you before ( leveneTest(residuals(fit), data$treatment) ) – mariachi Jan 12 '17 at 07:02
  • @mariachi How does you exact code look like? Let's assume it's: `mod – Stefan Jan 12 '17 at 15:34