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

Exploramos los datos de la muestra:

boxplot(insectos ~ colores, col = c("yellow", "blue", "white","green"), ylab = "Número de insectos atrapados") #el grafico de cajs para ANOVA, Y VER MEDIAS, #de diagrams de la variable factor,buscamos un anova pra ver q son estasiticamente diferntes

tapply(insectos, colores, mean)#aplicar media a una data frame o un conjuento de vectores`#aul y blano se puede hacer en el msmog rupo pero
## 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) ) #continua en funcion de

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 de la h nula)

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)#prueba de pos thoc"qf"funcion,se toma el cuerpo de la densi
## [1] 3.68232
#f valor amyor que teorico, estaditicamente difernetes

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) #prueab T, depende del grado de confianza
sp <- sqrt(46)  #desviación típica de la varianza muestral común,46 esti de la sd en anov  arri
ee  <- valor_t * (sp/ sqrt(6))  #error de estimación , 6 n de cada grupo
media#PREDICIION DEL NUMERO DE INSECTOS DE UN COLOR A APRTI DE LA MEUSTRA DADA
## [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 amarillas

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?

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) # si cercano a 0 son diferentse, si son aprox a 1 no son tan diferentes

Si un GRUPO CONTIENE AL CEO SON ESTASDITIACAMENT 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?

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) #Si lso redifuos son buenos, las predicciones son buenas, mientra mas disperso mejor, porque no esstaria correlacioonados

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)#homoselaticidad =varianza

hist(fm$residuals)

qqnorm(fm$residuals) 
qqline(fm$residuals)#outlets, no hay datos atipos, siguen dist 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) #si el p valor es emnor o igual a 0.05, son normaoles, depende 
## 
##  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"))  #residual es la diferencai del dato real y dato ajustado, la varianza de lso grupos e sobtienen a traves d elos residuos, mas dispersion en el verde

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) #no hay hocedaticidad, las varinzas no osn iguales    
## [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) #no hay evidencia suficiente para rechazar l hiop nula,mayor a 0,05, rechamos la hiptesis nula de quee las varinzas son iguales xq es 0.15, mayor a 0,05 no hay homocedaticidad
## 
##  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?

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)#ANOVA mulñtivariable
## 
##  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. si es mayor a 0,05 aceptamos la ha.

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)#numero de grupos -1
## [1] 5.991465
#este es el valor teorico y el valor de la muestra de chi cuadrao es 16.975, rechamos h0 porque es mayor 1u3 5, 999

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) #log aplicamos para quitar la dispersion
## 
##  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")#rayitas para las qeu se deb considerar en uno solo, por al media y la varaianza. Esa libreria hace todo el anlissi anterior
## 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 HACER ESTOS EJERCICIOS COMO DE DEBER====

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:

ESTABLESCO VARIABLES

frecuencia_sex<- 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)
factor_titul<- as.factor(c(rep(c("T1","T2","T3"),each=40)))
boxplot(frecuencia_sex ~ factor_titul, col = c("yellow", "blue", "white"), ylab = "Frecuencia sexual")

tapply(frecuencia_sex, factor_titul, mean)
##    T1    T2    T3 
## 11.60  6.90  5.45

ANOVA APLICADO

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

HIPÓTESISI NULA=todas las varianzas son iguales HIPÓTESISI ALTERNATIVA=hay una varinza diferente, se acepta la Ha porque Pr(>F)=1.07e-13 menor 0.05

F_value: 38.98 grados de libertad: 2 grados de libertad residuales: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. Así obtenemos el cuantil buscado:

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

Valores del estadístico > 3.96347 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.

sp <- sqrt(10.6)
sp
## [1] 3.255764

Intervalo de confianza de la media de lA FRECUENCIA SEXUAL para T1, con un nivel de confianza del 95%:

