Apuntes de R

Apunte 09: ANOVA

Author

David Torres Irribarra e Isidora Naranjo López

Published

May 16, 2024

La distribución F

La prueba ANOVA sigue una lógica similar a las pruebas revisadas anteriormente (Z y t). En ANOVA evaluamos la probabilidad de observar diferencias de medias observadas bajo la hipótesis nula usando la distribución F.

La distribución F depende de dos parámetros, ambos grados de libertad: el primero asociado al numerador (variabilidad inter grupo) y el segundo asociado al denominador (variabilidad intra grupo).

\[ F = \frac{ \frac{\text{Suma de cuadrados inter}}{\text{Grados de libertad inter}} }{ \frac{\text{Suma de cuadrados intra}}{\text{Grados de libertad intra}} } \]

Una distribución F solo puede tomar valores de 0 o más y su forma cambia considerablemente dependiendo de los distintos grados de libertad.

Por ejemplo, una distribución F con 2 grados de libertad inter y 5 grados de libertad intra es así:

Mientras que una distribución F con 6 grados de libertad inter y 30 grados de libertad intra es así:

Y una distribución F con 40 grados de libertad inter y 200 grados de libertad intra es así:

Al realizar una prueba ANOVA, R automáticamente calcula los grados de libertad correspondientes, calcula el estadígrafo F y obtiene el valor p correspondiente a la probabilidad de observar valores F mayores e igual en la distribución correspondiente. Por ejemplo, si obtenemos un valor F igual a 1,9 con 5 grados de libertad inter y 80 grados de libertad intra, el valor p se calcularía de esta manera:

Lo que nos indica que el observar una variabilidad inter grupo al menos 1,9 veces más grande que la variabilidad intragrupo en una distribución \(F_{(5,80)}\) ocurre en un 10,4% de las veces.

Ejemplo de ANOVA: Autoestima entre distintos cursos

Usamos ANOVA para evaluar la hipótesis nula de que las medias poblacionales desde las que surgen los grupos son iguales. Si rechazamos la hipótesis nula, entonces tenemos evidencia que, de manera global, la membresía en los distintos grupos se relaciona con el valor de la variable cuantitativa que estamos estudiando.

Queremos ver si existen diferencias significativas entre como se comportan los promedios en la variable autoestima entre distintos cursos.

#Cargamos la base de datos
poli <- read.csv("https://david-ti.github.io/introstats/data/datos_poli_2017.csv")

Recordamos que la variable se compone de 10 preguntas:

  • AU1 Siento que soy una persona valiosa, al menos igual que los demás

  • AU2 Siento que tengo cualidades positivas

  • AU3 En general, tiendo a sentir que soy un fracaso

  • AU4 Soy capaz de hacer las cosas tan bien como la mayoría de las otras personas

  • AU5 Siento que no tengo mucho de lo que sentirme orgulloso

  • AU6 Tengo una actitud positiva hacia mí mismo

  • AU7 Considerando todas las cosas, estoy satisfecho conmigo mismo

  • AU8 Me gustaría tener más respeto conmigo mismo

  • AU9 Me siento inútil a veces

  • AU10 Algunas veces pienso que no soy bueno en absoluto

#Invertimos los puntajes de aquellas preguntas en donde corresponde
poli$au3i <- (5 + 1)-poli$au3
poli$au5i <- (5 + 1)-poli$au5
poli$au8i <- (5 + 1)-poli$au8
poli$au9i <- (5 + 1)-poli$au9
poli$au10i <- (5 + 1)-poli$au10

# Creamos un objeto para almacenar las 10 preguntas
poli$au <- NA
poli$au <- poli[,c("au1","au2","au3i","au4","au5i","au6","au7","au8i","au9i","au10i")]

# Calculamos el puntaje de autoestima para cada persona de la base de datos
poli$puntaje_au <- rowSums(poli$au)

Ahora, veamos los grupos que compararemos en nuestro análisis, es decir, la variable curso.

# Para facilitar el trabajo con la variable la convertimos en factor y asignamos nombres a las categorías
poli$curso = factor(poli$curso, 
                    levels=c("1","2","3","4","5"), 
                    labels=c("Septimo", "Octavo", "PrimeroM", "SegundoM", "TerceroM"))
table(poli$curso)

 Septimo   Octavo PrimeroM SegundoM TerceroM 
    4435     4418     3706     3648     3477 
barplot(table(poli$curso), 
        ylab = 'Frecuencia', 
        xlab = 'Cursos')

Veremos que pasos son relevantes para llevar a cabo el análisis:

  1. Ya que hemos generado la variable autoestima, ahora queremos examinar si los estudiantes en distintos cursos tienen distintos niveles de autoestima. Para esto podemos comenzar usando describeBy()
# Activamos la libraría
library(psych)

# describeBy nos muestra los datos descriptivos estadísticos para cada uno de los cursos en relación a la variable autoestima
describeBy(poli$puntaje_au, poli$curso)

 Descriptive statistics by group 
group: Septimo
   vars    n  mean  sd median trimmed  mad min max range  skew kurtosis   se
X1    1 3958 34.65 7.4     35   34.87 7.41  11  50    39 -0.29    -0.11 0.12
------------------------------------------------------------ 
group: Octavo
   vars    n  mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 3986 34.53 7.66     35   34.78 7.41  10  50    40 -0.31    -0.09 0.12
