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.

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("yellow", "blue", "white","green"), ylab = "Número de insectos atrapados")

tapply(insectos, colores, mean)
## 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) )

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

el valor de F y Pr se complementan

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"

Cuando se lica se generan 13 elemtos, los demas son argumenrtos con los que se genrera el modelo y dan los mismso que una regresion lienal *** significa que es significativo y que se acerca a 0

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)
## [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. el valor teorico se usa para compararlo con la tabla ANOVA

¿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%: la prueba t se le da una valor de alpha, esta a los extremos y se divide para 2 en los extremos,

media <- mean(insectos[colores =="amarillo"]) 
valor_t <- pt(0.05/2, 18 - 3)  #depende del grado de confianza que tengamos 
sp <- sqrt(46)  #desviación típica de la varianza muestral común, el 46 sale del ANOVA 
ee  <- valor_t * (sp/ sqrt(6))  #error de estimación, nos ayuda a sacar el intervalo de confianza 
media
## [1] 47.16667

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 amarills

media - ee 
## [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? si son cercanos a cero son estadisticamente iguales y los que no son estadisticamente diferentes

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

vValidació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)

residuales la grafica mientras mas dispersa es mejor, por lo que se encuentran menos correlacionados ya que los residuos no deben estar correlacionados y si estos estuvieran correlacionados el ANOVA no esta dando informacion estadisticamente creible

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

para que los residuos sean normales aplicamos graficos descriptivos

boxplot(fm$residuals)

se puede ver que tenemos un dato que no se ajusta y los cuantiles son del mismo ancho entonces hay una misma variabilidad, homocelasticidad significa la misma varianza

hist(fm$residuals)

el hisitograma nos dira si se aproxma a una campana de Gauss, estaria mal si estuviera en una recta

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) Nunca se debe elimianr un test de shapiro-wilk p valor mayor o igual a 0,75 son normales sino no lo son

shapiro.test(fm$residuals)
## 
##  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: grafuico de caja a los residuales da, el dato real y el dato ajustado se lo ve como la varianza a partir de los residuos la carianza de los grupos se la obtiene atravess de los residuos

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

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)

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), la hpotesis nula en esta son iguales

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

Kruskal-Wallis y pruebas post-hoc

anova de mas de una variable ¿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.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)
## [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? para aplicar un modelo solemos aplicar transformaciones, por lo que se suelen aplicar los logartmos para quitar la dispersion

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?

usamos la libreria PMCMR

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

la rayita implica que no son significativos entre ellos por lo que se debe considerar como uno solo, hace todo el analisis por detras y confirma que al usar todos los anaisis de ANOVA llega a esa conclusion

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:

Relaciones <- 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)
carreras <- as.factor(c(rep(c("T1", "T2", "T3"), each =40)))
boxplot(Relaciones ~ carreras, col = c("pink", "purple", "orange"), ylab = "Frecuencia de Relaciones Sexuales")

tapply(Relaciones, carreras, mean)
##    T1    T2    T3 
## 11.60  6.90  5.45

Se observa que la media en el grupo de la primera carrera tiene una media mayor de relaciones sexuales en relacion a las carreras 2 y 3.

Hipotesis nula: La frecuencia de relaciones sexuales = tipo de carrera Hipotesis alternativa: La frecuencia de relaciones sexuales es diferente en relacion al tipo de carrera

fm = aov( lm(Relaciones ~ carreras) )
summary(fm)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## carreras      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

Df: son los grados de libertad

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

El valor estadistico de contraste F: 38.98 grados de libertad del factor: 2 que corresponde a las carreras grados de libertad residuales: 117 p-valor:1.07e-13 Como el p-valor es muy pequeño se concluye que hay diferencias muy significativas entre las carreras.

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

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

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.

Para la Carrea T3

media <- mean(Relaciones[carreras =="T3"]) 
valor_t <- pt(0.05/2, 120 - 2)  #depende del grado de confianza que tengamos 
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común, el 46 sale del ANOVA 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación, nos ayuda a sacar el intervalo de confianza 
media
## [1] 5.45

El limite superior para el analis de Frecuencia de Relaciones Sexuales en carrera T3

media + ee 
## [1] 5.712514

El limite iferior para el analsis de Frecuencia de Relaciones Sexuales en la carrera T3

media - ee 
## [1] 5.187486

Para la Carrera T2

media2 <- mean(Relaciones[carreras =="T2"]) 
valor_t <- pt(0.05/2, 120 - 2)  #depende del grado de confianza que tengamos 
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común, el 46 sale del ANOVA 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación, nos ayuda a sacar el intervalo de confianza 
media2
## [1] 6.9

El limite superior para el analis de Frecuencia de Relaciones Sexuales en carrera T2

media2 + ee 
## [1] 7.162514

El limite iferior para el analsis de Frecuencia de Relaciones Sexuales en la carrera T2

media2 - ee 
## [1] 6.637486

Para la Carrera T1

media1 <- mean(Relaciones[carreras =="T1"]) 
valor_t <- pt(0.05/2, 120 - 2)  #depende del grado de confianza que tengamos 
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común, el 46 sale del ANOVA 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación, nos ayuda a sacar el intervalo de confianza 
media1
## [1] 11.6

El limite superior para el analis de Frecuencia de Relaciones Sexuales en carrera T1

media1 + ee 
## [1] 11.86251

El limite iferior para el analsis de Frecuencia de Relaciones Sexuales en la carrera T1

media1-ee
## [1] 11.33749
intervals = TukeyHSD(fm)
intervals
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm(Relaciones ~ carreras))
## 
## $carreras
##        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)

Los dos intervalos T3-T1 y T2-T1 no contienen un cero por lo que son significativamente diversas, al existir un cero en T3-T2 estas no son significativamente diversas.

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

para que los residuos sean normales aplicamos graficos descriptivos

boxplot(fm$residuals)

Como se observa en la grafica existe una misma variablilidad entre los cuartiles.

hist(fm$residuals)

La grafica presenta un histograma adecuado ya que se presenta una variación entre los tamaños de las barras.Este se aproxima a una campana de Gauss

qqnorm(fm$residuals) 
qqline(fm$residuals)

La grafica muestra una dispersion de datos normal

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

Al realizar el analisis del test de shapiro-wilk Se rechazará la hipótesis nula de normalidad si el estadístico W es menor que el valor crítico proporcionado por la tabla, por lo que existe presencia de normalidad que se valida con la grafica antes mostrada. ### Homocedasticidad Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos: grafico de caja a los residuales da, el dato real y el dato ajustado se lo ve como la varianza a partir de los residuos la carianza de los grupos se la obtiene atravess de los residuos

boxplot(fm$residuals~carreras,col =c("pink", "orange", "purple"))  

desviaciones <- tapply(fm$residuals, carreras,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

En base a la comparacion >2, se ha obtenido presencia de homocedasticidad

  • Si NO se verificaran estas condiciones, ¿hay alguna prueba no paramétrica para abordar los datos? Pruébalo y comenta los resultados obtenidos.

Las pruebas no parametricas aplicables son Kruskal-Wallis y pruebas post-hoc

L as hipotesis que se plantean

Ho: la frecuencia de Relaciones Sexuales es la misma en todas las carreras estudiadas

Ha: la frecuencia de relaciones Sexuales mayor en la carrera T1

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

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? para aplicar un modelo solemos aplicar transformaciones, por lo que se suelen aplicar los logartmos para quitar la dispersion

kruskal.test(log(Relaciones), carreras) 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  log(Relaciones) and carreras
## 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