Apuntes de R

Apunte 09: ANOVA

Author

David Torres Irribarra, Isidora Naranjo López y Tomás Vargas Silva

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:

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 24.65 7.4     25   24.87 7.41   1  40    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 25.53 7.66     26   25.78 7.41   1  41    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 24.23 7.71     25   24.46 7.41   1  40    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 26.12 7.78     27   26.43 7.41   1  41    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 26.81 7.67     27   27.11 7.41   1  41    40 -0.37    -0.15 0.14

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')

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

Tipos de ANOVA

Hasta ahora hemos revisado la lógica básica detrás de la prueba ANOVA y la hemos aplicado a un ejemplo.

Sin embargo, existen varios tipos de pruebas ANOVA, como el One-Way ANOVA, el Two-Way ANOVA, el ANOVA factorial, el ANCOVA, y otros tipos más complejos que no revisaremos en este curso

Veamos uno por uno en qué se diferencian:

One-Way ANOVA

Este tipo de ANOVA, también conocido como ANOVA unifactorial o de un solo factor, compara 2 o más grupos en una variable de respuesta en base a una única variable predictora.

Si te fijas, este ha sido el caso que hemos revisado hasta ahora en este apunte.

ANOVA Factorial y Two-Way ANOVA

El ANOVA factorial estudia las diferencias entre 2 o más grupos en una variable de respuesta en base a 2 o más variables predictoras, y puede incluir su interacción.

El Two-Way ANOVA, o ANOVA de dos factores, es un caso de ANOVA Factorial en que existen 2 variables predictoras y su posible interacción.

Aplicándolo al ejemplo

En R, cambia un poco el código para realizar un ANOVA Factorial. Revisémoslo tomando el ejemplo anterior, pero esta vez intentaremos ver si hay diferencias entre los puntajes de autoestima en base al curso y a la comuna donde viven los estudiantes.

Primero que nada, revisemos la variable comu que agregaremos al análisis.

table(poli$comu)


  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
193  51  91  31  26 194 427  20 541  71  14  26 404  67  49 146  25 180  25  32 
 21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
121  27  60  87  62  28  26  59  19 168  27  29  18  28  29 195  30  57  92 126 
 41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
 50  35 227 228 111  61  58  62  61 327  45 135  27  36 142  34  34  29  20  51 
 61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
 28  19  32  57  60  61  86 176  67 251  30 289  31 145 101  49 267  53  32  33 
 81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
 25  26 145  72  27  86 179  30  35  66  21  68  31  41 254  16  30 124  36  13 
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
471  35  31  28 125  29  42  26  60  32  60  33  61  27 194 333 217  56  65  55 
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
 60 104  89  28 141  55  26  32  34 128  29 171  25  94 558 474  29  28 521  27 
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
 90  33  73 272  36 116  61 434 223  65  50  31  26  28  32  49  88 354  26  31 
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
 83 140  63  92 145  57  21 143  57  99  60 713  30  91 307 146  61 403  72  33 
181 182 183 184 185 186 187 188 189 190 191 192 193 
 87  66 243 115 189  17  59 133  95 160  26  34  27

Como podemos ver en la tabla, hay muchas comunas en este estudio, así que, para el propósito educativo de este apunte, vamos a reducir la variable a solo 5 comunas.

#Para propósitos educativos, tomaremos las 5 comunas con más diferencias en autoestima, que en estos datos corresponden a las comunas 104, 186, 37, 5, y 91
categorias_interes <- c(104, 186, 37, 5, 91)
poli$comu_reducida <- ifelse(
  poli$comu %in% categorias_interes, #Condición = estar en las categorías de interés
  as.character(poli$comu), #Si está, mantener el valor
  NA_character_) #Si no está, poner NA


#Podemos revisar la variable y sus ptjes en autoestima graficando boxplots

boxplot(poli$puntaje_au ~ poli$comu_reducida, 
        xlab = 'Comuna', 
        ylab = 'Distribución de Autoestima')

Ahora realicemos el ANOVA factorial con nuestra variable reducida

#Agregamos el segundo factor mediante un signo de multiplicación
#aov(dependiente ~ factor1 * factor2, data = datos)

anova2 <- aov(puntaje_au ~ curso  * comu_reducida, data = poli)
summary(anova2)

                    Df Sum Sq Mean Sq F value   Pr(>F)    