------------------------------------------------------------ 
group: PrimeroM
   vars    n  mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 3378 34.23 7.71     35   34.46 7.41  11  50    39 -0.26    -0.33 0.13
------------------------------------------------------------ 
group: SegundoM
   vars    n  mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 3403 35.12 7.78     36   35.43 7.41  10  50    40 -0.39    -0.12 0.13
------------------------------------------------------------ 
group: TerceroM
   vars    n  mean   sd median trimmed  mad min max range  skew kurtosis   se
X1    1 3216 35.81 7.67     36   36.11 7.41  10  50    40 -0.37    -0.15 0.14
  1. Podemos también generar un boxplot con los mismos datos y visualizar las distribuciones para cada curso.
boxplot(poli$puntaje_au ~ poli$curso, 
        xlab = 'Curso', 
        ylab = 'Distribuación de Autoestima')

  1. Con los datos anteriores entonces podemos ilustrar el uso de ANOVA, a través de esta veremos si existen diferencias estadísticamente significativas entre uno o más de estos grupos.

A diferencia de la prueba t, donde solo eran dos grupos que comparamos, acá evaluamos de manera simultánea la hipótesis nula de que no hay diferencias entre ninguno de los grupos.

La hipótesis nula de esta prueba es:

\[ H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 = \mu_5 \]

\[ H_a: \text{No todos los promedios poblacionales son iguales} \]

# Se genera un objeto, y se utiliza aov(y~x, data=basededatos)
anova1 = aov(puntaje_au~curso, data=poli)

# Ahora para saber el reporte del anova usamos summary del objeto (anova1)
summary(anova1)
               Df  Sum Sq Mean Sq F value Pr(>F)    
curso           4    5070  1267.4   21.72 <2e-16 ***
Residuals   17936 1046680    58.4                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1743 observations deleted due to missingness
# Recuerden que los números muy pequeños son reportados en R con notación científica
format(2e-16, scientific = FALSE)
[1] "0.0000000000000002"
# Cómo calcular directamente el valor p en base al F y los grados de libertad
1 - pf(106.3,df1 = 4, df2 = 44)
[1] 0

Al igual que en la prueba t o en la prueba z, contrastaremos el p-value (pr(>F)) con nuestro alfa determinado por el grado de confianza de nuestra prueba estadística.

  • Si el p-value es mayor o igual que nuestro alfa entonces no rechazamos la hipótesis nula.

  • Si el valor p es menor que nuestro alfa entonces rechazamos la hipótesis nula.

Si queríamos realizar nuestra prueba al 99% de confianza, vale decir, con un alfa de 0.01, entonces el resultado de esta prueba ANOVA muestra que el valor p es menor que nuestro alfa, por lo que rechazaríamos la hipótesis nula.

El rechazo de la hipótesis nula significa que tenemos evidencia de que no todas las carreras se originan en una misma distribución, ya que ANOVA es solo un análisis global. Por ende, no nos dice entre cuáles grupos hay diferencias.

Tamaño del efecto

También podemos calcular el tamaño del efecto utilizando el mismo ejemplo anterior.

#Calculamos el tamaño del efecto dividiendo la suma de cuadrados intergrupal por la suma de cuadrados total (inter + intra)
Tefecto = 5070 / (5070 + 1046680)
Tefecto
[1] 0.004820537

El valor del tamaño del efecto en este caso corresponde a 0.0048 aproximadamente, lo cual nos indica que un 0.48% de la varianza de las personas en la variable autoestima puede ser explicada por la pertenencia a las distintos cursos.

Pruebas Post Hoc

Si es que queremos conocer entre qué grupos existencias, debemos hacer uso de análisis adicionales conocidos como pruebas Post Hoc. Hay muchas pruebas distintas que pueden ser usadas para identificar qué diferencias específicas existen entre grupos. En este curso vamos a usar la prueba de Tukey.

TukeyHSD(anova1, conf.level = .99)
  Tukey multiple comparisons of means
    99% family-wise confidence level

Fit: aov(formula = puntaje_au ~ curso, data = poli)

$curso
                        diff         lwr       upr     p adj
Octavo-Septimo    -0.1250093 -0.68300661 0.4329880 0.9497578
PrimeroM-Septimo  -0.4245700 -1.00704999 0.1579099 0.1228423
SegundoM-Septimo   0.4676682 -0.11365628 1.0489926 0.0670160
TerceroM-Septimo   1.1547282  0.56438602 1.7450704 0.0000000
PrimeroM-Octavo   -0.2995607 -0.88109789 0.2819764 0.4483951
SegundoM-Octavo    0.5926774  0.01229769 1.1730572 0.0079137
TerceroM-Octavo    1.2797375  0.69032553 1.8691494 0.0000000
SegundoM-PrimeroM  0.8922382  0.28828245 1.4961939 0.0000151
TerceroM-PrimeroM  1.5792982  0.96665777 2.1919387 0.0000000
TerceroM-SegundoM  0.6870600  0.07551811 1.2986020 0.0023760

Acá entonces podemos visualizar los contrastes pareados entre todos los grupos. Podemos identificar aquellos contrastes que son significativos a través de un p-value < 0.05 (p adj). En este ejemplo vemos que hay diferencias estadísticamente significativas entre los siguientes grupos:

  • TerceroM - Septimo

  • SegundoM - Octavo

  • TerceroM - Octavo

  • SegundoM - PrimeroM

  • TerceroM - PrimeroM

  • TerceroM - SegundoM