ANOVA

library(PMCMR)

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.

v1<-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)
v2<-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)
v3<-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)
frec_as<-c(v1,v2,v3)
titulaciones<- 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?

Exploramos los datos de la muestra:

boxplot(frec_as ~ titulaciones, col = c("yellow", "blue", "white"), ylab = "Frecuencia Relaciones Sexuales")

Luego las medias son:

tapply(frec_as, titulaciones, mean)

   T1    T2    T3 
11.60  6.90  5.45 

Con el gráfico anterior podemos notar que las varianzas son diferentes, además se podría decir que las 3 medias son diferentes, pero entre “T1” y “T2”, parece que tienen una media estadísticamente semejante.

---

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

Vamos a probar si \[H_0:\mu_{T1}=\mu_{T2}=\mu_{T3}\], es decir si las tres medias so iguales.

anova1 = aov( lm(frec_as~ titulaciones) )#Realizamos el ANOVA

Pedimos un resumen de la tabla del ANOVA

summary(anova1)

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

Elementos generados en el ANOVA:

names(anova1)

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

Debido a que el valor p es 1.07e-13, que es menor que el nivel de significancia de 0.05, puede rechazar la hipótesis nula y concluir que algunas de las titulaciones tienen diferentes medias.

---

* 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, 40*3-2, lower.tail = F)

[1] 3.921478

---

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

Grados de libertad factor:2 Grados de libertad Residuales:117 Valor P:1.07e-13

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.

---

* ¿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, se utiliza en la obtención de los intervalos de confianza de las medias en cada uno de los grupos de interés.

[1] 10.6

---

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

Se rechaza la hipotesis nula y se concluye que existe diferencia entre la actividad sexual de los estudiantes de cada titulación.

---

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

Intervalo de confianza de la media de las relaciones sexuales para la titulación T1, con un nivel de confianza del 95%:

media1 <- mean(frec_as[titulaciones=="T1"]) 
valor_t1 <- pt(0.05/2, 117 - 2) 
sp1 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee1  <- valor_t1 * (sp1/ sqrt(6))  #error de estimación 
media1

[1] 11.6

límite superior del intervalo de confianza de la media

media1 + ee1

[1] 12.27781

límite inferior del intervalo de confianza de la media

media1 - ee1

[1] 10.92219

Intervalo de confianza de la media de las relaciones sexuales para la titulación T2, con un nivel de confianza del 95%:

media2 <- mean(frec_as[titulaciones=="T2"]) 
valor_t2 <- pt(0.05/2, 117 - 2) 
sp2 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee2  <- valor_t2 * (sp2/ sqrt(6))  #error de estimación 
media2

[1] 6.9

límite superior del intervalo de confianza de la media

media2 + ee2

[1] 7.577806

límite inferior del intervalo de confianza de la media

media2 - ee2

[1] 6.222194

Intervalo de confianza de la media de las relaciones sexuales para la titulación T3, con un nivel de confianza del 95%:

media3 <- mean(frec_as[titulaciones=="T3"]) 
valor_t3 <- pt(0.05/2, 117 - 2) 
sp3 <- sqrt(10.6)  #desviación típica de la varianza muestral común
ee3  <- valor_t3 * (sp3/ sqrt(6))  #error de estimación 
media3

[1] 5.45

límite superior del intervalo de confianza de la media

media3 + ee3

[1] 6.127806

límite inferior del intervalo de confianza de la media

media3 - ee3

[1] 4.772194

---

* 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(anova1)
intervals1

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = lm(frec_as ~ titulaciones))

$`titulaciones`
       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(intervals1)

Se puede ver que son muy diferentes entre todos los grupos.

---

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

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(anova1$residuals)

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

summary(anova1$residuals)

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

boxplot(anova1$residuals)

hist(anova1$residuals)

qqnorm(anova1$residuals) 
qqline(anova1$residuals)

El test de Shapiro-Wilk indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (normalidad de los residuos)

shapiro.test(anova1$residuals)

    Shapiro-Wilk normality test

data:  anova1$residuals
W = 0.98473, p-value = 0.1945

Como obtenemos un valor p mayor a 0.05, podemos concluir que los datos se distribuyen normalmente.

Homocedasticidad

boxplot(anova1$residuals~titulaciones, col = c("yellow", "blue", "white")) 

* 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

kruskal.test(frec_as, titulaciones)

    Kruskal-Wallis rank sum test

data:  frec_as and titulaciones
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. Por lo tanto se concluye que provienen de poblaciones diferentes.