Análisis de Varianza

El Análisis de Varianza ANOVA permite comparar las medias de distintas poblaciones. El Análisis ANOVA prueba la igualdad de medias por medio de la comparación de distintas varianzas(variaciones).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. Para saber cual o cuales medias son distintas se realizan las pruebas PostHoc.

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:

Introduce los datos en R creando las 2 variables: una que incluya las frecuencias de actividad sexual y otra que será un factor, que nos proporcionará información sobre la titulación universitaria de cada uno de los estudiantes.

fsexual<-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)
titulacion <- as.factor(c(rep(c("T1", "T2", "T3"), each =40)))

Explora los datos de la muestra mediante gráficos y descriptivos. ¿Observamos diferencias en los valores promedios y de variabilidad por grupos?

boxplot(fsexual ~ titulacion, col = heat.colors(3), ylab = "Frecuencia sexual",main="Diagrama de Caja",col.main="darkorange2")

Hay diferencias en los valores promedios o medias en los tres grupos, las varianzas y el ancho de cada caja es similar en los tres grupos, además no existen datos atípicos.

tapply(fsexual, titulacion, mean)

##    T1    T2    T3 
## 11.60  6.90  5.45

Aquí se observan los valores de las medias para las tres titulaciones que son la varibale factor, las medias son distintas lo cual demuestra que no se pueden agrupar los factores.

Realiza un test F (ANOVA) para comparar las medias de las 3 poblaciones. ¿Cuáles serían las hipótesis nula y alternativa?

Ho: las medias son iguales en las tres titulaciones. Ha: hay alguna media distinta.

Pedir un anova

sr = aov( lm(fsexual ~ titulacion) )

El modelo lineal busca la relación entre las variables, fsexual es una variable de tipo continua y titulación variable tipo factor. Cuando se aplica el modelo de regresión lineal lm se genera el análisis anova(aov)

Resumen tabla ANOVA

summary(sr)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## titulacion    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

#Cuando hay *** la variable es significativa en cuanto a las medias, es decir no son iguales y el valor se acerca a cero (P).

La variable titulacion es significativa ya que Pr(>F) tiende a cero, es decir las medias no son iguales.

Elementos generados en el ANOVA:

names(sr)

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

#Muestra los argumentos con los que se  genera el modelo.

Describe los resultados obtenidos indicando quién es el valor del estadístico de contraste (F), los grados de libertad del factor(titulación), los grados de libertad residuales y el valor de P.

Estadístico de contaste: 38.98 G.L de Titulación: 2 G.L.Residuales: 117 Valor P:1.07e-13

También indica quién sería el valor crítico de F bajo la hipótesis nula, que nos proporcionará la definición de una región de aceptación y rechazo (consideramos un nivel de significación alfa = 0.05).

qf(0.05, 2-1, 80-2, lower.tail = F)

## [1] 3.963472

#Usamos qf para tener la prueba(N.C, I-1, n-1)
#Si F<Fcrit: Ho no se rechaza
#Si F>Fcrit: Ho se rechaza

Con un nivel de confianza del 95% podemos deducir que la frecuencia sexual de los estudiantes de las tres titulaciones es distinta, se rechaza la hipótesis nula ya que F=38.98 > Fcrit= 3.96.

¿Qué valor obtenemos para la estimación de la varianza común de los datos? Mean sq=10.6
Tras evaluar la tabla ANOVA, ¿cuál sería tu conclusión en el contexto del problema? Al realizar el análisis ANOVA se conluye que las medias son distintas con un niverl de significancia del 5% donde se rechaza la hipótesis nula. Ahora necesitamos saber por que se rechaza y esto lo podemos lograr mediante las pruebas postHoc para saber cuales medias son distintas.
Determina los intervalos de confianza para las medias de frecuencia sexual en cada uno de las titulaciones descritas. Titulacion 1 Intervalo de confianza de la media de estudiantes de la titulacion 1, con un nivel de confianza del 95%:

media1 <- mean(fsexual[titulacion =="T1"]) 
valor_t1 <- pt(0.05/2, 80 - 2) 
sp1 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee1  <- valor_t1 * (sp1/ sqrt(40))  #error de estimación (6) es el n de cada grupo.
media1

## [1] 11.6

límite superior del intervalo de confianza de la media de estudiantes de la titulacion 1

media1 + ee1

## [1] 11.86251

límite inferior del intervalo de confianza de la media de estudiantes de la titulacion 1

media1 - ee1

## [1] 11.33749

Titulacion 2 Intervalo de confianza de la media de estudiantes de la titulacion 2, con un nivel de confianza del 95%:

media2 <- mean(fsexual[titulacion =="T2"]) 
valor_t2 <- pt(0.05/2, 80 - 2) 
sp2 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee2  <- valor_t2 * (sp2/ sqrt(40))  #error de estimación (6) es el n de cada grupo.
media2

## [1] 6.9

límite superior del intervalo de confianza de la media de estudiantes de la titulacion 2

media2 + ee2

## [1] 7.162508

límite inferior del intervalo de confianza de la media de estudiantes de la titulacion 2

media2 - ee2

## [1] 6.637492

Titulacion 3

Intervalo de confianza de la media de estudiantes de la titulacion 3, con un nivel de confianza del 95%:

media3 <- mean(fsexual[titulacion =="T3"]) 
valor_t3 <- pt(0.05/2, 80 - 2) 
sp3 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee3  <- valor_t3 * (sp3/ sqrt(40))  #error de estimación (6) es el n de cada grupo.
media3

## [1] 5.45

límite superior del intervalo de confianza de la media de estudiantes de la titulacion 3

media3 + ee3

## [1] 5.712508

límite inferior del intervalo de confianza de la media de estudiantes de la titulacion 3

