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))) #each= cada 6 corresponde al color

Exploramos los datos de la muestra: Generalmente en ANOVA se usa el diagrama de caja y bigotes

boxplot(insectos ~ colores, col = c("yellow", "blue", "white","green"), ylab = "Número de insectos atrapados")

#insectos en funcion de colores, en funcion de la variable tipo FACTOR 
tapply(insectos, colores, mean) #aplicar una funcion a varios vectores
## amarillo     azul   blanco    verde 
## 47.16667 14.83333 15.66667 31.50000
#mapply a la matriz
#lapply a la lista

ANOVA y pruebas post-hoc

Esta es la forma de pedir un ANOVA en R:

fm = aov( lm(insectos ~ colores) ) #insectos es continua colores es de tipo factor
#aov es la funcion anova
fm
## Call:
##    aov(formula = lm(insectos ~ colores))
## 
## Terms:
##                  colores Residuals
## Sum of Squares  4218.458   920.500
## Deg. of Freedom        3        20
## 
## Residual standard error: 6.784173
## Estimated effects may be unbalanced

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
#como tiene los asteriscos y la probabilidad de que no sea significativa es baja, entonces colores es significativa , esto dice que el color si influye en que los insectos se acerquen.

#Evalue es el chi^2 experimental que debe ser comparado con el valor teorico

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)
## [1] 3.68232
#nivel de significancia, n-1, I-1, se especifica si tiene cola 
#sigue una chi^2

#este valor que se dale de aqui es el valor teorico, debe compararse con el experimental

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) #se quiere con el 0.95 de confianza, se divide 0.05/2 ... tamaño de la muestra - nº de grupos

sp <- sqrt(46)  #desviación típica de la varianza muestral común, el 46 es del anova del mean sq

ee  <- valor_t * (sp/ sqrt(6))  #error de estimación 
#el 6 es el nº de elementos de cada grupo

media
## [1] 47.16667
#intervalo de confianza necesita de la media, la varianza y los parametros de arriba

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

media + ee  #media + error de estiamcion (para el liminite superior)
## [1] 48.57826

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

media - ee  #media - error de estiamcion (para el liminite 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)
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
#entre la media d azul y amarillo el p es 0
#SON DIFERENTES ESTADISTICAMENTE SI EL P ES CERCANO A 0
# SON IGUALES ESTADISTICAMENTE SI EL P ES CERCANO A 1
plot(intervals)

#si el intervalo contiene al 0 estos grupos son estadisticamente iguales

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)

#Mas dispersa la grafica es mejor, esto significa que los residuos estan menos relacionados, por lo tanto se cumple la condicion de independencia

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)

#si se aproxima a una campana de gauss
qqnorm(fm$residuals) 
qqline(fm$residuals)

#Es normal cuando casi todos estan sbre la linea normal, que no se generen los numeros indicando la posicion del punto, significa que es normal

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.97337, p-value = 0.75
#si el pvalor es > a 0.05 se puede afirmar que es normal la distribucion

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

#LOS RESIDUALES DAN LA DIFERENCIA ENTRE DEL DATO REAL Y EL AJUSTADO
#Aqui se debe fijar el ancho de la caja
#el verde tiene mayor dispersion por los bigotes

#la varianza de los grupos se oobtiene a traves de los residuos
desviaciones <- tapply(fm$residuals, colores, sd)
#Desviacion estandar por grupo, sd se obtiene a partir de los residuos <- desviacion estandar

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
#si la diferencia del maximo de la sd de los residuos para el minico de la sd de los residuos es >2 se dice que no hay homocedasticidad

#En este caso no hay homocedasticidad

El test de Bartlett para la homocedasticidad-> indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas 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
#si el pvalor es >0.05 se rechaza la hipotesis nula de que son iguales, es decir son diferentes
#no hay homocedasticidad

Kruskal-Wallis y pruebas post-hoc

