ANOVA
Nataly Cuichán
28 de septiembre de 2018
Idea intuitiva del ANOVA
La técnica de análisis de varianza (ANOVA) también conocida como análisis factorial y desarrollada por Fisher en 1930, constituye la herramienta básica para el estudio del efecto de uno o más factores (cada uno con dos o más niveles) sobre la media de una variable continua.
Es por lo tanto el test estadístico a emplear cuando se desea comparar las medias de dos o más grupos. Esta técnica puede generalizarse también para estudiar los posibles efectos de los factores sobre la varianza de una variable.
La hipótesis nula de la que parten los diferentes tipos de ANOVA es que la media de la variable estudiada es la misma en los diferentes grupos, en contraposición a la hipótesis alternativa de que al menos dos medias difieren de forma significativa.
ANOVA permite comparar múltiples medias, pero lo hace mediante el estudio de las varianzas. 
El funcionamiento básico de un ANOVA consiste en calcular la media de cada uno de los grupos para a continuación comparar la varianza de estas medias (varianza explicada por la variable grupo, intervarianza) frente a la varianza promedio dentro de los grupos (la no explicada por la variable grupo, intravarianza).
Bajo la hipótesis nula de que las observaciones de los distintos grupos proceden todas la misma población (tienen la misma media y varianza), la varianza ponderada entre grupos será la misma que la varianza promedio dentro de los grupos.
Conforme las medias de los grupos estén más alejadas las unas de las otras, la varianza entre medias se incrementará y dejará de ser igual a la varianza promedio dentro de los grupos.
Análisis de la Varianza (ANOVA)
Ya que los procedimientos ANOVA de dos muestras permiten comparar las medias de dos poblaciones o las respuestas medias a dos tratamientos de un experimento. Sin embargo, en ocasiones necesitamos comparar más de 2 grupos. El modelo del Análisis de la Varianza (ANOVA), nos permitirá abordar este tipo de situaciones. Lo vemos con un ejemplo:
Estamos interesados en conocer si hay colores más atractivos para los insectos. Para ello se diseñaron trampas con los siguientes colores: amarillo, azul, blanco y verde. Se cuantificó el número de insectos que quedaban atrapados:
- Azul: 16 11 20 21 14 7
- Verde: 37 32 15 25 39 41
- Blanco: 21 12 14 17 13 17
- Amarillo: 45 59 48 46 38 47
Organización de los datos
Generamos dos variables: insectos es la variable respuesta y colores es la variable factor (establece los grupos de interés):
insectos <- c(16,11,20,21,14,7,37,32,15,25,39,41,21,12,14,17,13,17,45,59,48,46,38,47)
colores <- as.factor(c(rep(c("azul", "verde", "blanco", "amarillo"), each =6)))Exploramos los datos de la muestra:
boxplot(insectos ~ colores, col = c("yellow", "blue", "white","green"), ylab = "Número de insectos atrapados") #amarillo casi no tiene varianza pero si dos datos atipicos
tapply(insectos, colores, mean)## amarillo azul blanco verde
## 47.16667 14.83333 15.66667 31.50000ANOVA y pruebas post-hoc
Esta es la forma de pedir un ANOVA en R:
#asi se inicia un ANOVA (continua, categorica)
#si hacemos (continua, continua) no tienen los mismos resultados
fm = aov( lm(insectos ~ colores) )Pedimos un resumen de la tabla del ANOVA
summary(fm)## Df Sum Sq Mean Sq F value Pr(>F)
## colores 3 4218 1406 30.55 1.15e-07 ***
## Residuals 20 921 46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#ya que es estadistic signific podemos concluir que las medias no son iguales y por lo tanto los insectos diferencian el colorElementos generados en el ANOVA:
names(fm)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"Identifica en la tabla ANOVA los grados de libertad del factor, los grados de libertad residuales, la suma de cuadrados de los grupos, la suma de cuadrados del error, las medias correspondientes de las sumas de cuadrados de los grupos y del error, el valor del estadístico F. Describe cómo obtenemos cada uno de ellos.
¿Cuál es el valor crítico de F bajo la hipótesis nula con un nivel de significación alfa = 0.05? (Este valor nos delimitará la región de aceptación y rechazo)
Bajo la Ho el estadístico de contraste F se distribuye como una F de grados de libertad (I-1), (n-I) donde I es el número de grupos que disponemos y n el tamaño total de la muestral. Así obtenemos el cuantil buscado:
#Asi se saca el valor en tablas
#qf es la prueba de (nivel de signici, I-1,n-1,si tiene cola inferior)
#Xi cuadrado
qf(0.05, 3-1, 18-3, lower.tail = F)## [1] 3.68232Valores del estadístico > 3.68232 estarán incluidos en la región de rechazo. En nuetro caso 30.55 es mucho mayor que el valor crítico obtenido.
¿Qué valor de la tabla ANOVA nos proporciona la varianza muestral común (pooled variance en inglés)? ¿Para qué es útil?
La raíz cuadrada de la media de los cuadrados del error, además de proporcionarnos una estimación de la varianza muestral de todos los datos, se utiliza en la obtención de los intervalos de confianza de las medias en cada uno de los grupos de interés.
Por ejemplo, este sería el intervalo de confianza de la media de los insectos capturados para las trampas amarillas, con un nivel de confianza del 95%:
media <- mean(insectos[colores =="amarillo"])
valor_t <- pt(0.05/2, 18 - 3) #0.05 son los grados de confianza
sp <- sqrt(46) #desviación tÃpica de la varianza muestral común
ee <- valor_t * (sp/ sqrt(6)) #error de estimación
media## [1] 47.16667límite superior del intervalo de confianza de la media de insectos capturados para las trampas amarillas
media + ee ## [1] 48.57826límite inferior del intervalo de confianza de la media de insectos capturados para las trampas amarillas
media - ee ## [1] 45.75507Si hemos detectado diferencias significativas entre las medias de las poblaciones. ¿Será posible saber cuáles son los grupos que generan estas diferencias?
intervals = TukeyHSD(fm)
intervals## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm(insectos ~ colores))
##
## $colores
## diff lwr upr p adj
## azul-amarillo -32.3333333 -43.296330 -21.37034 0.0000004
## blanco-amarillo -31.5000000 -42.462996 -20.53700 0.0000006
## verde-amarillo -15.6666667 -26.629663 -4.70367 0.0036170
## blanco-azul 0.8333333 -10.129663 11.79633 0.9964823
## verde-azul 16.6666667 5.703670 27.62966 0.0020222
## verde-blanco 15.8333333 4.870337 26.79633 0.0032835#Si en los intervalos contienen al cero son estad iguales si no no lo son,
#en este ejemplo blanco azul si loocntiene
plot(intervals)
Explica las diferencias existentes por parejas de trampas según el color. ¿Algunas de estas diferencias son significativas? Si el objetivo es atrapar un mayor número de insectos, ¿con qué tipo de trampas te quedarás?
Validación del modelo ANOVA
A partir de los residuos del modelo comprobaremos si el modelo ANOVA es adecuado. Los supuestos que se deben cumplir son tres: independencia, homocedasticidad y normalidad.
Independencia
#Mientras mas dispersos menos correlacionados
plot(fm$residuals)
Normalidad
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:
summary(fm$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -16.5000 -2.9167 0.1667 0.0000 5.2083 11.8333boxplot(fm$residuals)
hist(fm$residuals)
#no hay datos atipicos y casi estan sobre la linea normal
#ya que no dan outlets nos da un buen indicio de que son normales
qqnorm(fm$residuals)
qqline(fm$residuals)
El test de Shapiro-Wilk indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (normalidad de los residuos)
shapiro.test(fm$residuals)##
## Shapiro-Wilk normality test
##
## data: fm$residuals
## W = 0.97337, p-value = 0.75#Hasta 0.05 son normales
#este si es normalHomocedasticidad
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:
#En este grafico nos fijamos en el ancho no en la media
#el que mas dispersion tiene es el verde
#la varianza de los grupos se obtiene a traves de los residuos pq no esta mezclada con la media
boxplot(fm$residuals~colores, col = c("yellow", "blue", "white","green")) 
desviaciones <- tapply(fm$residuals, colores, sd)
#sd desviacion estandarComparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente)
#como es estricto mayor a 2 podemos decir que son de varianzas diferentes
#no hay evidencia suficiente para rechazar la H0
max(desviaciones) / min(desviaciones) ## [1] 2.980357El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales)
#Afirmamos con un test que las varianzas son diferentes
bartlett.test(fm$residuals ~ colores)##
## Bartlett test of homogeneity of variances
##
## data: fm$residuals by colores
## Bartlett's K-squared = 5.2628, df = 3, p-value = 0.1535Kruskal-Wallis y pruebas post-hoc
¿Qué hipótesis contrasta el test ANOVA?
Ho: las medias son iguales en todas las poblaciones
Ha: hay alguna media distinta
¿Qué hipótesis contrasta la prueba de Kruskal-Wallis?
Ho: la variable respuesta es la misma en todas las poblaciones valoradas
Ha: la variable respuesta es mayor en alguna de las poblaciones
Cuando no se cumplen las hipótesis exigidas por el modelo ANOVA, es posible utilizar la prueba no paramétrica Kruskal-Wallis: ¿hay diferencias significativas entre las poblaciones?
kruskal.test(insectos, colores)##
## Kruskal-Wallis rank sum test
##
## data: insectos and colores
## Kruskal-Wallis chi-squared = 16.975, df = 3, p-value = 0.000715Indica cuál es el estadístico de contraste, los grados de libertad, el p-valor correspondiente y cuál sería el valor crítico que definiría las regiones de aceptación y rechazo con un nivel de significación alfa = 0.05.
Bajo la Ho el estadístico de contraste H del test de Kruskal-Wallis se distribuye como una Chi-cuadrado de grados de libertad (I-1) (donde I es el número de grupos que disponemos). Así obtenemos el cuantil buscado:
qchisq(0.05, 3-1, lower.tail = F)## [1] 5.991465Valores del estadístico > 5.991465 estarán incluidos en la región de rechazo. En nuetro caso 16.9755 es mucho mayor que el valor crítico obtenido.
Si transformáramos los datos de la variable respuesta, utilizando logaritmos y después aplicáramos el test de KrusKal-Wallis al logaritmo del número de insectos atrapados, ¿variarían los resultados del test estadístico?
kruskal.test(log(insectos), colores) ##
## Kruskal-Wallis rank sum test
##
## data: log(insectos) and colores
## Kruskal-Wallis chi-squared = 16.975, df = 3, p-value = 0.000715Los resultados son exactamente los mismos. No se producen variaciones porque el test de Kruskal-Wallis trabaja sobre rangos, es decir, sobre ordenaciones de los valores de la variable en cada uno de los grupos. Aunque realicemos una transformación logarítmica, el orden entre los valores de la variable se mantiene y por lo tanto la transformación no afecta a los resultados del test.
Si hemos detectado diferencias significativas en la variable respuesta para las distintas poblaciones. ¿Sería posible saber cuáles son los grupos que generan estas diferencias?
library(PMCMR)## Warning: package 'PMCMR' was built under R version 3.4.4## PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.posthoc.kruskal.nemenyi.test(insectos, colores, method = "Chisq")## Warning in posthoc.kruskal.nemenyi.test.default(insectos, colores, method =
## "Chisq"): Ties are present, p-values are not corrected.##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: insectos and colores
##
## amarillo azul blanco
## azul 0.0022 - -
## blanco 0.0039 0.9985 -
## verde 0.4068 0.1878 0.2559
##
## P value adjustment method: none#donde nos salga raya se considera que son similares como en azul y blancoEjercicio
Relaciones sexuales entre universitarios
El departamento de Psicología de una Universidad de Castilla-La Mancha ha realizado un estudio sobre hábitos, preferencias y satisfacción sexual en estudiantes universitarios. Hemos utilizado los datos que recogieron en sus encuestas y queremos conocer si existen diferencias entre la frecuencia mensual de relaciones sexuales de estudiantes universitarios pertenecientes a tres titulaciones universitarias diferentes:
T1: 11 14 7 15 11 13 11 16 10 15 18 12 9 9 10 10 15 10 14 10 10 12 14 12 15 7 13 6 10 15 20 10 13 10 6 14 8 10 8 11
T2: 13 10 12 7 5 10 10 16 9 7 7 2 6 9 9 8 8 10 3 6 5 2 9 3 4 5 10 8 5 9 10 8 13 10 0 2 1 1 0 4
T3: 6 7 3 5 9 6 1 6 0 2 5 6 11 6 7 0 5 7 5 4 7 4 2 8 9 6 1 4 7 7 8 9 7 5 1 6 9 4 7 6
Contesta las siguientes preguntas:
- Introduce los datos en R creando las 2 variables: una que incluya las frecuencias de actividad sexual y otra que será un factor, que nos proporcionará información sobre la titulación universitaria de cada uno de los estudiantes.
actividsex <- c(11, 14, 7, 15, 11, 13, 11, 16, 10, 15, 18, 12, 9, 9, 10, 10, 15, 10, 14, 10, 10, 12, 14, 12, 15, 7, 13, 6, 10, 15, 20, 10, 13, 10, 6, 14, 8, 10, 8, 11,13, 10, 12, 7, 5, 10, 10, 16, 9, 7, 7, 2, 6, 9, 9, 8, 8, 10, 3, 6, 5, 2, 9, 3, 4, 5, 10, 8, 5, 9, 10, 8, 13, 10, 0, 2, 1, 1, 0, 4,6, 7, 3, 5, 9, 6, 1, 6, 0, 2, 5, 6, 11, 6, 7, 0, 5, 7, 5, 4, 7, 4, 2, 8, 9, 6, 1, 4, 7, 7, 8, 9, 7, 5, 1, 6, 9, 4, 7, 6)
titulacion <- as.factor(c(rep(c("T1", "T2", "T3"), each =40)))- Explora los datos de la muestra mediante gráficos y descriptivos. ¿Observamos diferencias en los valores promedios y de variabilidad por grupos?
Media por grupos
## T1 T2 T3
## 11.60 6.90 5.45Desviación estándar por grupos
## T1 T2 T3
## 3.160656 3.848410 2.650109Notamos por el gráfico de caja y bigotes que la variabilidad por grupo es similar, sobre todo entre T1 y T3, además, la media en cada una de las titulaciones es diferente. Los dos datos podemos comprobarlos en los métodos descriptivos.
- Realiza un test F (ANOVA) para comparar las medias de las 3 poblaciones. ¿Cuáles serán las hipótesis nula y alternativa?
Presentamos la tabla ANOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## titulacion 2 826.9 413.4 38.98 1.07e-13 ***
## Residuals 117 1241.1 10.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para decidir sobre el estadístico F, consultamos el valor en tabla:
## [1] 3.073763Por el estadístico F, 38.98>3.073, y por el valor p=1.07e-13, concluimos que la diferencia en media entre cada titulación es estadísticamente significativa. Es decir, rechazamos la hipótesis nula y aceptamos la alternativa: \[ H0: media_{T1}=media_{T2}=media_{T3} \]
\[ H1: media_{Ti} \neq media_{Tj} (\exists i,j \in (1,2,3) ) \]
- Describe los resultados obtenidos indicando quién es el valor del estad;istico de contraste (F), los grados de libertad del factor, los grados de libertad residuales y el valor de P.
Como ya se mencionó, el valor del estadístico F es igual a 39.98, los grados de libertad del factor es 2, del residual es 117 y el valor de P es 1.07e-13.
- También indica quién será el valor crítico de F bajo la hipótesis nula, que nos proporcionará la definición de una región de aceptación y rechazo (consideramos un nivel de significancia alfa = 0.05). El valor es:
## [1] 3.073763- ¿Qué valor obtenemos para la estimación de la varianza común de los datos? Por grupos, obtenemos la varianza:
## T1 T2 T3
## 3.160656 3.848410 2.650109Tras evaluar la tabla ANOVA, ¿cuál será tu conclusión en el contexto del problema? Tras analizar los resultados de la tabla Anova, podemos concluir que existe evidencia de diferencias entre las medias de la frecuencia mensual de cada titulación.
Determina los intervalos de confianza para las medias de frecuencia sexual en cada uno de las titulaciones descritas.
Los intervalos de confianza para cada grupo son:
Titulación T1
Límite inferior
## T1
## 11.34516Límite superior
## T1
## 11.85484Titulación T2
Límite inferior
## T2
## 6.589701Límite superior
## T2
## 7.210299Titulación T3
Límite inferior
## T3
## 5.236321Límite superior
## T3
## 5.663679- Si se han obtenido diferencias significativas entre los grupos, determina dónde se muestran esas diferencias utilizando el test HSD de Tukey. Representa gráficamente las diferencias encontradas e intrepreta los resultados obtenidos.
intervals = TukeyHSD(fm)
intervals## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm(actividsex ~ titulacion))
##
## $titulacion
## diff lwr upr p adj
## T2-T1 -4.70 -6.428861 -2.9711395 0.0000000
## T3-T1 -6.15 -7.878861 -4.4211395 0.0000000
## T3-T2 -1.45 -3.178861 0.2788605 0.1189269
Así podemos ver que hay diferencias en la media entre las tres titulaciones.
- A partir de los residuos obtenidos y mediante procedimientos gráficos, confirma la validez del modelo comprobando el cumplimiento de las 3 condiciones asumidas al hacer el ANOVA. ¿Crees que el modelo ANOVA utilizado es válido para estos datos?
Independencia
Normalidad
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en cada una de las titulaciones:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.900 -1.900 0.400 0.000 2.175 9.100
Por el test de Shapiro-Wilk, no tenemos evidencia suficiente para rechazar la normalidad de los residuos.
##
## Shapiro-Wilk normality test
##
## data: fm$residuals
## W = 0.98473, p-value = 0.1945Homocedasticidad
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en las titulaciones:
Comparando la desviación máxima con la mínima obtenemos una orientación sobre la homocedasticidad
## [1] 1.45217El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales)
##
## Bartlett test of homogeneity of variances
##
## data: fm$residuals by titulacion
## Bartlett's K-squared = 5.3467, df = 2, p-value = 0.06902- Si NO se verificaran estas condiciones, ¿hay alguna prueba no paramétrica para abordar los datos? Pruébalo y comenta los resultados obtenidos.
Kruskal-Wallis y pruebas post-hoc
Ho: la variable respuesta es la misma en todas las poblaciones valoradas
Ha: la variable respuesta es mayor en alguna de las poblaciones
Utilizaremos la prueba no paramétrica de Kruskal-Wallis
kruskal.test(actividsex, titulacion)##
## Kruskal-Wallis rank sum test
##
## data: actividsex and titulacion
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12Obtenemos el cuantil buscado:
qchisq(0.05, 3-1, lower.tail = F)## [1] 5.991465Valores del estadístico > 5.991465 estarán incluidos en la región de rechazo. En nuetro caso 51.504 es mucho mayor que el valor crítico obtenido.
¿Sería posible saber cuáles son los grupos que generan estas diferencias?
posthoc.kruskal.nemenyi.test(actividsex, titulacion, method = "Chisq")##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: actividsex and titulacion
##
## T1 T2
## T2 1.5e-06 -
## T3 1.4e-11 0.14
##
## P value adjustment method: none#donde nos salga raya se considera que son similares como en azul y blancoAsí, concluímos que hay evidencia estadísticamente significativa para asegurar que las medias son distintas en cada una de las titulaciones.