ANOVA

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

T1<-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)
T2<-c(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<-c(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)
length(T1)

## [1] 40

length(T2)

## [1] 40

length(T3)

## [1] 40

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.

actividad <- 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)))
titulacion

##   [1] T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1
##  [24] T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T1 T2 T2 T2 T2 T2 T2
##  [47] T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2
##  [70] T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T2 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3
##  [93] T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3 T3
## [116] T3 T3 T3 T3 T3
## Levels: T1 T2 T3

actividad

##   [1] 11 14  7 15 11 13 11 16 10 15 18 12  9  9 10 10 15 10 14 10 10 12 14
##  [24] 12 15  7 13  6 10 15 20 10 13 10  6 14  8 10  8 11 13 10 12  7  5 10
##  [47] 10 16  9  7  7  2  6  9  9  8  8 10  3  6  5  2  9  3  4  5 10  8  5
##  [70]  9 10  8 13 10  0  2  1  1  0  4  6  7  3  5  9  6  1  6  0  2  5  6
##  [93] 11  6  7  0  5  7  5  4  7  4  2  8  9  6  1  4  7  7  8  9  7  5  1
## [116]  6  9  4  7  6

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

boxplot(actividad ~ titulacion, col = c("yellow", "blue","green"), ylab = "Número de Relaciones Sexuales ")

tapply(actividad,titulacion,mean)

##    T1    T2    T3 
## 11.60  6.90  5.45

SI SE EVIDENCIA DIFERENCIA ENTRE LOS VALORES PROMEDIOS Y SU VARIABILIDAD.

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

h0:Las medias de cada proyecto de titulacion son iguales h1:Las medias entre las titulaciones difieren

fm = aov( lm(actividad ~ titulacion) )

Pedimos un resumen de la tabla del ANOVA

summary(fm)

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

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

El valor de F será:

summary(fm)[[1]][1,4]

## [1] 38.97486

Los grados de libertad del factor:

#fm$rank-1
summary(fm)[[1]][1,1]

## [1] 2

Los grados de libertad del residuo:

#fm$df.residual
summary(fm)[[1]][2,1]

## [1] 117

El valor de P

summary(fm)[[1]][1,5]

## [1] 1.067578e-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).

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.963472 estarán incluidos en la región de rechazo. En nuestro caso 38.98 es mucho mayor que el valor crítico obtenido.

• ¿Qué valor obtenemos para la estimación de la varianza común de los datos?

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.

media<- summary(fm)[[1]][2,3]
sp <- sqrt(media)  #desviación típica de la varianza muestral común
sp

## [1] 3.256945

• Tras evaluar la tabla ANOVA, ¿cuál sería tu conclusión en el contexto del problema?

TRAS LA APLICACION DEL ANALISIS ANOVA, SE INFIERE QUE: NO EXISTE EVIDENCIA ESTADISTICA SUFICIENTE PARA ACEPTAR LA HIPOTESIS NULA. ES DECIR SE SUGIERE RECHAZAR LA HIPOTESIS NULA LO QUE IMPLICA QUE EXISTE DIFERENCIA ENTRE LAS MEDIAS DE LAS TITULACIONES

• Determina los intervalos de confianza para las medias de frecuencia sexual en cada uno de las titulaciones descritas.

La raíz cuadrada de la media de los cuadrados del error, además se utiliza en la obtención de los intervalos de confianza de las medias en cada uno de los grupos de interés.

media <- mean(actividad[titulacion =="T1"]) 
valor_t <- pt(0.05/2, 80-2) 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación 
media

## [1] 11.6

límite superior del intervalo de confianza de la media de actividad sexual para la titulacion T1

media + ee 

## [1] 11.8626

límite inferior del intervalo de confianza de la media de actividad sexual para la titulacion T1

media - ee 

## [1] 11.3374

media <- mean(actividad[titulacion =="T2"]) 
valor_t <- pt(0.05/2, 80-2) 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación 
media

## [1] 6.9

límite superior del intervalo de confianza de la media de actividad sexual para la titulacion T2

media + ee 

## [1] 7.162603

límite inferior del intervalo de confianza de la media de actividad sexual para la titulacion T2

media - ee 

## [1] 6.637397

media <- mean(actividad[titulacion =="T3"]) 
valor_t <- pt(0.05/2, 80-2) 
ee  <- valor_t * (sp/ sqrt(40))  #error de estimación 
media

## [1] 5.45

límite superior del intervalo de confianza de la media de actividad sexual para la titulacion T3

media + ee 

## [1] 5.712603

límite inferior del intervalo de confianza de la media de actividad sexual para la titulacion T3

media - ee 

## [1] 5.187397

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

intervals = TukeyHSD(fm)
intervals

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

plot(intervals)

SE EVIDENCIA LAS DIFERENCIAS ENTRE LAS TITULACIONES, ES DECIR DIFERENCIAS ENTRE LAS MEDIAS

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

Los supuestos que se deben cumplir son tres: independencia, homocedasticidad y normalidad.

Independencia

plot(fm$residuals)

Homocedasticidad

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

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

desviaciones <- tapply(fm$residuals, titulacion, 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 ~ titulacion)

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

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

boxplot(fm$residuals)

hist(fm$residuals)

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)

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

COMO SE CUMPLEN LAS 3 CONDICIONES, ES VALIDA LA APLICACION DEL ESTUDIO ANOVA

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

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

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

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)

## [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 de la actividad sexual, ¿variarían los resultados del test estadístico?

kruskal.test(log(actividad), titulacion) 

## 
##  Kruskal-Wallis rank sum test
## 
## data:  log(actividad) 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.