ANOVA

Hoy vamos a estudiar una de las pruebas de hipótesis más populares: ANOVA (ANalysis Of VAriance).

Esta es una prueba que utilizamos cuando la variable dependiente es continua y la variable independiente es categórica con tres o más niveles. Recuerden que si tenemos dos categóricas utilizamos ji-cuadrada y si tenemos continua con categórica de 2 niveles utilizamos t-test.

Existen múltiples tipos de ANOVA según las comparaciones que se pueden realizar pero nosotros nos quedaremos con la versión más simple: one-way ANOVA

Seguiremos con nuestro análisis de la escala de depresión en la ENCOVID. Sólo que ahora nuestra pregunta es distinta:

¿Tenemos diferencias significativas en depresión por nivel de escolaridad? ¿Quién dirían que reporta mayor depresión alta o baja educación?

Vamos a explorar nuestros datos

library(tidyverse)

There were 22 warnings (use warnings() to see them)

library(skimr)
library(infer)
library(ggridges)
library(car)

library(readr)
base_covid_abril2020_final <- read_csv("BD/base_covid_abril2020_final.csv")

abril <-base_covid_abril2020_final


#sexo
abril <- abril %>% 
  mutate(mujer = ordered(sexo, levels = c("0", "1"), labels= c("Varones", "Mujeres")))

abril <- abril %>% 
  mutate(educ= case_when(
    p60_socio==1 ~ 1, #sin
    p60_socio==2 ~ 1, #Preescolar
    p60_socio==3 ~ 1, #Primaria
    p60_socio==4 ~ 2, #Secundaria
    p60_socio==5 ~ 3, #Preparatoria
    p60_socio==6 ~ 4, #Superior
    p60_socio==7 ~ 4, #Maestria
    p60_socio==8 ~ 4, #Doctorado 
    p60_socio==9 ~ NA_real_, #NS/NR
  ))


abril <- abril %>% 
  mutate(educ= ordered(educ, levels=c("1", "2", "3", "4"), labels= c("Primaria", "Secundaria", "Preparatoria", "Superior")))

table(abril$educ, useNA = "ifany")

    Primaria   Secundaria Preparatoria     Superior         <NA> 
         203          246          135          236           13 

# CESD-7

abril <- abril %>% 
  rename(
    d_tristeza="p40_depresion",
    d_concentracion="p41_depresion",
    d_depresion="p42_depresion",  
    d_esfuerzo="p43_depresion",
    d_dormir="p44_depresion",
    d_disfrutar="p45_depresion",
    d_triste="p46_depresion")


abril <- abril %>% 
  mutate(d_total = d_tristeza + d_concentracion + d_depresion + d_esfuerzo + d_dormir + d_disfrutar + d_triste)

table(abril$d_total)

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
153  60  56  85  42  50  53  33  26  37  21  21  20  21  20  10   8  17  17   6   5  10 

summary(abril$d_total)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  0.000   1.000   4.000   5.792   9.000  21.000      62 

abril <- abril %>% 
  mutate(dep_level= if_else(d_total >=9, 1, 0))  #Punto de corte en >=9


abril <- abril %>% 
  mutate(dep_level2= ordered(dep_level, levels=c("0", "1"), labels= c("Sin Síntomas", "Con Síntomas")))

Ahora que tenemos datos definimos:

variable dependiente: d_total (min=0; max=21)

variable independiente: educ (4 niveles)

Vamos a primero entender el problema. Veamos la distribución del puntaje de depresión con su media= 5.79

Otra forma de ver los mismos datos. Aquí en lugar de contar, R calcula la proporción de respuestas como proporción de área bajo la curva.

abril %>% 
  ggplot(aes(x= d_total))+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8) +
  geom_vline(xintercept = 5.79, color= "red")+
  theme_minimal()

Y ahora las stats por grupo.

Ahora tenemos 4 medias. ¿Qué grupo tiene mayor depresión? La pregunta, claro está, es si esta diferencia es estadísticamente significativa.

Veamos otra representación del problema, porque no sólo estamos comprando 4 medias, sino 4 distribuciones

abril %>%
  dplyr::select(d_total,educ) %>% 
  drop_na() %>% 
  ggplot(aes(x = d_total, y = educ, fill = educ)) +
  stat_density_ridges(quantile_lines = TRUE,
                      quantile_fun=function(x,...)mean(x), # mean
                      quantiles = 2, scale = 3, color = "black",
                      alpha = .7)+
  xlim(0, 22)+
  scale_x_continuous(breaks=seq(0,22,1))+
  geom_vline(xintercept = 5.79, color= "red")+
  labs(x = "Puntaje depresion", y = "nivel educativo")+
  scale_fill_brewer(palette = "Set2")+
  theme_minimal()+
  theme(legend.position = "none")

