ANOVA y Kruskal-Wallis
junio 24, 2022
1 Autor
Autor: Álvaro Alonso Fernández
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
2 Nuestros datos
Queremos comparar si hay diferencias en el valor de una variable dependiente entre tres o más grupos (=tratamientos). Por ejemplo, una aplicación de un test ANOVA sería la comparación estadística entre la altura de los habitantes de tres ciudades (Variable dependiente=Altura, grupos=la ciudad que tiene tres niveles Ciudad A, B y C). Otro ejemplo sería evaluar el efecto de diferentes fertilizantes sobre el crecimiento de una especie de planta.
Cargamos nuestra matriz de datos, que tendrá el siguiente aspecto:
setwd(dir = "F:/R/MARKDOWN/anovaykwoneway")
datos<-read.table("ANOVA.csv", sep = ";", header = TRUE, dec = ".")
datos## Treatment Var1 Var2
## 1 1 2.2 0.20
## 2 1 2.3 0.30
## 3 1 2.4 0.40
## 4 1 4.5 0.50
## 5 1 5.4 0.20
## 6 1 5.4 0.40
## 7 1 5.4 0.40
## 8 1 4.0 1.00
## 9 1 4.0 0.90
## 10 1 4.4 0.80
## 11 2 0.0 0.90
## 12 2 0.0 0.80
## 13 2 4.0 0.90
## 14 2 5.0 1.20
## 15 2 6.0 1.30
## 16 2 8.0 1.30
## 17 2 9.0 1.20
## 18 2 8.0 1.40
## 19 2 0.0 2.00
## 20 2 9.0 2.00
## 21 3 9.0 3.20
## 22 3 8.0 3.40
## 23 3 7.8 4.50
## 24 3 8.9 5.60
## 25 3 12.0 6.50
## 26 3 2.0 5.40
## 27 3 12.0 4.50
## 28 3 11.0 6.40
## 29 3 14.0 6.70
## 30 3 3.0 6.90
## 31 4 45.0 9.00
## 32 4 54.0 9.80
## 33 4 46.0 9.80
## 34 4 56.0 9.80
## 35 4 76.0 9.90
## 36 4 54.0 9.00
## 37 4 54.0 7.00
## 38 4 34.0 8.00
## 39 4 45.0 9.00
## 40 4 44.0 6.00
## 41 5 22.0 6.80
## 42 5 33.0 9.00
## 43 5 33.0 9.00
## 44 5 33.0 9.90
## 45 5 33.0 9.89
## 46 5 33.0 9.90
## 47 5 32.0 6.90
## 48 5 32.0 9.90
## 49 5 43.0 9.99
## 50 5 55.0 9.80str(datos)## 'data.frame': 50 obs. of 3 variables:
## $ Treatment: int 1 1 1 1 1 1 1 1 1 1 ...
## $ Var1 : num 2.2 2.3 2.4 4.5 5.4 5.4 5.4 4 4 4.4 ...
## $ Var2 : num 0.2 0.3 0.4 0.5 0.2 0.4 0.4 1 0.9 0.8 ...En este caso tenemos cinco grupos y queremos comparar dos variables dependientes (Var1 y Var2).
3 ¿Qué pinta tienen mis datos?
Hacemos un Boxplot con nuestros datos:
boxplot(Var1 ~ Treatment, data=datos, col="green", cex.axis=0.7,las = 2, ylab="var1", xlab="Treatments", cex.lab=0.75)
boxplot(Var2 ~ Treatment, data=datos, col="red", cex.axis=0.7,las = 2, ylab="var2", xlab="Treatments", cex.lab=0.75)
4 ANOVA
Hacemos un modelo ANOVA para nuestra Variable 1 (se haría lo mismo para la Variable 2):
nuestraanova <- lm(Var1~Treatment, data=datos)
nuestromodelo <- anova(nuestraanova)
nuestromodelo## Analysis of Variance Table
##
## Response: Var1
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 1 11599 11599.3 67.536 1.033e-10 ***
## Residuals 48 8244 171.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p valor es muy bajo, nos indica diferencias en el valor de la variable entre los tratamientos (en al menos un grupo).
5 ¿Son nuestros datos normales y homogéneos?
El test ANOVA requiere normalidad de los residuos y homogeneidad en las varianzas.
Representamos visualmente los residuos:
par(mfrow=c(2,2))
plot(nuestraanova)
La figura Q-Q plot muestra residuos alejados de la l?nea, esto no pinta muy bien…..Deberían aparecer cercanos o encima de la línea en caso de ser normales o aproximarse a la normalidad.
Hagamos un test para comprobar la normalidad:
shapiro.test(residuals(nuestraanova))##
## Shapiro-Wilk normality test
##
## data: residuals(nuestraanova)
## W = 0.91811, p-value = 0.002015Nuestros residuos no son normales (p<0.05).
fligner.test(Var1 ~ Treatment, data = datos)##
## Fligner-Killeen test of homogeneity of variances
##
## data: Var1 by Treatment
## Fligner-Killeen:med chi-squared = 13.058, df = 4, p-value = 0.011Tampoco se cumple el requisito de la homogeneidad de varianzas.
Si el test cumpliese normalidad y homogeneidad, el siguiente paso sería comparar el valor de la variable dependiente entre tratamientos de dos en dos. Para ello utilizaríamos un t test con una corrección de Bonferroni:
pairwise.t.test(x = datos$Var1, g = datos$Treatment,
p.adjust.method = c("bonferroni"))##
## Pairwise comparisons using t tests with pooled SD
##
## data: datos$Var1 and datos$Treatment
##
## 1 2 3 4
## 2 1 - - -
## 3 1 1 - -
## 4 < 2e-16 < 2e-16 < 2e-16 -
## 5 2.6e-12 6.7e-12 4.1e-10 3.9e-05
##
## P value adjustment method: bonferroni6 Alternativa no paramétrica: Kruskal-Wallis
Hacemos un test Kruskal-Wallis:
kruskal.test(Var1 ~ Treatment, data = datos)##
## Kruskal-Wallis rank sum test
##
## data: Var1 by Treatment
## Kruskal-Wallis chi-squared = 39.621, df = 4, p-value = 5.185e-08Hay diferencias en la variable dependiente Var1 entre los tratamientos.
¿Dónde se encuentran esas diferencias? Para responder a esta pregunta comparamos dos a dos los tratamientos:
pairwise.wilcox.test(datos$Var1,datos$Treatment,
p.adjust.method = c("bonferroni"),paired=F,exact=F)##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: datos$Var1 and datos$Treatment
##
## 1 2 3 4
## 2 1.0000 - - -
## 3 0.1380 0.5273 - -
## 4 0.0017 0.0017 0.0018 -
## 5 0.0016 0.0016 0.0016 0.0199
##
## P value adjustment method: bonferroniEn la tabla podemos ver los tratamientos que difieren entre si por parejas, todos los que tengan p<0.05 son distintos entre si.
7 Representación gráfica de resultados
Podemos representar nuestros resultados en una figura de barras indicando con un asterisco las diferencias con el control. Primero calculamos la media y la desviación estándar en función del tratamiento:
mediasvar1<-aggregate(datos$Var1~datos$Treatment, FUN=mean)
sdvar1<-aggregate(datos$Var1~datos$Treatment, FUN=sd)
mediasvar1## datos$Treatment datos$Var1
## 1 1 4.00
## 2 2 4.90
## 3 3 8.77
## 4 4 50.80
## 5 5 34.90sdvar1## datos$Treatment datos$Var1
## 1 1 1.290133
## 2 2 3.754997
## 3 3 3.857475
## 4 4 11.113555
## 5 5 8.633912library(plyr)
medvar1<-ddply(datos,.(Treatment), summarize, mean=mean(Var1))
sdvar1<-ddply(datos,.(Treatment), summarize, sd=sd(Var1))
str(medvar1)## 'data.frame': 5 obs. of 2 variables:
## $ Treatment: int 1 2 3 4 5
## $ mean : num 4 4.9 8.77 50.8 34.9str(sdvar1)## 'data.frame': 5 obs. of 2 variables:
## $ Treatment: int 1 2 3 4 5
## $ sd : num 1.29 3.75 3.86 11.11 8.63Ahora hacemos una figura de barras sin los ejes (se los pondremos más abajo con axis):
BARRAS<-barplot(medvar1$mean, axes=FALSE,axisname=FALSE, ylim=c(0,100),
col=c('white', 'gray','azure4','dimgray','black'),main="Var1",
xlab="Treatments", ylab="Var1")
axis(1,labels=c("control", "T1","T2","T3","T4"), at=BARRAS)
axis(2,at=seq(0,100,by=10))
Y ahora ponemos las barras de error segments:
BARRAS<-barplot(medvar1$mean, axes=FALSE,axisname=FALSE, ylim=c(0,100),
col=c('white', 'gray','azure4','dimgray','black'),main="Var1",
xlab="Treatments", ylab="Var1")
axis(1,labels=c("control", "T1","T2","T3","T4"), at=BARRAS)
axis(2,at=seq(0,100,by=10))
segments(BARRAS-0.1,medvar1$mean-sdvar1$sd,BARRAS+0.1,medvar1$mean-sdvar1$sd,lwd=2)
segments(BARRAS-0.1,medvar1$mean+sdvar1$sd,BARRAS+0.1,medvar1$mean+sdvar1$sd,lwd=2)
segments(BARRAS,medvar1$mean-sdvar1$sd,BARRAS,medvar1$mean+sdvar1$sd,lwd=2)
Ahora ponemos unos asteriscos test para indicar las diferencias, en este caso con el control:
BARRAS<-barplot(medvar1$mean, axes=FALSE,axisname=FALSE, ylim=c(0,100),
col=c('white', 'gray','azure4','dimgray','black'),main="Var1",
xlab="Treatments", ylab="Var1")
axis(1,labels=c("control", "T1","T2","T3","T4"), at=BARRAS)
axis(2,at=seq(0,100,by=10))
segments(BARRAS-0.1,medvar1$mean-sdvar1$sd,BARRAS+0.1,medvar1$mean-sdvar1$sd,lwd=2)
segments(BARRAS-0.1,medvar1$mean+sdvar1$sd,BARRAS+0.1,medvar1$mean+sdvar1$sd,lwd=2)
segments(BARRAS,medvar1$mean-sdvar1$sd,BARRAS,medvar1$mean+sdvar1$sd,lwd=2)
text(4.3,69,labels="*",cex=6)
text(5.5,55,labels="*",cex=6)