ANOVA
Daniela Gutierrez
28 de septiembre de 2018
Idea intuitiva del ANOVA
Ejercicios
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. variables: actividad (frecuancia de actividad sexual) y titulacion (titulación universitaria de cada uno de los estudiantes).
actividad <- 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? Sí, 11.60 6.90 5.45 para TI, T2 y T3 respectivamente.
boxplot(actividad ~ titulacion, col = c("yellow", "green", "white","green"), ylab = "Frecuencia mensual de relaciones sexuales")
tapply(actividad, titulacion, mean)## T1 T2 T3
## 11.60 6.90 5.45- Realiza un test F (ANOVA) para comparar las medias de las 3 poblaciones. ¿Cuáles serían las hipótesis nula y alternativa? nula: medias de las poblaciones son iguales, T1=T2=T3.
alternativa: al menos una de las poblaciones es diferente (T1).
ANOVA y pruebas post-hoc
ANOVA en R:
fm = aov( lm(actividad ~ titulacion) )Resumen de la tabla del ANOVA
summary(fm)## 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 ' ' 1Elementos 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"Describe los resultados obtenidos indicando quién es el valor del estadístico de contraste (F), los grados de libertad del factor, los grados de libertad residuales y el valor de P. F= 38.98 grados de libertad (n-I) 117 grado de libertad (I-1): 3-1=2 P= 1.07e-13
También indica quién sería 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 significación alfa = 0.05).
qf(0.05, 2-1, 80-2, lower.tail = F)## [1] 3.963472##########################################Valores del estadístico > 3.089203 estarán incluidos en la región de rechazo. En nuetro caso 38.98 es mucho mayor que el valor crítico obtenido.
¿Qué valor obtenemos para la estimación de la varianza común de los datos?
Tras evaluar la tabla ANOVA, ¿cuál sería tu conclusión en el contexto del problema? Se cumple la hipótesis alternativa, al menos una de las muestras es diferente.
Determina los intervalos de confianza para las medias de frecuencia sexual en cada uno de las titulaciones descritas.
Intervalo de confianza de la media de T1, con un nivel de confianza del 95%:
media <- mean(actividad[titulacion =="T1"])
valor_t <- pt(0.05/2, 80 - 2)
sp <- sqrt(10.6) #desviación típica de la varianza muestral común
ee <- valor_t * (sp/ sqrt(40)) #error de estimación
media## [1] 11.6límite superior del intervalo de confianza
media + ee ## [1] 11.86251límite inferior del intervalo de confianza
media - ee ## [1] 11.33749T2:
media <- mean(actividad[titulacion =="T2"])
valor_t <- pt(0.05/2, 80 - 2)
sp <- sqrt(10.6) #desviación típica de la varianza muestral común
ee <- valor_t * (sp/ sqrt(40)) #error de estimación
media## [1] 6.9límite superior del intervalo de confianza de la media
media + ee ## [1] 7.162508límite inferior del intervalo de confianza de la media
media - ee ## [1] 6.637492T3:
media <- mean(actividad[titulacion =="T3"])
valor_t <- pt(0.05/2, 80 - 2)
sp <- sqrt(10.6) #desviación típica de la varianza muestral común
ee <- valor_t * (sp/ sqrt(40)) #error de estimación
media## [1] 5.45límite superior del intervalo de confianza de la media
media + ee ## [1] 5.712508límite inferior del intervalo de confianza de la media
media - ee ## [1] 5.187492- 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(actividad ~ 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.1189269plot(intervals)
- 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
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.
## -6.900 -1.900 0.400 0.000 2.175 9.100boxplot(fm$residuals)
hist(fm$residuals)
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.98473, p-value = 0.1945Homocedasticidad
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:
boxplot(fm$residuals~titulacion, col = c("yellow", "green", "white")) 
- Si NO se verificaran estas condiciones, ¿hay alguna prueba no paramétrica para abordar los datos? Pruébalo y comenta los resultados obtenidos.
Prueba no paramétrica Kruskal-Wallis: ¿hay diferencias significativas entre las poblaciones?
kruskal.test(actividad, titulacion)##
## Kruskal-Wallis rank sum test
##
## data: actividad and titulacion
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12Indica 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 51.504 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(actividad), titulacion) ##
## Kruskal-Wallis rank sum test
##
## data: log(actividad) and titulacion
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12Los 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(actividad, titulacion, method = "Chisq")## Warning in posthoc.kruskal.nemenyi.test.default(actividad, titulacion,
## 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: actividad and titulacion
##
## T1 T2
## T2 1.5e-06 -
## T3 1.4e-11 0.14
##
## P value adjustment method: none