El tema con la ANOVA es que estamos haciendo dos comparaciones simultáneas. Primero vemos si hay diferencias en la varianza entre-grupos (between groups) y luego intra-grupos.

• Varianza entre grupos: tenemos que ver la distancia que hay entre la gran media (5.79) y las medias de cada grupo (7.4, 6.2, 5.5, 4.2).

• Varianza intra grupos: también vemos la distancia que hay entre la media grupal (7.4) y los puntajes de cada persona con primaria o menos.

Lo que hace ANOVA es dividir between/ within. Si la varianza between es mayor (> 1), entonces las diferencias entre grupos son significativas. Si la varianza within es mayor (< 1), entonces no podemos afirmar que hay diferencias entre grupos. Y claro que después enfrentamos el problema de definir qué tan grande debe de ser el resultado de esa división para concluir que no es un resultado atribuible al azar.

La ANOVA se caracteriza por apoyarse en la distribución F. Esta distribución permite examinar como stat la razón entre between y within.

Como en toda prueba de hipótesis, definimos la nula como punto de partida

• H0: no hay diferencias entre grupos
• HA: al menos un grupo difiere de los demás grupos

Noten que la HA no es concluyente porque basta con que alguno sea diferente, pero no sabremos cuál de ellos.

Hagamos primero la forma computacional con infer.

abril1 <- abril %>%
    select(d_total,educ) %>% 
  drop_na() %>% 
  mutate(educ1= as.character(educ))


f_observado <- abril1 %>%
  specify(d_total ~ educ1) %>%
  calculate(stat = "F")


null_dist <- abril1 %>%
  specify(d_total ~ educ1) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "F")

null_dist

NA

Y ahora comparemos con la estat observada

null_dist %>%
  visualize() + 
  shade_p_value(f_observado,
                direction = "greater")

y el valor_p

p_value <- null_dist %>%
  get_p_value(obs_stat = f_observado,
              direction = "greater")

Please be cautious in reporting a p-value of 0. This result is an approximation based on the number of `reps` chosen in the `generate()` step. See `?get_p_value()` for more information.

p_value

¿Qué interpretamos?

Para ANOVA tomemos la ruta tradicional porque queremos explorar más detalles.

Les recomiendo mucho este tutorial

Lo primero es que las pruebas de hipótesis tienen supuestos.

• Independencia: ¿hay relación entre lo que responde una persona con primaria y lo que responde con secundaria? No, no tiene nada que ver: son independientes. Diseños pre-post serían dependientes o pareadas

• Normalidad: las distribuciones de cada grupo aproximan una distribución normal. No del todo pero más o menos asemejan una curva; aunque hay muchos ceros.

• Equivalencia de varianza: al interior de cada grupo, la varianza es similar (within)

• Outliers: no tenemos observaciones extremas en ningún grupo.

Vamos a ver con otra gráfica los últimos dos supuestos.

abril %>%
  select(d_total,educ) %>% 
  drop_na() %>%
  ggplot( aes(x=educ, y=d_total, fill=educ)) +
    geom_boxplot() +
    geom_jitter(color="grey", size=0.7, alpha=0.5) +
  stat_summary(fun=mean, geom="point", shape=23, size=2, fill= "red")+
  theme_minimal()

Aquí hay muchas cosas sucediendo. Primero, interpretamos los boxplots

• ¿Qué nos dicen las cajas? límite inferior, superior y raya en medio
• Yo le añadí ese diamante rojo: ¿qué creen que es?
• ¿Qué son los puntitos grises claro - los copos de nieve?
• ¿Qué son los puntos negros arriba?
• ¿Qué nos dice esta gráfica sobre la homogeneidad de varianza?

A mi los outliers no me preocupan. Pero la varianza realmente no lo sé gráficamente.

Vemos formalmente con una prueba de Levene de homogeneidad de varianza

leveneTest(d_total ~ educ, data = abril) # noten que volví a abril normal

Levene's Test for Homogeneity of Variance (center = median)
       Df F value    Pr(>F)    
group   3  6.2477 0.0003425 ***
      756                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

El valor-p es 0.00034. Si la H0 es que hay homogeneidad de varianza, entonces este resultado indica que nuestros grupos tienen heterogenidad de varianza (HA).

Lo correcto al correr una ANOVA, dado esta violación a un supuesto, es lo siguiente:

oneway.test(d_total ~ educ,
  data = abril,
  var.equal = FALSE) # non-equal variances

    One-way analysis of means (not assuming equal variances)

data:  d_total and educ
F = 11.861, num df = 3.00, denom df = 374.53, p-value = 1.945e-07