Es un 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
#como 16.975 es mayor que el teorico 5.99, se rechaza la H0

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
#si se aplica logaritmos no deben variar los datos
#Se aplican los logartimos para minimizar la escala de los datos

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
#cuando salen las rayias es que no hay diferencias significativas, se compruba homocesaticidad, etc
#cuales son iguales o diferentes

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:

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)
titulo <- as.factor(c(rep(c("fisico","quimico","sociales"), each =40)))
boxplot(frecuencia ~ titulo, col = c("violet", "light blue", "light green"), ylab = "Frecuencia sexual por mes")

tapply(frecuencia, titulo, mean)
##   fisico  quimico sociales 
##    11.60     6.90     5.45

H0:Las medias de frecuencia sexual para los 3 grupos de titulación son iguales H1: Al menos una de las medias de frecuencia sexual para los 3 grupos de titulación es diferente

anova = aov( lm(frecuencia ~ titulo) )
anova
## Call:
##    aov(formula = lm(frecuencia ~ titulo))
## 
## Terms:
##                    titulo Residuals
## Sum of Squares   826.8667 1241.1000
## Deg. of Freedom         2       117
## 
## Residual standard error: 3.256945
## Estimated effects may be unbalanced
summary(anova)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## titulo        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

Valor F teórico

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

El valor F empírico es mayor que el valor F teórico, por lo que se rechaza la H0.

valor_t <- pt(0.05/2, 80 - 2)
valor_t
## [1] 0.5099406
sp <- sqrt(10.6)
sp
## [1] 3.255764
ee  <- valor_t * (sp/ sqrt(40))
ee
## [1] 0.262508
  1. Límite superior del intervalo de confianza de la media de frecuencia sexual para los estudiantes fisicos
media1 <- mean(frecuencia[titulo =="fisico"])
media1 + ee
## [1] 11.86251
  1. Límite inferior del intervalo de confianza de la media de frecuencia sexual para los estudiantes fisicos
media1 - ee
## [1] 11.33749
  1. Límite superior del intervalo de confianza de la media de frecuencia sexual para los estudiantes químicos
media2 <- mean(frecuencia[titulo =="quimico"])
media2 + ee
## [1] 7.162508
  1. Límite inferior del intervalo de confianza de la media de frecuencia sexual para los estudiantes químicos
media2 - ee
## [1] 6.637492
  1. Límite superior del intervalo de confianza de la media de frecuencia sexual para los estudiantes sociales
media3 <- mean(frecuencia[titulo =="sociales"])
media3 + ee
## [1] 5.712508
  1. Límite inferior del intervalo de confianza de la media de frecuencia sexual para los estudiantes químicos
media3- ee
## [1] 5.187492

El valor F empírico es mayor que el valor F teórico, por lo que se rechaza la H0.Por lo tanto, las medias de la frencuencia sexual para los grupos de titulacion son significativamente diferentes.

intervals = TukeyHSD(anova)
intervals
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm(frecuencia ~ titulo))
## 
## $titulo
##                   diff       lwr        upr     p adj
## quimico-fisico   -4.70 -6.428861 -2.9711395 0.0000000
## sociales-fisico  -6.15 -7.878861 -4.4211395 0.0000000
## sociales-quimico -1.45 -3.178861  0.2788605 0.1189269
#SON DIFERENTES ESTADISTICAMENTE SI EL P ES CERCANO A 0
# SON IGUALES ESTADISTICAMENTE SI EL P ES CERCANO A 1
plot(intervals)

Independencia

plot(anova$residuals)

Se cumple la condicion de independencia debido a la dispersión de puntos

Normalidad

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

hist(anova$residuals)

qqnorm(anova$residuals) 
qqline(anova$residuals)

No existe dispersión de puntos, se cumple con la condición de normalidad

Homocedasticidad

Shapiro Test

shapiro.test(anova$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  anova$residuals
## W = 0.98473, p-value = 0.1945
boxplot(anova$residuals~titulo, col = c("violet","light blue", "light green"))

Determinación de normalidad

desviaciones <- tapply(anova$residuals,titulo, sd)
desviaciones
##   fisico  quimico sociales 
## 3.160656 3.848410 2.650109
max(desviaciones) / min(desviaciones)   
## [1] 1.45217

Diferencia entre desviaciones es menor a 2, por lo tanto, existe homocedasticidad