curso                4  309.4    77.3   2.136  0.08431 .  
comu_reducida        4 1300.1   325.0   8.976 5.28e-06 ***
curso:comu_reducida  9 1112.3   123.6   3.413  0.00141 ** 
Residuals           77 2788.1    36.2                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
19589 observations deleted due to missingness

En este caso, manteniendo nuestro 99% de confianza (alfa = .01), encontramos que:

No hay diferencias significativas en la autoestima según el curso (p > .01)
Sí hay diferencias significativas según la comuna (p < .01)
La relación entre curso y autoestima SÍ varía significativamente según la comuna (p < .01)

¿Por qué una variable que antes era significativa, ya no lo es?

Es normal que una variable que era significativa en un One-Way ANOVA deje de serlo cuando se agregan más variables en un ANOVA factorial.

Ello puede ocurrir por varias razones:

Puede ser que la nueva variable explique parte de la diferencia que antes se atribuía a la primera variable
Puede que la interacción entre las variables sea importante y haga que el efecto de la primera variable no sea tan claro
Puede que sean efectos pequeños, y al complejizar el modelo hayan dejado de ser estadísticamente significativos.

Apliquémoslo a nuestro ejemplo:

En un primer momento, encontramos que había diferencias significativas en el puntaje de autoestima de los estudiantes según el curso en el que van.

Sin embargo, al incluir la comuna de la que provienen, este efecto se vuelve más pequeño (p = 0.08), de forma que no es estadísticamente significativo según nuestro nivel de confianza (alfa = 0.01).

Ahora bien, encontramos que sí hay diferencias significativas según la comuna, y también encontramos que la relación que había entre el autoestima y el curso se ve afectado significativamente por la interacción con la comuna de procedencia.

Es decir, al complejizar el modelo encontramos que la relación inicial entre el curso y los puntajes de autoestima no es tan fuerte como creíamos y se ve afectada por la comuna, la cual también está asociada signficativamente con dichos puntajes.

Prueba Post-Hoc de Tukey

Nuevamente, el ANOVA nos indica que sí hay diferencias entre los grupos, pero no entre cuáles grupos están las diferencias.

Con la prueba Post-Hoc de Tukey podemos analizar eso:

#Creamos la prueba de Tukey
Tukey_anova2 <- TukeyHSD(anova2)

##Como son muchísimas comparaciones, filtraremos aquellas significativas

### Filtrar comparaciones significativas para 'curso'
sig_curso <- subset(as.data.frame(Tukey_anova2$curso), `p adj` < 0.01)

### Filtrar comparaciones significativas para 'comu_reducida'
sig_comu <- subset(as.data.frame(Tukey_anova2$comu_reducida), `p adj` < 0.01)

### Filtrar comparaciones significativas para la interacción
sig_interaccion <- subset(as.data.frame(Tukey_anova2$`curso:comu_reducida`), `p adj` < 0.01)

En concordancia con el análisis anterior, vemos que:

No arroja resultados significativos según el curso

sig_curso

[1] diff  lwr   upr   p adj
<0 rows> (or 0-length row.names)

Arroja que hay diferencias significativas entre las comunas 5 y 104, entre 5 y 37, y entre 91 y 5

sig_comu

           diff        lwr       upr        p adj
5-104  9.193535   3.903660 14.483410 5.956527e-05
5-37   6.975014   1.931104 12.018924 2.106753e-03
91-5  -9.505999 -15.190954 -3.821045 1.195083e-04

Arroja que hay interacciones significativas entre algunos grupos particulares

sig_interaccion

                             diff        lwr       upr        p adj
TerceroM:37-Septimo:104  15.33333   3.825568 26.841099 7.079807e-04
Septimo:5-Septimo:104    17.90476   6.396997 29.412527 2.484892e-05
SegundoM:91-Septimo:5   -16.57143 -30.884031 -2.258826 7.633014e-03

ANCOVA

El ANCOVA, o análisis de covarianza, es una extensión del ANOVA para poder examinar diferencias controlando por una o más variables continuas (covariable(s)).

Es decir, busca estimar si hay diferencias significativas eliminando el efecto de la(s) covariable(s) para ver mejor el efecto del factor principal