Tenemos un valor p= 0.00000001945. Es bajo. Nos quedamos con HA.

Siento que la varianza no era tan diferente. Y la ANOVA es robusta a violaciones de supuestos.

Vamos a ver qué pasa si lo reptimos con una ANOVA normal.

oneway.test(d_total ~ educ,
  data = abril,
  var.equal = TRUE) # equal variances

    One-way analysis of means

data:  d_total and educ
F = 11.516, num df = 3, denom df = 756, p-value = 2.179e-07

Noten cómo el estadístico F es casi igual. Y el valor-p sigue siendo MUY bajo. La decisión no cambia

Como los resultados se ven muy parecidos voy a usar otro comando que hace la misma ANOVA pero también puede hacer unos trucos extra. Vean que la sintaxix es casi igual. Sólo no nos deja especificar el supuesto de varianza.

anova_educ <- aov(d_total ~ educ,
  data = abril)

anova_educ # veamos qué hay adentro. Es un objeto de tipo modelo estadístico (list)

Call:
   aov(formula = d_total ~ educ, data = abril)

Terms:
                    educ Residuals
Sum of Squares   1023.15  22389.15
Deg. of Freedom        3       756

Residual standard error: 5.441993
Estimated effects may be unbalanced
73 observations deleted due to missingness

Lo interesante de estos objetos es que tienen MUCHAS cosas dentro que podemos extraer.

summary(anova_educ)

             Df Sum Sq Mean Sq F value   Pr(>F)    
educ          3   1023   341.1   11.52 2.18e-07 ***
Residuals   756  22389    29.6                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
73 observations deleted due to missingness

Ahí esta el mismo estadístico F. Y un valor muy pequeño. Misma decisión (HA)

Vean algo padre si tienen dudas de cómo interpretar esto. Esas son las cosas que van surgiendo en R. Les debo el concepto de tamaño del efecto.

(Recuerden que la primera vez lo deben instalar. Solo una vez. Borren el gato previo a install y corranlo. Luego vuelvan a poner el # para que ignore ese comando.)

#install.packages("remotes")
#remotes::install_github("easystats/report")

library(report)

report(anova_educ)

The ANOVA (formula: d_total ~ educ) suggests that:

  - The main effect of educ is statistically significant and small (F(3, 756) = 11.52, p < .001; Eta2 = 0.04, 95% CI [0.02, 1.00])

Effect sizes were labelled following Field's (2013) recommendations.

De tarea exploren esta página; cap 6 ANOVA

Sigamos. Vamos a verificar si violamos el supuesto de normalidad. Si fuera nomal, todas se acostarían sobre la diagonal. Si se separan mucho de la diagonal, se viola.

plot(anova_educ, which = 2)

Otra prueba similar, pero con IC. Sí se sale de los intervalos…

library(car) #instalenlo

qqPlot(anova_educ$residuals,
  id = FALSE)

La prueba formal. H0 los residuales del modelo tienen una distribución normal. HA no es una distribución normal y se viola el supuesto.

shapiro.test(anova_educ$residuals)

    Shapiro-Wilk normality test

data:  anova_educ$residuals
W = 0.91915, p-value < 2.2e-16

El valor p es muuuuy pequeño. así que esos 0 de depresión parece que sí pesan.

Debemos usar una prueba no-paramétrica: kruskal wallis. Comparemos mejor medianas (misma sintaxis, misma H0)

kruskal.test(d_total ~ educ, data = abril)

    Kruskal-Wallis rank sum test

data:  d_total by educ
Kruskal-Wallis chi-squared = 33.379, df = 3, p-value = 2.679e-07

Otro valor p muy pequeño. Otra vez nos quedamos con HA

Nuestra conclusión se va fortaleciendo. Varias pruebas, misma decisión

Vean cómo a pesar de violar el supuesto la ANOVA normal nos dice bastante.

Por eso me voy a quedar con el objeto que me da más cositas.

Vamos ahora a algo aún más interesante: ¿Qué grupos difieren entre sí?

Para eso necesitamos hacer muchas pruebas simultaneas.

• Primaria vs secundaria
• Primaria vs prepratoria
• Primaria vs superior

… y así con las demás…

El problema es que al hacer muchas comparaciones es fácil encontrar falsos positivos por casualidad. Lo bueno es que esto se corrige si se hacen todas las pruebas simulataneamente.

Vamos a usar los que se conocen como pruebas post-hoc de Tukey (famoso estadístico, inventó los boxplots)

TukeyHSD(anova_educ) #volví a mi objeto anova

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

Fit: aov(formula = d_total ~ educ, data = abril)

$educ
                             diff       lwr        upr     p adj
Secundaria-Primaria     -1.179745 -2.573886  0.2143958 0.1300643
Preparatoria-Primaria   -1.832112 -3.455971 -0.2082541 0.0197337
Superior-Primaria       -3.106716 -4.503611 -1.7098199 0.0000001
Preparatoria-Secundaria -0.652367 -2.208988  0.9042539 0.7024513
Superior-Secundaria     -1.926970 -3.245102 -0.6088386 0.0010325
Superior-Preparatoria   -1.274603 -2.833692  0.2844852 0.1524115

Ahora sí tenemos resultados más interesantes…

1. Vean la columna abajo de $educ. Son las comparaciones 1:1. Empecemos por el primer renglón. ¿son diferentes los puntajes de depresión entre secundaria y primaria?

2. La segunda columna es el estadístico “diff”. Sí, es una prueba de diferencia de medias (t-test). Primaria tiene, en promedio, 1.17 más puntos en depresión que secundaria.

3. Las siguientes dos columnas son nuestros queridos intervalos de confianza. Noten cómo el estadísticos de diff tiene un intervalo de -2.57 hasta 0.21. ¿Notan algo sospechoso? Sí, el intervalo incluye el 0. Esto significa que si la muestra se hubiera repetido 100 veces, en más de 5 ocasiones no hubiera habido diferencia entre estos niveles. De hecho, en varias de ellas secundaria habría sido mayor.

4. Nuestro famoso valor-p. Está ajustado porque Tukey está cuidando el problema de falsos positivos cuando se hacen varias pruebas simultáeas. El valor-p es 0.13. Significa que hay un 13% de probabilidad de que H= ocurra. En otras palabras, si repetimos la muestra 100 veces, en 13 no encontraríamos diferencias. Esas son muchas. Más de las 5 que estamos dispuestos a tolerar. Y ese 0.13 ciertamente es mayor al .05 que definimos a priori como frontera del 95% de confianza. Así que no rechazamos H0: no hay diferencia entre secundaria y primaria.

¿En qué niveles sí hay diferencias?

Vean las combinaciones en las que el valor-p es menor a 0.05.

¿Les hace sentido teóricamente? ¿Dónde marcan la frontera? ¿Cuál es la zona gris?

Vamos a verlo de otra forma. No es la más atractiva, pero sí útil. Están en el mismo orden que arriba.

¿qué interpretamos?

plot(TukeyHSD(anova_educ))

Si el ala del avión del IC cruza el 0, entonces esa comparación no es estadísticamente significativa.

A veces queremos hacer menos comparaciones para tener mayor poder estadístico. Y usamos una categoría de referencia. Voy a ver todos contra Primaria, la del puntaje más alto. En ese caso uso las prubeas de Dunett.

library(multcomp) #intalar

dunnett <- glht(anova_educ,
  linfct = mcp(educ = "Dunnett"))

summary(dunnett)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = d_total ~ educ, data = abril)