media3 - ee3

## [1] 5.187492

Si se han obtenido diferencias significativas entre los grupos, determina dónde se muestran esas diferencias utilizando el test HSD de Tukey. Representa gráficamente las diferencias encontradas e intrepreta los resultados obtenidos.

intervals1 = TukeyHSD(sr)
intervals1

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

#Nos ayuda a ver entre cuales hay diferencia.
#son diferentes si p adj es cercano a cero, si es cercano a uno son iguales.

plot(intervals1, col="purple2")

#Si un intervalo contiene al cero son estadisticamente iguales.

Al rechazar la Hipótesis Nula (Ho) se aplican las pruebas de comparación multiple, como la prueba de Tukey que permite determinar si dos medias poblacionales son significativamente distintas entre si. Al realizar esta prueba vemos que T3-T2 son estadísticamente iguales, mientras que T3-T1 y T2-T1 son significativamente distintas ya que el p valor es cercano a cero.

A partir de los residuos obtenidos y mediante procedimientos gráficos, confirma la validez del modelo comprobando el cumplimiento de las 3 condiciones asumidas al hacer el ANOVA. ¿Crees que el modelo ANOVA utilizado es válido para estos datos?

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

color_1 <-colorRampPalette(c("yellow ", "blue", "yellow"))

plot(sr$residuals, main = "Prueba de independencia" ,pch=20,cex = 2, col=color_1(120), ylab = "Residuos", xlab = " ")

#fm  guardamos el anova.
#La grafica mas dispersa mejor, hay menor correlacion, cumplen con uno de los supuestos de los residuos.

Podemos observar en la gráfica que la mayoría de datos se encuentran dispersos, pero hay algunos que están agrupados, entonces se puede concluir que la correlación de los residuos no es tán buena.

Normalidad

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos: Ya sabemos que las medias son diferentes. Ahora la varianza es diferente?

summary(sr$residuals)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -6.900  -1.900   0.400   0.000   2.175   9.100

#Hay un dato que no se ajusta, los cuantiles son casi del mismo ancho.
#Homocedastidad: igual varianza.

boxplot(sr$residuals, main="Diagrama de Caja Residuos",col="goldenrod1" )

color_2 <-colorRampPalette(c("orchid1", "darkviolet", "dodgerblue"))
hist(sr$residuals, main = "Histograma de Residuos", ylab = "Frecuencia", xlab = "Residuales", col=color_2(14))

En la primera gráfica, del diagrama de caja se identifica un valor atípico y los cuartiles son casi del mismo ancho. En la segunda gráfica se observa que los datos siguen una distribución normal ya que la gráfica se aproxima a una campana de Gauss.

qqnorm(sr$residuals) 
qqline(sr$residuals)

#Grafico cuantil cuantil, cuando en verdad son normales. No hay datos atípicos y casi están sobre la linea normal.

shapiro.test(sr$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  sr$residuals
## W = 0.98473, p-value = 0.1945

#Se aplica el test de Shapiro los residuos si p valora =>0.05 son normales y si <=0.05 no son normales.

El test de Shapiro-Wilk indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (normalidad de los residuos), el p valor es 0.1945 es mayor a 0.05, hay evidencia de normalidad.

Homocedasticidad

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

boxplot(sr$residuals~titulacion, col = c("olivedrab1", "olivedrab2", "olivedrab3"))

#Diferencia, la varianza a partir de los residuos. 
#nos fijamos en el ancho de las cajas. 
#La varianza de los grupos se obtiene a través de los residuos.

desviaciones1 <- tapply(sr$residuals, titulacion, sd)
#Desv std por grupo.

Comparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente)

max(desviaciones1) / min(desviaciones1)

## [1] 1.45217

#si este valor es >2 no hay homocedasticidad, la varianza no es igual.

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

bartlett.test(sr$residuals ~ titulacion)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  sr$residuals by titulacion
## Bartlett's K-squared = 5.3467, df = 2, p-value = 0.06902

#Test de Bartlett para homocedasticidad. 
#p valor es > a 0.05 se rechaza la hipotesis nula, no hay homocedasticidad.

Con la prueba de Barlett determinamos la falta de homocedasticidad con un p valor de 0.069 mayor a 0.05 las varianzas no son iguales.

Se puede concluir que el modelo ANOVA no es el adecuado ya que solo se cumplen dos supuestos: independencia y normalidad.

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

Kruskal-Wallis y pruebas post-hoc (Test Multivariante)

La alternativa paramétrica del ANOVA es comparar distribuciones de varias muestras independientes.

¿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(fsexual, titulacion)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  fsexual and titulacion
## Kruskal-Wallis chi-squared = 51.504, df = 2, p-value = 6.547e-12

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, 2-1, lower.tail = F)#Valor teorico

## [1] 3.841459

#chi cuadrado es 16.975 >5.99 se rechaza la Ho.

El valor del estadístico es mayor que el valor teórico y esta incluido en la región de rechazo. En nuestro caso 51.504 es mucho mayor que el valor crítico obtenido(3.48). La frecuencia sexual en las tres titulaciones es distinta.

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(fsexual), titulacion) #log para reducir escala de dispersión

## 
##  Kruskal-Wallis rank sum test
## 
## data:  log(fsexual) and titulacion
## 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) #Sirve para comparar. 
posthoc.kruskal.nemenyi.test(fsexual, titulacion, method = "Chisq")# Hace todo lo de arriba pero internamente.

## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  fsexual and titulacion 
## 
##    T1      T2  
## T2 1.5e-06 -   
## T3 1.4e-11 0.14
## 
## P value adjustment method: none

#sale una - cuando se debe considerar que deben ser uno solo. Como grupos tanto la var y med son iguales.
#Comprueba homocedasticidad, los residuos.