Aplicándolo al ejemplo

En R también cambia el código ligeramente. Intentemos aplicarlo a nuestro ejemplo.

En este caso, podríamos querer saber si hay diferencias entre los puntajes de autoestima según la comuna, pero controlando por la edad

#La sintaxis es la siguiente:
#aov(dependiente ~ factor + covariable, data = datos)

#Aplicado a nuestro ejemplo
ancova <- aov(puntaje_au ~ comu_reducida + edad, data = poli)
summary(ancova)

              Df Sum Sq Mean Sq F value   Pr(>F)    
comu_reducida  4   1348   336.9   7.187 4.69e-05 ***
edad           1      9     8.7   0.186    0.668    
Residuals     88   4126    46.9                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
19590 observations deleted due to missingness

En este caso, encontramos que sí hay diferencias significativas en el puntaje de autoestima entre las distintas comunas (de las 5 seleccionadas; p < .01), aún cuando controlamos por la edad, y encontramos que esta variable no afecta la relación (p > .01)

Prueba Post-Hoc de Tukey

Nuevamente, podemos ver entre qué grupos específicamente existe la diferencia con una prueba post-hoc.

Sin embargo, la función TukeyHSD() no considera las covariables, por lo que es necesario usar otro método. Por ejemplo, existen las pruebas llamadas Estimated Marginal Means, que comparan las medias entre grupos específicos controlando la covariable, y que integran la corrección por múltiples comparaciones propia de la prueba de Tukey.

# Comparaciones post-hoc ajustadas con Tukey
library(emmeans)

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

emmeans(ancova, pairwise ~ comu_reducida, adjust = "tukey")

$emmeans
 comu_reducida emmean   SE df lower.CL upper.CL
 104             31.0 1.57 88     27.9     34.2
 186             39.5 2.42 88     34.7     44.3
 37              34.6 1.41 88     31.8     37.4
 5               40.6 1.62 88     37.3     43.8
 91              31.1 1.66 88     27.8     34.4

Confidence level used: 0.95 

$contrasts
 contrast                            estimate   SE df t.ratio p.value
 comu_reducida104 - comu_reducida186   -8.510 2.87 88  -2.964  0.0312
 comu_reducida104 - comu_reducida37    -3.538 2.28 88  -1.553  0.5312
 comu_reducida104 - comu_reducida5     -9.548 2.22 88  -4.310  0.0004
 comu_reducida104 - comu_reducida91    -0.119 2.27 88  -0.052  1.0000
 comu_reducida186 - comu_reducida37     4.973 2.82 88   1.762  0.4021
 comu_reducida186 - comu_reducida5     -1.038 2.91 88  -0.357  0.9965
 comu_reducida186 - comu_reducida91     8.391 2.94 88   2.858  0.0413
 comu_reducida37 - comu_reducida5      -6.011 2.20 88  -2.733  0.0571
 comu_reducida37 - comu_reducida91      3.419 2.20 88   1.555  0.5297
 comu_reducida5 - comu_reducida91       9.429 2.32 88   4.067  0.0010

P value adjustment: tukey method for comparing a family of 5 estimates

En este caso, encontramos que hay diferencias significativas controlando por la edad entre las siguientes comunas, con un 99% de confianza:

Entre las comunas 104 y 5 (p < .01)
Entre las comunas 5 y 91 (p < .01)

BONUS: ¿Por qué el orden importa en ANOVA?

Las pruebas ANOVA siguen una lógica secuencial. Es decir, cada variable intenta explicar la variabilidad NO explicada por la(s) anterior(es).

De hecho, es por eso que la covariable en ANCOVA va al final, pues el propósito es controlar por esa variable.

La regresión, por otro lado, sigue una lógica no secuencial, que busca ver cuánto estiman todas las variables predictoras juntas.

Ahora bien, se puede hacer una regresión con lógica secuencial si se realiza un modelo anidado, el cual busca explicar la variabilidad por capas (contenido que va más allá del alcance de este curso).

De la misma forma, se puede hacer que un ANOVA siga una lógica no-secuencial, utilizando otros tipos de sumas de cuadrados (por defecto se usa Tipo I, que es secuencial, pero existen los Tipo II, no secuencial sin interacción, y Tipo III, no secuencial con interacción).