media <- mean(frecuencia_sex[factor_titul =="T1"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueab T, depende del grado de confianza
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación , 40 n de cada grupo
media#PREDICIION DE LA FRECUENCIA SEX DE UNa titulacion A APRTI DE LA MEUSTRA DADA
## [1] 11.6
#LIMITE SUPERIOR
lS1<-media + ee 
lS1
## [1] 11.86251
#LIMITE INFERIOR
lI1<-media - ee
lI1
## [1] 11.33749

Intervalo de confianza de la media de lA FRECUENCIA SEXUAL para T2, con un nivel de confianza del 95%:

media <- mean(frecuencia_sex[factor_titul =="T2"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueab T, depende del grado de confianza
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación , 40 n de cada grupo
media#PREDICIION DE LA FRECUENCIA SEX DE UNa titulacion A APRTI DE LA MEUSTRA DADA
## [1] 6.9
#LIMITE SUPERIOR
lS1<-media + ee 
lS1
## [1] 7.162508
#LIMITE INFERIOR
lI1<-media - ee
lI1
## [1] 6.637492

Intervalo de confianza de la media de lA FRECUENCIA SEXUAL para T3, con un nivel de confianza del 95

media <- mean(frecuencia_sex[factor_titul =="T3"]) 
valor_t <- pt(0.05/2, 80 - 2) #prueab T, depende del grado de confianza
sp <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación , 40 n de cada grupo
media#PREDICIION DE LA FRECUENCIA SEX DE UNa titulacion A APRTI DE LA MEUSTRA DADA
## [1] 5.45
#LIMITE SUPERIOR
lS1<-media + ee 
lS1
## [1] 5.712508
#LIMITE INFERIOR
lI1<-media - ee
lI1
## [1] 5.187492

GRUPOS QUE GENERAN LAS DIFERENCIAS

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

GRAFICO DE DIFERENCIAS

plot(intervals)

INTERPRETACION: Si el valor de p adjunto se aproxima a 0 los grupos son diferentes, y si el valor se aproxima a 1 no son tan diferentes, y se puede consiedera cono un solo grupo. T2 y T1 son diferente T3 y T1 son difernetes T2 y T3 son semejantes estadísticamente

A partir de los residuos del modelo comprobaremos si el modelo ANOVA es adecuado. Los supuestos que se deben cumplir son tres: independencia, normalidad y homocedasticidad.

INDEPENDENCIA

plot(fm$residuals) #Si lso redifuos son buenos, las predicciones son buenas, mientra mas disperso mejor, porque no esstaria correlacioonados

NORMALIDAD

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

#HISTOGRAMA
hist(fm$residuals)

#VERIFICACION DE DATOS ATIPICOS
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)  #si el p valor es mEnor o igual a 0.05, son normales.
## 
##  Shapiro-Wilk normality test
## 
## data:  fm$residuals
## W = 0.98473, p-value = 0.1945

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

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

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

El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales)

bartlett.test(fm$residuals ~ factor_titul) #no hay evidencia suficiente para rechazar l hiop nula,mayor a 0,05, rechamos la hiptesis nula de quee las varinzas son iguales xq es 0.15, mayor a 0,05 no hay homocedaticidad
## 
##  Bartlett test of homogeneity of variances
## 
## data:  fm$residuals by factor_titul
## Bartlett's K-squared = 5.3467, df = 2, p-value = 0.06902

Si, la prueba de Kruskall-Wallis:

kruskal.test(frecuencia_sex, factor_titul)#ANOVA multivariable
## 
##  Kruskal-Wallis rank sum test
## 
## data:  frecuencia_sex and factor_titul
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12
#
qchisq(0.05, 2-1, lower.tail = F)
## [1] 3.841459

Este es el valor teorico y el valor de la muestra de chi cuadrao es 51.505, rechazamos h0 porque,el valor experimental es mayor que el teórico.

EJERCICIO 2 ====

¿Hay diferencias de salario en función del lugar de trabajo?

Nos gustaría saber si el salario varía en alguna de las 3 provincias del Ecuador. Para ello, se realizó un estudio con 50 personas por provincia a las que se preguntó su salario en euros por semana.

Quito: 299 313 300 321 308 312 300 310 281 308 309 300 303 303 311 308 291 298 276 290 310 308 295 310 286 295 289 293 291 297 297 287 297 302 298 301 313 290 306 313 294 308 295 303 316 299 313 296 290 299

Cuenca: 252 248 232 229 256 233 240 237 248 232 230 246 236 250 238 243 245 241 235 249 238 231 230 239 261 243 242 245 249 258 245 236 244 242 229 246 244 244 255 247 236 252 237 259 248 237 236 252 236 239

Guayaquil: 272 268 285 274 278 287 297 275 269 281 270 284 282 281 280 286 265 283 281 272 269 286 268 288 284 282 304 280 283 281 281 286 287 288 278 272 268 287 269 272 270 271 291 265 280 280 275 294 269 277

Contesta las siguientes preguntas: