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

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


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.

• In addition, power is better with weights for the unequal group sizes (and also apply weights to correct for the unequal variance). If left uncorrected then a larger group size (if not balanced) may actually increase the type I or type II error rate. ie take $10$ samples from group A with $\mathcal{N}(1,0.5)$ and $k$ samples from group B $\mathcal{N}(0,2)$ then $(\mu_A=\mu_B)$ is less likely rejected for larger $k$. Similarly take $10$ samples from group A with $\mathcal{N}(0,2)$ and $k$ samples from group B $\mathcal{N}(0,0.5)$ then $(\mu_A=\mu_B)$ is more likely rejected for larger $k$. – Martijn Weterings Jul 6 '17 at 12:52
• Yes, I have independent observations :).... Thank you! – Timo Jul 6 '17 at 13:04

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).

> 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