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. En todo el esquema de las cosas, somos realmente tan diferentes?

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. Mas varianza mas diferente la muestra

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:

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("blue", "green", "white","yellow"), ylab = "Número de insectos atrapados") #el signo de la ñ porque es factor #amarillo tiene dos datos atipicos

#
tapply(insectos, colores, mean) #aplique la media en tales cosas = tapply
## amarillo     azul   blanco    verde 
## 47.16667 14.83333 15.66667 31.50000

ANOVA y pruebas post-hoc

Esta es la forma de pedir un ANOVA en R:

fm = aov( lm(insectos ~ colores) ) #modelo lineal que busca la relacion entre las variables (tipo continua y tipo factor)

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

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) #qf nivel de significacia n-1 y el i-1 y si tiene cola inferior se pone F 
## [1] 3.68232

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 (nivel de confianza.. como tiene dos colas dividimos para dos, y el grado de libertad)
sp <- sqrt(46)  #desviación t??pica de la varianza muestral común #46 media
ee  <- valor_t * (sp/ sqrt(6))  #error de estimación #6 es el n de cada grupo
media
## [1] 47.16667

El intervalo de confianza me sirve para hacer predicciones

l?mite superior del intervalo de confianza de la media de insectos capturados para las trampas amarillas

media + ee 
## [1] 48.57826

l?mite inferior del intervalo de confianza de la media de insectos capturados para las trampas amarillas

media - ee 
## [1] 45.75507

Si 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) #Prueba de Turkey  #cercano a 1 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?

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)

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.8333
boxplot(fm$residuals)

hist(fm$residuals)#para ver si se paproxima a la campana ed gauuss

qqnorm(fm$residuals)

El test de Shapiro-Wilk ca que no tenemos evidencia suficiente para rechazar la hip?tesis nula (normalidad de los residuos)

shapiro.test(fm$residuals)  #p valor es mayor que 0,05 es normal 
## 
##  Shapiro-Wilk normality test
## 
## data:  fm$residuals
## W = 0.97337, p-value = 0.75

Homocedasticidad

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","green"))  #La varianza de los grupos se la obtiene a traves de los residuos

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

Comparando la desviaci?n m?xima con la m?nima obtenemos una orientaci?n sobre la falta de homocedasticidad (>2 aproximadamente) # extrictamente >2 no hay hocedasticidad

max(desviaciones) / min(desviaciones)    
## [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) #p > 0,05 varianzas diferentes.. se rechaza la hipotesis 
## 
##  Bartlett test of homogeneity of variances
## 
## data:  fm$residuals by colores
## Bartlett's K-squared = 5.2628, df = 3, p-value = 0.1535

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? anova multivariante?

Ho: la variable respuesta es la misma en todas las poblaciones valoradas (esto es cuando tenemos mas de un variable independiente)

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) #acepto la alternativa
## 
##  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. el p valor se contrasta con el valor teorico es con el chi-cuadrado

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.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. aqui el valor empirico es chi-cuadrado y el teorico 5.991 se rechaza la hipotesis nula

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")
## 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

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:

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

 titulaciones <- as.factor(c(rep(c("T1", "T2", "T3"), each =40)))
boxplot(freq ~ titulaciones,col = c("yellow", "blue","green"), ylab = "Frecuencia mensual de relaciones sexuales") 

tapply(freq, titulaciones, mean) 
##    T1    T2    T3 
## 11.60  6.90  5.45
anov = aov( lm(freq ~ titulaciones) )
summary(anov)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## titulaciones   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
anov
## Call:
##    aov(formula = lm(freq ~ titulaciones))
## 
## Terms:
##                 titulaciones Residuals
## Sum of Squares      826.8667 1241.1000
## Deg. of Freedom            2       117
## 
## Residual standard error: 3.256945
## Estimated effects may be unbalanced

Ho: Las medias son iguales H1: Las medias son difeentes

F value = 38.98 gl factor = 2 gl residuals = 117 Pr(>F) = 1.07e-13

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. Asi obtenemos el cuantil buscado:

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

Valores del estadistico > 4.098172 estar?n incluidos en la regi?n de rechazo. En nuetro caso 38.98 es mucho mayor que el valor cr?tico obtenido. Se rechaza la Hip?tesis nula.

media <- mean(freq[titulaciones =="T1"]) 
valor_t <- pt(0.05/2, 80-2) 
sp <- sqrt(10.6)  
ee  <- valor_t * (sp/ sqrt(40))
media
## [1] 11.6
media +ee 
## [1] 11.86251
media -ee 
## [1] 11.33749
media <- mean(freq[titulaciones =="T2"]) 
valor_t <- pt(0.05/2, 80-2) 
sp <- sqrt(10.6)  
ee  <- valor_t * (sp/ sqrt(40))
media
## [1] 6.9
media +ee 
## [1] 7.162508
media -ee 
## [1] 6.637492
media <- mean(freq[titulaciones =="T3"]) 
valor_t <- pt(0.05/2, 80-2) 
sp <- sqrt(10.6)  
ee  <- valor_t * (sp/ sqrt(40))
media
## [1] 5.45
media +ee 
## [1] 5.712508
media -ee 
## [1] 5.187492

T1: lim sup 11.86251 lim inf 11.33749 T2: lim sup 7.162508 lim inf 6.637492 T3: lim sup 5.712508 lim inf 5.187492

intervals = TukeyHSD(anov) #Prueba de Turkey  #cercano a 1 estadisticamente iguales
intervals
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm(freq ~ titulaciones))
## 
## $titulaciones
##        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)

El cero ‘’0’’ se encuentra dentro de la comparacioon T3 y T2 por tanto son similares.

Independencia

plot(anov$residuals)

### Normalidad

summary(anov$residuals)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -6.900  -1.900   0.400   0.000   2.175   9.100
boxplot(anov$residuals)

hist(anov$residuals)

qqnorm(anov$residuals) 
qqline(anov$residuals)

shapiro.test(anov$residuals) 
## 
##  Shapiro-Wilk normality test
## 
## data:  anov$residuals
## W = 0.98473, p-value = 0.1945

No tenemos evidencia suficiente para rechazar la hip?tesis nula ### Homocedasticidad

boxplot(anov$residuals~titulaciones, col = c("red", "black", "pink")) 

desviaciones <- tapply(anov$residuals, titulaciones, sd)
max(desviaciones) / min(desviaciones) 
## [1] 1.45217

Si hay homocedasticidad.

bartlett.test(anov$residuals ~ titulaciones)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  anov$residuals by titulaciones
## 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? Pruebalo y comenta los resultados obtenidos.
kruskal.test(freq, titulaciones) 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  freq and titulaciones
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12
library(PMCMR)
posthoc.kruskal.nemenyi.test(freq, titulaciones, method = "Chisq")
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  freq and titulaciones 
## 
##    T1      T2  
## T2 1.5e-06 -   
## T3 1.4e-11 0.14
## 
## P value adjustment method: none
qchisq(0.05, 2-1, lower.tail = F)
## [1] 3.841459

51.504 > 3.841 S e rechaza la hip?tesis nula.