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:

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:

Organización de los datos

Generamos dos variables: frecuencia_sexual es la variable respuesta y universidades es la variable factor (establece los grupos de interés):

frecuencia <- 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)))

Exploramos los datos de la muestra:

Realiazamos digramas de caja y bigote para las 3 universidades

boxplot(frecuencia ~ universidades, col = c("red", "purple", "white"), ylab = "Frecuencia mensual de relaciones sexuales")

El valor de las medias de las Universidades es:

tapply(frecuencia, universidades, mean)
##    T1    T2    T3 
## 11.60  6.90  5.45

Los datos obtenidos nos permiten establecer que la mayor distribución de datos en la universidad T1 esta en la parte superior, tomando en cuenta que su media es de 11.60 relaciones sexuales cada mes. En tonor a la universidad T2 se puede observar que la dispersion esta en un sentido casi neutral es decir existe una distribucion similar en tanto arriba de la media como bajo de la misma, es importante recalcar que la media es de 6.90 relaciones sexules cada mes y esta es inferior a la media de la primera universidad. La ultima grafica muestra que la distribucion de la mayoria de los datos se encuentra por debajo de la media, la media es de 5.45 relaciones sexuales cada mes y es la menor de las tres universidades.

ANOVA y pruebas post-hoc

Esta es la forma de pedir un ANOVA en R:

fm = aov( lm(frecuencia ~ universidades) )#busca la relacion entre las variables

Resumen de la tabla del ANOVA

summary(fm)#estadisticamente significativo
##                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 el ANOVA:

names(fm)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

H0: Las medias de las tres universidades son iguales

H1: Al menos una media de las universidades es diferente

Valor estadistico de constraste: 38.98

Grados de libertad: 2

Grados de libertad residuales: 117

Estadístico de contraste F

qf(0.05, 2-1, 80-2, lower.tail = F)
## [1] 3.963472

Debido a que el valor de contraste es mayor el Estadistico de contraste, se dice que son estadisticamente iguales.

El intervalo de confianza de la media

Universidad T1, con un nivel de confianza del 95%:

mediaT1 <- mean(frecuencia[universidades =="T1"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueba de dos colas
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40)) 
mediaT1
## [1] 11.6
mediaT1 + ee 
## [1] 11.86251
mediaT1 - ee 
## [1] 11.33749

Universidad T2, con un nivel de confianza del 95%:

mediaT2 <- mean(frecuencia[universidades =="T2"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueba de dos colas
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40)) 
mediaT2
## [1] 6.9
## [1] 7.162508
mediaT2 - ee 
## [1] 6.637492

Universidad T3, con un nivel de confianza del 95%:

mediaT3 <- mean(frecuencia[universidades =="T3"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueba de dos colas
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40)) 
mediaT3
## [1] 5.45
mediaT3 + ee 
## [1] 5.712508
mediaT2 - ee 
## [1] 6.637492

Si hemos detectado diferencias no tan significativas entre las medias de las poblaciones. Por eso para asegurar unos resultados veraces se procede a hacer:

Prueba de Turkey

intervals = TukeyHSD(fm)
intervals 
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm(frecuencia ~ 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
plot(intervals)

Debido a que los valores son cercanos a 0 se puede interpretar que son valores diferentes, es decir las medias de las tres universidades son diferentes entre si.

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

plot(fm$residuals)

Ya que los datos estan muy dispersos se puede establecer que no estan corelacionados.

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.100

Homocedasticidad = igual varianza porque la distancia de los cuartiles es semenjante aunque hay un dato tipico

boxplot(fm$residuals)

Existe un solo valor atipico en la parte superior.

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.1945

Debido a que el p-valor es mayot a 0.05, se dice que existe una normalidad en los residuos.

Homocedasticidad

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:

boxplot(fm$residuals~universidades, col = c("blue", "white","green"))  

desviaciones <- tapply(fm$residuals, universidades, sd)

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)   
## [1] 1.45217

Debido a que el valor calculado no es ayor que dos, se dice que existe Homocedasticidad entre las varibles, relacionando al grafico, entre T2 y T3 esta la mayor cantidad de relacion entre sus varianzas.

Test de Bartlett

Indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales)

bartlett.test(fm$residuals ~ universidades)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  fm$residuals by universidades
## Bartlett's K-squared = 5.3467, df = 2, p-value = 0.06902

El resultado del test de Barlet es que no existe homocedasticidad, debido a que el p-valor es un poco mayor a 0.05, pero debemos analisar que en el anterior se dijo que es homocedastico, esto es debido a que la relacion entre la universidad T1 y T2 difiere por un poco y genera esta contradiccion.

Kruskal-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(frecuencia, universidades)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  frecuencia and universidades
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12

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, 2-1, lower.tail = F)
## [1] 3.841459

Valores del estadístico > 3.841459 estarán incluidos en la región de rechazo. En nuetro caso 51.504 es mucho mayor que el valor crítico obtenido.Acepto la hipotesis alternativa

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-Wallis rank sum test
## 
## data:  log(frecuencia) and universidades
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12

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)
posthoc.kruskal.nemenyi.test(frecuencia, universidades, method = "Chisq") #ayuda a saber cual se debe considerar como un solo (-) = no hay diferentes significativa se le podria considerar como un solo (azul, blanco). Libreria muestra cuales son iguales y cuales diferentes/analisis por detras de todos los analisis de anova
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  frecuencia and universidades 
## 
##    T1      T2  
## T2 1.5e-06 -   
## T3 1.4e-11 0.14
## 
## P value adjustment method: none