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.
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:
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") #para tipo factor, se generan diagramas de caja que especifico la variable factor
tapply(insectos, colores, mean) #una funcion a varios vectores usamos tapply
## amarillo azul blanco verde
## 47.16667 14.83333 15.66667 31.50000
Esta es la forma de pedir un ANOVA en R:
fm = aov( lm(insectos ~ colores) ) #lo guardamos un un objeto, aplicamos el modelo lineal, propio de R
Pedimos un resumen de la tabla del ANOVA
summary(fm) #si fueran iguales me darian cercano a uno pero como da casi cero no son iguales, si son diferentes las cantidades de insectos que van por color
## 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
#Fvalue: valor de la muestra
Elementos 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:
qf(0.05, 3-1, 18-3, lower.tail = F) #para tener la prueba de post hoc, si realmente las medias son diferentes, especicamos le nivel de significacion, el I-1, y si hay valores en la cola
## [1] 3.68232
#ggrados de libertad, no se toma el grupo 1
#y 18-3, porque no se toma un grupo
#valor teorico que obtengo, comparo con el F value
Valores 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) #prueba t, dos colar se divide le intervalo de confianza, 18 es le n total, n-1-(I-1)
sp <- sqrt(46) #desviación típica de la varianza muestral común, el 46 me da la tabla del anova en sq, residuales
ee <- valor_t * (sp/ sqrt(6)) #error de estimación, 6 es el n de cada grupo
media #intervalos de confianza con la media y estadistico P
## [1] 47.16667
límite superior del intervalo de confianza de la media de insectos capturados para las trampas amarillas
media + ee #solo superior
## [1] 48.57826
límite inferior del intervalo de confianza de la media de insectos capturados para las trampas amarillas
media - ee #solo inferior
## [1] 45.75507
Si hemos detectado diferencias significativas entre las medias de las poblaciones. ¿Sería posible saber cuáles son los grupos que generan estas diferencias?
intervals = TukeyHSD(fm) #son diferentes estadisticamente si son cercanos a 0 y si son a uno son estadisticamente iguales
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
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ías?
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.
plot(fm$residuals) # residuo buenos predicciones buenas, si hay un patron correlacionados pero si estan dispersos no estan correlacionados
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.8333
boxplot(fm$residuals) #supondria que hay homocedastidad porque la distancia a los cuartiles es semejante, pero hay un dato atipico
hist(fm$residuals) # que se asemeje a una campana de gauss
qqnorm(fm$residuals)
qqline(fm$residuals) #dato quantil-quantil
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) #test a los residuos y me da un pa valor, si el pvalor >0.05 son normales, si es <0.05 no son normales
##
## Shapiro-Wilk normality test
##
## data: fm$residuals
## W = 0.97337, p-value = 0.75
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:
boxplot(fm$residuals~colores, col = c("yellow", "blue", "white","pink")) #residuales me dan las diferencia entre el datos real y el ajustado, la varianza de los grupos se la obtiene a traves de los residuos
desviaciones <- tapply(fm$residuals, colores, sd) #sacar desviacion por grupo y ocupamos tapply
Comparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente)
max(desviaciones) / min(desviaciones) #maximo de la desviacion estandar, >2 las varianzas diferentes, si es <2 si hay homocedasticidad, varianzas iguales
## [1] 2.980357
El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales)
bartlett.test(fm$residuals ~ colores) #nos da un pvalor>0.05 no hay homocedasticidad, si pvalor<0.05 si hay homocedasticidad
##
## Bartlett test of homogeneity of variances
##
## data: fm$residuals by colores
## Bartlett's K-squared = 5.2628, df = 3, p-value = 0.1535
¿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) #analisis multivariante
##
## Kruskal-Wallis rank sum test
##
## data: insectos and colores
## Kruskal-Wallis chi-squared = 16.975, df = 3, p-value = 0.000715
Indica 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) #16.975 es mayor a 5.99 por lo cual se rechaza la hipotesis nula y se acepta la hipotesis alternativa
## [1] 5.991465
Valores 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.000715
Los 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)
## 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") #ayuda a saber cual se debe considerar como un solo (-) =no hay diferencia siginificativa se le podria considerar como uno solo (azul,blanco)
## 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
#libreria muestra cuales son iguales y cuales diferentes/ analisis por detras de todos el analisis de anova
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:
frecuencias<-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)
universidades<-as.factor(c(rep(c("T1","T2","T3"),each=40)))
Diagrama de caja
Media de las frecuencias
## T1 T2 T3
## 11.60 6.90 5.45
Ho: Las medias de las frecuencias son iguales
H1: Al menos una media de las frecuencias es diferente
ANOVA
Tabla se resumen ANOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## universidades 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 ' ' 1
Elementos generados en ANOVA
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "contrasts" "xlevels" "call" "terms"
## [13] "model"
Valor estadistico de contraste: 38.98
Grados de libertad: 2
Grados de libertad residuales: 117
Valor de P: 0.1945
## [1] 3.963472
Valores del estadístico > 3.963472 estarán incluidos en la región de rechazo. En nuetro caso 38.98 es mucho mayor que el valor crítico obtenido, por lo tanto se rechaza la Ho y se acepta H1.
## [1] 134.56
Se puede concluir que las medias de las frecuencias de actividad sexual de los estudiantes de las diferentes titulaciones, son diferentes. Lo que se comprueba al comparar el valor F obtenido en la tabla de ANOVA con el obtenido.
Titulacion 1
media1 <- mean(frecuencias[universidades =="T1"])
valor_t1 <- pt(0.05/2, 80 - 2)
sp1 <- sqrt(10.6)
ee1 <- valor_t * (sp1/ sqrt(40))
media1
## [1] 11.6
Limite superior del intervalo de confianza
## [1] 11.86244
Limite inferior del intervalo de confianza
## [1] 11.33756
Titulacion 2
media2<- mean(frecuencias[universidades =="T2"])
valor_t2 <- pt(0.05/2, 80 - 2)
sp2 <- sqrt(10.6)
ee2 <- valor_t * (sp2/ sqrt(40))
media2
## [1] 6.9
Limite superior del intervalo de confianza
## [1] 7.16244
Limite inferior del intervalo de confianza
## [1] 6.63756
Titulacion 3
media3<- mean(frecuencias[universidades =="T3"])
valor_t3 <- pt(0.05/2, 80 - 2)
sp3 <- sqrt(10.6)
ee3 <- valor_t * (sp3/ sqrt(40))
media3
## [1] 5.45
Limite superior del intervalo de confianza
## [1] 5.71244
Limite inferior del intervalo de confianza
## [1] 5.18756
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm(frecuencias ~ universidades))
##
## $universidades
## 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
Grafico de los intervalos
Validacion del modelo ANOVA
1.Independencia
2.Normalidad
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.900 -1.900 0.400 0.000 2.175 9.100
Diagrama de caja
Histograma
Grafico Q-Q Normal
Test Shapiro-Wilk
##
## Shapiro-Wilk normality test
##
## data: ej$residuals
## W = 0.98473, p-value = 0.1945
3. Homocedasticidad
Diagrama de caja
Desviaciones
## [1] 1.45217
Como 1.45217<2, existe homocedasticidad (varianzas iguales)
Test Barlett
Ho=Las varianzas son iguales (Homocedasticiidad)
H1=Las varianzas son diferentes
##
## Bartlett test of homogeneity of variances
##
## data: ej$residuals by universidades
## Bartlett's K-squared = 5.3467, df = 2, p-value = 0.06902
Como 0.065902>0.05 se rechaza las hipotesis alternativa por lo tanto no hay homocedasticidad (varianzas diferentes)
Si, es posible utilizar la prueba no paramétrica Kruskal-Wallis: ¿hay diferencias significativas entre las poblaciones
Ho=Las muestras de las poblaciones provienen de la misma poblacion
H1=Las muestras de las poblaciones provienen de diferentes problaciones
##
## Kruskal-Wallis rank sum test
##
## data: frecuencias and universidades
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12
## [1] 3.841459
Como el valor 51.504>3.841459, se rechaza Ho y se acepta H1, por lo tanto las muestras provienen de poblaciones diferentes.
##
## Kruskal-Wallis rank sum test
##
## data: log(frecuencias) and universidades
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12
##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: frecuencias and universidades
##
## T1 T2
## T2 1.5e-06 -
## T3 1.4e-11 0.14
##
## P value adjustment method: none