How to run 3-way ANOVA on data with no equality of variance in R?

Question

The plan was to run a simple THREE-WAY-ANOVA on my data to analyze the effects of 3 factors on 1 dependent variable.

I have 8 treatments in total, which results from 2 different manipulations of each factor. The Shapiro-Wilk test showed, that I have a normal distribution, which should be fine. The Problem comes up with the second requirement: The equality of variance. The levene-test sadly showed significance.

As you can see below, the biggest variance (group 5) is 2.9 times bigger than the lowest variance (group 1).

Is that ok enough, so that I can just use the "simple" ANOVA (Anova(aov(...), type=3), or is that not allowed? If not, what else could I do? I already read, that I could try: (Anova(aov(...), type=3), white.adjust=True). What´s about that?

Thank you so much in advance!

 Descriptive statistics by group 

group: 1
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 50 5.58 0.81    5.6    5.58 0.89   4   7     3 -0.03    -0.73 0.11
------------------------------------------------------------------- 
group: 2
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 43 5.59 0.84    5.6    5.62 0.89 3.8   7   3.2 -0.12    -0.66 0.13
------------------------------------------------------------------- 
group: 3
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 48 5.04 0.98      5    5.07 1.19 1.8   7   5.2 -0.53     0.73 0.14
------------------------------------------------------------------- 
group: 4
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 41 4.52 1.03    4.6    4.56 1.19 2.4 6.4     4 -0.25    -0.94 0.16
------------------------------------------------------------------- 
group: 5
   vars  n mean   sd median trimmed  mad min max range skew kurtosis  se
X1    1 47 4.29 1.38    4.4    4.28 1.48 1.6 6.8   5.2 0.03    -0.95 0.2
------------------------------------------------------------------- 
group: 6
   vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
X1    1 40 3.53 1.18    3.6    3.47 1.33 1.6 6.4   4.8 0.39    -0.48 0.19
------------------------------------------------------------------- 
group: 7
   vars  n mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 40 4.25 1.13    4.1    4.26 1.19 1.6 6.2   4.6 -0.05    -0.89 0.18
------------------------------------------------------------------- 
group: 8
   vars  n mean  sd median trimmed  mad min max range  skew kurtosis   se
X1    1 41 2.89 0.9    2.8     2.9 0.89   1 4.6   3.6 -0.09    -0.61 0.14

Answer 1

Group sizes are quite large, distributions not too far away from ideal normality and the standard deviations are comparable (Levène's test does not tell you if the deviation from equal variance is problematic). Thus I would expect classic F-tests to provide approximately correct p values. As a sensitivity analysis, adding p values from White-corrected standard errors would be nice. If they lead to different decisions, this will be worth mentioning.

This assumes independent observations, so no repeated measures, hidden clusters etc.

Answer 2

There are some models around that allow for unequal variance.

Below is an example with 'gls' from the 'nlme' package

some random example data

set.seed(1)

# effect sizes    
e0 = 5
e1 = 0.05
e2 = 0.15
e3 = 0.12

# variable  variances
d0 = 0.1
d1 = 0.0
d2 = 0.1
d3 = 0.5


# experiments
com <- matrix(as.numeric(sapply(c(1,2,4),function(x) gl(2,x,8,labels=c(0,1)))),8)
n = c(50,50,50,50,50,50,50,50)

# properties of distributions
m = e0+e1*com[,1]+e2*com[,2]+e3*com[,3]
d = d0+d1*com[,1]+d2*com[,2]+d3*com[,3]

# data matrix generation
v <- c(0,0,0,0)
for (i in 1:8) {
  v <- rbind(v,cbind(matrix(rep(com[i,],n[i]),n[i],byrow=1),rnorm(n[i],m[i],d[i])))
}
v <- v[-1,]
rownames(v) <- 1:sum(n)
colnames(v) <- c("e1","e2","e3","y")

#modelling allowing variable variance
library(nlme)
model1 <- gls(y~e1+e2+e3,data=as.data.frame(v),weights=varIdent(form=~1|e1*e2*e3))
#modelling with fixed variance
model2 <- lm(y~e1+e2+e3,data=as.data.frame(v))

summary(model1)
summary(model2)

confint(model1)
confint(model2)

comparison of output

Coefficients:
  Value  Std.Error  t-value p-value
(Intercept) 5.006409 0.01120563 446.7761  0.0000
e1          0.059496 0.01607292   3.7016  0.0002
e2          0.130702 0.02047058   6.3849  0.0000
e3          0.145693 0.04703072   3.0978  0.0021


Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.99673    0.04775 104.646  < 2e-16 ***
e1           0.05993    0.04775   1.255  0.21016    
e2           0.14993    0.04775   3.140  0.00182 ** 
e3           0.14530    0.04775   3.043  0.00250 ** 

Notice especially the change of the p-value for the e1 factor (similar behavior if you look at anova) which is higher in the second case. This relates to the lower associated variance / standard error in this class (and also the other classes).

See also the changes in ANOVA results

> Anova(model1,type=3)
Analysis of Deviance Table (Type III tests)

Response: y
Df      Chisq Pr(>Chisq)    
(Intercept)  1 156107.951  < 2.2e-16 ***
  e1           1     20.728  5.293e-06 ***
  e2           1     35.887  2.091e-09 ***
  e3           1     10.038   0.001533 ** 
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(model2,type=3)
Anova Table (Type III tests)

Response: y
Sum Sq  Df    F value    Pr(>F)    
(Intercept) 2483.69   1 12177.1800 < 2.2e-16 ***
  e1             0.40   1     1.9451  0.163893    
  e2             2.18   1    10.6999  0.001165 ** 
  e3             2.14   1    10.4823  0.001307 ** 
  Residuals     80.77 396                         
---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1