ANOVA
Oscar Ramirez Alan
10 de octubre de 2015
Análisis de Varianza (ANDEVA o ANOVA)
Los Análisis de Varianza (ANDEVA o ANOVA) permiten evidenciar la influencia de un determinado factor o grupo de factores (variables nominales) sobre una variable respuesta (variable continua). Para ello, se comparan efectos que las distintas dosis o niveles del factor producen en la respuesta.
El Anova se puede utilizar en las situaciones en las que nos interesa analizar una respuesta cuantitativa, llamada habitualmente variable dependiente, medida bajo ciertas condiciones experimentales identificadas por una o más variables categóricas (por ejemplo tratamiento, sexo), llamadas variables independientes. Cuando hay una sola variable que proporciona condiciones experimentales distintas, el análisis recibe el nombre de Anova de un factor.
Cuando el efecto producido por un nivel de un factor es influenciado por los niveles de otros factores, se puede afirmar la existencia de una iteracción entre dichos factores, es decir, no actúan de forma independiente.
Cuando tomamos en consideración un único factor para realizar el ANOVA se denomina unifactorial. Caso contrario, cuando se contemplan varios factores se le denomina multifactorial o factorial.
Pruebas Paramétricas
Para realizar un ANOVA con prueba F (de Fischer) se debe cumplir algunos supuestos fundamentales tales como:
Aleatoriedad (dependiendo del Modelo)
Homogeneidad de los residuos
Normalidad de los residuos
El ANOVA descompone la variabilidad de la variable de respuesta entre los diferentes factores. Dependiendo del tipo de análisis, puede ser importante determinar: (a) los factores que tienen un efecto significativo sobre la respuesta, o (b) la cantidad de la variabilidad en la variable de respuesta es atribuible a cada factor.
One-Way ANOVA
Un análisis de varianza unidireccional se utiliza cuando los datos se dividen en grupos de acuerdo con solo un factor. Este análisis se presente cuando se desea saber si: (a) existe una diferencia significativa entre los grupos, y (b) si esto es así, que grupos son significativamente diferentes de las que otros. Algunas pruebas estadísticas se pueden utilizar para comparar las medias de los grupos, las medianas de grupo y desviaciones estándar de grupo. Cuando se comparan las medias, se utilizan múltiples pruebas de gama, el más popular de los cuales es el procedimiento de Tukey (para pruebas paramétricas). Para muestras de igual tamaño, las diferencias significativas por grupo se pueden determinar mediante la interpretación del “Means Plot”. Es muy recomendable que el estudiante se entrene en la identificación de los intervalos y aprenda a reconocer si existe traslape entre ellos o no.
Ejemplo
Se ha realizado un estudio sobre el efecto de tres dietas distintas sobre el peso en individuos de prueba. Los datos obtenidos son:
#Introducimos nuestros datos
dieta1<-c(198,211,240,189,178)
dieta2<-c(214,200,259,194,188)
dieta3<-c(174,176,213,201,158)#Agrupamos las dietas a un solo vector
dietas<- c(dieta1,dieta2,dieta3)
dietas## [1] 198 211 240 189 178 214 200 259 194 188 174 176 213 201 158#Distribuimos los valores, segun los tratamientos
trats <- gl(3,5,label=c("dieta1","dieta2","dieta3")) #3 tratamientos, 5 valores que corresponden a la variabnle del "peso" de la dieta.
length(trats)==length(dietas) #verificamos que tanto los facores como variables dependientes sean del mismo tama??o.## [1] TRUEis.factor(trats) #Confirmamos si tratamientos es un factor## [1] TRUEanova1 <- aov(dietas~trats) #Aplicamos el ANOVA
summary(anova1) #Nos brinda un resumen del ANOVA## Df Sum Sq Mean Sq F value Pr(>F)
## trats 2 1870 934.9 1.496 0.263
## Residuals 12 7500 625.0Forma gráfica utilizando “ggplots”
#para este grafica, se requiere los datos en un data.frame. Además, me permite renombrar los ejes.
di<-data.frame(dietas, trats)
library(ggplot2)
p <- ggplot(di, aes(factor(trats), dietas))
p +scale_y_continuous(name = "Peso (libras)") + scale_x_discrete(name="Tratamientos")+ geom_boxplot()
Tratamientos vs dietas.
Supuestos en el ANOVA
par(mfrow=c(2,2)) #Particiona mi ventana grafica
plot(anova1) #Contrasta los supuestos del ANOVA, de forma grafica
Normalidad
#Deben ser evaluados a traves de los residuos
reg<-lm(dietas~trats) #permite obtener el modelo de regresión lineal
residuos<-residuals(reg) #Estima los residuos
#Normalidad
shapiro.test(residuos) ##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.90598, p-value = 0.1175Homogeneidad de las varianzas
## Bartlett Test de Homogeneidad de Varianzas
bartlett.test(dietas~trats)##
## Bartlett test of homogeneity of variances
##
## data: dietas by trats
## Bartlett's K-squared = 0.24597, df = 2, p-value = 0.8843# Figner-Killeen Test of Homogeneity of Variances, para datos no paramétricos
fligner.test(dietas~trats)##
## Fligner-Killeen test of homogeneity of variances
##
## data: dietas by trats
## Fligner-Killeen:med chi-squared = 0.0027684, df = 2, p-value =
## 0.9986## Levene Test de Homogeneidad de Varianzas
library(car)## Loading required package: carDataleveneTest(dietas~trats)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.0128 0.9873
## 12#cochran.test
library (outliers)
cochran.test(dietas~trats)##
## Cochran test for outlying variance
##
## data: dietas ~ trats
## C = 0.4336, df = 5, k = 3, p-value = 0.8443
## alternative hypothesis: Group dieta2 has outlying variance
## sample estimates:
## dieta1 dieta2 dieta3
## 569.7 813.0 492.3Test Post-Hoc
TukeyHSD: Tukey Honestly Significant Difference test
TukeyHSD(anova1) #Calcula la prueba de Tukey## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = dietas ~ trats)
##
## $trats
## diff lwr upr p adj
## dieta2-dieta1 7.8 -34.38263 49.98263 0.8758009
## dieta3-dieta1 -18.8 -60.98263 23.38263 0.4815568
## dieta3-dieta2 -26.6 -68.78263 15.58263 0.2514855plot(TukeyHSD(anova1)) #Genera un plot para Tukey
ANOVA NO Parametrico
Cuando no cumple con algún supuesto anterior. Aunque se recomienda utilizar la Prueba de Kruskall-Wallis. cuando las varianzas son homogéneas.
kruskal.test(dietas~trats)##
## Kruskal-Wallis rank sum test
##
## data: dietas by trats
## Kruskal-Wallis chi-squared = 2.34, df = 2, p-value = 0.3104Ejemplo1
Utilice los datos “InsectSprays”. Revise la información acerca de lo que tratan los datos.
data(InsectSprays)
attach(InsectSprays)
str(InsectSprays) #Permite ver las características de las variables## 'data.frame': 72 obs. of 2 variables:
## $ count: num 10 7 20 14 14 12 10 23 17 20 ...
## $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...tapply(count, spray, mean) #Permite ver la media de los factores (Tratamientos)## A B C D E F
## 14.500000 15.333333 2.083333 4.916667 3.500000 16.666667tapply(count, spray, sd) #Permite ver la desviación estándar de los factores.## A B C D E F
## 4.719399 4.271115 1.975225 2.503028 1.732051 6.213378tapply(count, spray, length) #Permite ver la cantidad de datos.## A B C D E F
## 12 12 12 12 12 12boxplot(count ~ spray, data=InsectSprays, col="16", xlab="Tipos de Spray", ylab="Cantidad de Insectos")
#o
library(ggplot2)
ggplot(InsectSprays,aes(spray, count, colour = spray ))+geom_point() + geom_boxplot()
ANOVA de una via
andeva2<-aov(count~spray, data=InsectSprays)
summary(andeva2)## Df Sum Sq Mean Sq F value Pr(>F)
## spray 5 2669 533.8 34.7 <2e-16 ***
## Residuals 66 1015 15.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparaciones planeadas:
Cuando tenemos sospecha de un posible resultado del análisis, es decir que entre grupos esperamos encontrar diferencias.
Ejm: cuando en un experimento hay un grupo control, a otro se le aplica un medicamento y otro grupo se le suministra mayor cantidad de un medicamento experimental. Posiblemente sospechemos que vamos a tener resultados diferentes (esto es lo que se entiende en el diseño de muestreo como un Experimento Planeado. Esto por tal debe de tratarse con la prueba de Fisher, en caso de los residuos cumplir con los supuestos paramétricos.
Contrastes no planeados (post-hoc)
Aquí hay múltiples técnicas, y lo que intenta es reducir posibilidad de errores tipo I (no acepta la hipótesis nula siendo esta verdadera), a costa de aumentar el error tipo II. Dicho de otro modo, en situaciones donde haya realmente diferencias entre grupos, las pruebas post-hoc no lo detectan. Por lo que se reconocen como los métodos más seguros, dado que detectan las verdaderas diferencias significativas.
Post Hoc test (Paramétrico)
TukeyHSD(andeva2)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = count ~ spray, data = InsectSprays)
##
## $spray
## diff lwr upr p adj
## B-A 0.8333333 -3.866075 5.532742 0.9951810
## C-A -12.4166667 -17.116075 -7.717258 0.0000000
## D-A -9.5833333 -14.282742 -4.883925 0.0000014
## E-A -11.0000000 -15.699409 -6.300591 0.0000000
## F-A 2.1666667 -2.532742 6.866075 0.7542147
## C-B -13.2500000 -17.949409 -8.550591 0.0000000
## D-B -10.4166667 -15.116075 -5.717258 0.0000002
## E-B -11.8333333 -16.532742 -7.133925 0.0000000
## F-B 1.3333333 -3.366075 6.032742 0.9603075
## D-C 2.8333333 -1.866075 7.532742 0.4920707
## E-C 1.4166667 -3.282742 6.116075 0.9488669
## F-C 14.5833333 9.883925 19.282742 0.0000000
## E-D -1.4166667 -6.116075 3.282742 0.9488669
## F-D 11.7500000 7.050591 16.449409 0.0000000
## F-E 13.1666667 8.467258 17.866075 0.0000000plot(TukeyHSD(andeva2), cex.axis=.7)
Compruebe los supuestos del andeva2
par(mfrow=c(2,2)) #Particiona mi ventana grafica
plot(andeva2)
shapiro.test(residuals(andeva2))##
## Shapiro-Wilk normality test
##
## data: residuals(andeva2)
## W = 0.96006, p-value = 0.02226#Homogeneidad de varianzas
library(car)
leveneTest(andeva2)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 3.8214 0.004223 **
## 66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dado que no cumplen con los supuestos paramétricos, la opción es aplicar es un test no paramétrico.
kruskal.test(count~spray, data=InsectSprays)##
## Kruskal-Wallis rank sum test
##
## data: count by spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10Kruskall-Wallis es muy limitados en cuanto al número de muestras, y se vuelve muy útil cuando el tamaño de muestra es pequeño.
Es importante entender que los test Post-Hoc se deben aplicar cuando en las pruebas se encuentran diferencias significativas. Sobre todo, para contrastar en que grupo es que están ocurriendo esas diferencias.
En cuanto a los contrastes No paramétricos se pueden ejecutar utilizando “pairwise.wilcox.test”, que sería la prueba alternativa al Tukey.
Post-hoc (No paramétricas)
Permiten realizar comparaciones múltiples para identificar los factores significativamente diferentes y así obtener una conclusión sobre qué tratamiento es el adecuado.
Las pruebas más comúnmente utilizadas son Bonferroni, Tukey y LSD (Todas estas deben asumir varianzas iguales)
#LSD test (Least Significant Difference) (Comparación Planeada)
pairwise.wilcox.test(count,spray, p.adj = "none", exact=F)##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: count and spray
##
## A B C D E
## B 0.5812 - - - -
## C 3.8e-05 3.8e-05 - - -
## D 7.8e-05 6.9e-05 0.0027 - -
## E 3.4e-05 3.4e-05 0.0526 0.1744 -
## F 0.4342 0.9078 3.4e-05 7.0e-05 3.4e-05
##
## P value adjustment method: none#corrección de Bonferroni es bastante estricta cuando el número de grupos es grande (más de 6).#Bonferroni (recomendado para tamaño de muestras iguales, medias entre trat. diferentes)
pairwise.wilcox.test(count,spray, p.adj = "bonf",exact=F)##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: count and spray
##
## A B C D E
## B 1.00000 - - - -
## C 0.00058 0.00058 - - -
## D 0.00117 0.00104 0.03977 - -
## E 0.00051 0.00051 0.78860 1.00000 -
## F 1.00000 1.00000 0.00052 0.00105 0.00052
##
## P value adjustment method: bonferroni#corrección de Bonferroni es bastante estricta cuando el número de grupos es grande (más de 6).Algunos test que adecuados cuando las varianzas son desiguales tenemos a:
T2 de Tamhane
T3 de Dunnett
Games-Howell
C de Dunnett
Ejemplo 2
Suponga que se mide en tres estratos diferentes una variable aleatoria, donde se aplica un tratamiento diferente en cada estrato. Se tiene un control que no se aplica ningun efecto. Se quiere medir si existe un crecimiento diferente es los tratamientos asociados al estrato (a,b,y c). Los datos son:
ta <- c(40, 48, 54, 42, 55, 38, 42, 49, 41, 50)
tb <- c(22, 24, 23, 29, 29, 27, 27, 29, 33, 30)
tc <- c(24, 23, 24, 36, 27, 22, 31, 33, 34, 29)
tctl <- c(39, 34, 39, 40, 38, 36, 36, 40, 39, 37)
crecimiento <- c(ta,tb,tc,tctl) #permite agrupar el crecimiento en un solo vector
crecimiento## [1] 40 48 54 42 55 38 42 49 41 50 22 24 23 29 29 27 27 29 33 30 24 23 24
## [24] 36 27 22 31 33 34 29 39 34 39 40 38 36 36 40 39 37trats <- gl(4,10,label=c("ta","tb","tc", "tctl")) #asocia los valores de crecimiento al los tratamientos
length(trats)==length(crecimiento) #confirmamos que tanto las variables independientes como los factores tengan la misma longitud, en caso no de tener la misma longitud, tendremos que utilizar NA, para el valor faltante.## [1] TRUEis.factor(trats) #Confirmamos que el vector es factor## [1] TRUE#Ejecutamos el ANOVA
anova3 <- aov(crecimiento~trats)
#Resumen
summary(anova3)## Df Sum Sq Mean Sq F value Pr(>F)
## trats 3 2307.1 769.0 39.51 1.76e-11 ***
## Residuals 36 700.7 19.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(crecimiento~trats)
#Pueba Post Hoc test
TukeyHSD(anova3)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = crecimiento ~ trats)
##
## $trats
## diff lwr upr p adj
## tb-ta -18.6 -23.91377 -13.28623 0.0000000
## tc-ta -17.6 -22.91377 -12.28623 0.0000000
## tctl-ta -8.1 -13.41377 -2.78623 0.0012165
## tc-tb 1.0 -4.31377 6.31377 0.9569349
## tctl-tb 10.5 5.18623 15.81377 0.0000322
## tctl-tc 9.5 4.18623 14.81377 0.0001498#Grafica de Tukey que mide los efectos
plot(TukeyHSD(anova3))
par(mfrow=c(2,2))
#Confirmamos los supuestos
plot(anova3)
Como redactar sus resultados
Se encontró que el crecimiento medio de los árboles, fue mayor en el estrato A (45.9±6.06, n=10), seguido de estrato control (37.8±1.98, n=10), estrato c(28.3±5.03, n=10) y estrato b (27.3±3.43, n=10)de manera significativa (F=39.51; gl=3,36; p<0.05).Así mismo se realizó una comparación no planeada Tukey HSD, y se encontraron diferencias significativas entre tb-ta (p<0.05), tc-ta (p<0.05), tctl-ta (p<0.05), tctl-tb (p<0.05), tctl-tc (p<0.05),pero no entre tc-tb (p>0.05).
plot(TukeyHSD(anova3))
data.frame(crecimiento,trats)->plantas
library(ggplot2)
p <- ggplot(plantas, aes(factor(trats), crecimiento, colour = trats))
p +scale_y_continuous(name = "Crecimiento") + scale_x_discrete(name="Tratamientos")+ geom_boxplot()
Fig. 9. Crecimiento de las plantas en diferentes tipos de suelos y su comparación de sus medias según ajuste no planeado de TukeyHSD.
ANOVA Factorial
Utlice los datos “ToothGrowth”.
data(ToothGrowth)
str(ToothGrowth) #Tenemos dos variables numericas y un factor## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...head(ToothGrowth)## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5andeva3<-aov(len ~ supp + dose, data=ToothGrowth)
summary(andeva3)## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 11.45 0.0013 **
## dose 1 2224.3 2224.3 123.99 6.31e-16 ***
## Residuals 57 1022.6 17.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Observamos dos efectos significantes "supp" y "dose"
#Revisamos supuestos
par(mfrow=c(2,2))
plot(andeva3)
Si estamos interesados en saber si hay una interaccion entre las variables, hacemos los siguiente.
andeva4<-aov(len ~ supp * dose, data=ToothGrowth) #Observar que se utiliza el "*", para encontrar la interaccion entre "supp"" y "dose".
summary(andeva4)## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 12.317 0.000894 ***
## dose 1 2224.3 2224.3 133.415 < 2e-16 ***
## supp:dose 1 88.9 88.9 5.333 0.024631 *
## Residuals 56 933.6 16.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Volvemos a observar los dos efectos significativos, así como la iteracción.
interaction.plot(ToothGrowth$dose,ToothGrowth$supp, ToothGrowth$len, main="Interaction Plot")
par(mfrow=c(1,2))
boxplot(len ~ supp,data=ToothGrowth)
boxplot(len ~ as.factor(dose),data=ToothGrowth)
#Coplot: Analice la información que se brinda y discútalo con su profesor
coplot(len ~ as.factor(dose) | supp, data=ToothGrowth, panel=panel.smooth)
#Valores numéricos del coplot pueden ser calculados de la siguiente forma
#Permite observar la variación por cada grupo independiente
with(ToothGrowth, tapply(len, list(supp,dose), mean))## 0.5 1 2
## OJ 13.23 22.70 26.06
## VC 7.98 16.77 26.14with(ToothGrowth, tapply(len, supp, mean))## OJ VC
## 20.66333 16.96333with(ToothGrowth, tapply(len, as.factor(dose), mean))## 0.5 1 2
## 10.605 19.735 26.100library(Rcmdr)## Loading required package: splines## Loading required package: RcmdrMisc## Loading required package: sandwich## Loading required package: effects## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## The Commander GUI is launched only in interactive sessions##
## Attaching package: 'Rcmdr'## The following object is masked from 'package:car':
##
## Confintwith(ToothGrowth, plotMeans(len, supp, error.bars="conf.int", level=0.95))
with(ToothGrowth, plotMeans(len, as.factor(dose), error.bars="conf.int", level=0.95))
#Revisamos supuestos
par(mfrow=c(2,2))
plot(andeva4)
Ejemplo 3
suponga que se mide los pesos de una variable aleatoria “x”, en cuatro sitios diferentes, y en dos temporadas (lluviosa y seca). Los datos son:
set.seed(12345)
pesos<-rnorm(40, 1.8,0.70)
sitios<-factor(rep(1:4,10))
sitios## [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3
## [36] 4 1 2 3 4
## Levels: 1 2 3 4#confirmamos la cantidad de datos para sitio
length(sitios)## [1] 40temporada<-factor(rep(c("lluvioso", "seco"), each=20))
temporada## [1] lluvioso lluvioso lluvioso lluvioso lluvioso lluvioso lluvioso
## [8] lluvioso lluvioso lluvioso lluvioso lluvioso lluvioso lluvioso
## [15] lluvioso lluvioso lluvioso lluvioso lluvioso lluvioso seco
## [22] seco seco seco seco seco seco seco
## [29] seco seco seco seco seco seco seco
## [36] seco seco seco seco seco
## Levels: lluvioso seco#confirmamos la cantidad de datos para temporada
length(temporada)## [1] 40#Comparamos que ambas tengan la misma cantidad de valores
length(pesos)==length(sitios)## [1] TRUEexperimento<-data.frame(sitios, temporada, pesos)
andeva5<-aov(pesos~sitios+temporada)
summary(andeva5)## Df Sum Sq Mean Sq F value Pr(>F)
## sitios 3 0.202 0.0674 0.119 0.949
## temporada 1 0.525 0.5250 0.924 0.343
## Residuals 35 19.880 0.5680#Buscamos si hay interaccion
andeva6<-aov(pesos~sitios*temporada, data=experimento)
summary(andeva6)## Df Sum Sq Mean Sq F value Pr(>F)
## sitios 3 0.202 0.0674 0.113 0.952
## temporada 1 0.525 0.5250 0.879 0.355
## sitios:temporada 3 0.769 0.2564 0.429 0.733
## Residuals 32 19.111 0.5972with(experimento, plotMeans(pesos, as.factor(sitios), error.bars="conf.int", level=0.95))
#Revisamos supuestos
par(mfrow=c(2,2))
plot(andeva6)