Linear Hypotheses:
                             Estimate Std. Error t value Pr(>|t|)    
Secundaria - Primaria == 0    -1.1797     0.5415  -2.179   0.0768 .  
Preparatoria - Primaria == 0  -1.8321     0.6307  -2.905   0.0107 *  
Superior - Primaria == 0      -3.1067     0.5425  -5.726   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Ahora la gráfica. De nuevo, sigan el orden del output previo.

plot(dunnett)

Yo creo que con estos resultados ya es más claro que hay dos grupos primaria-secundaria y luego preparatoria-superior. Estas ya son conslusiones más sólidas

Y ahora una de pilón… Una joyita más o menos reciente en R.

library(ggstatsplot)

package ‘ggstatsplot’ was built under R version 4.0.2You can cite this package as:
     Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
     Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

ggbetweenstats(
  data = abril,
  x = educ,
  y = d_total,
  type = "parametric", # ANOVA or Kruskal-Wallis
  var.equal = FALSE, # ANOVA or Welch ANOVA
  plot.type = "box",
  pairwise.comparisons = TRUE,
  pairwise.display = "significant",
  centrality.plotting = FALSE,
  bf.message = FALSE
)

Ahora con kruskall-Wallis. Vean Primaria vs Preparatoria

ggbetweenstats(
  data = abril,
  x = educ,
  y = d_total,
  type = "nonparametric", # ANOVA or Kruskal-Wallis
  var.equal = FALSE, # ANOVA or Welch ANOVA
  plot.type = "box",
  pairwise.comparisons = TRUE,
  pairwise.display = "significant",
  centrality.plotting = FALSE,
  bf.message = FALSE
)

Por último, les dejo un recurso donde pueden consultar las muchas pruebas existentes según necesiten.

Recuerden que la decisión depende de sus tipos de variables y si usan paramétrica o no.

Este es un buen resumen