27 de octubre de 2014

Análisis de la varianza

Es una técnica para analizar la relación entre una variable dependiente métrica (cuantitativa, medida en escala de intervalo) y una variable independiente categórica llamada “Factor” (cualitativa, nominal u ordinal)

  • Variable dependiente: Métrica (intervalo o razón)
  • Variable independiente (Factor): Categórica (nominal u ordinal)

Analiza la relación entre estas variables identificando si existen diferencias significativas entre las medias de la variable dependiente en diferentes niveles del factor independiente.

Ejemplo: Rendimiento académico y nivel educativo de los padres

¿El rendimiento educativo de un estudiante está asociado de alguna forma con el nivel educativo de sus padres?

  • Variable dependiente: Rendimiento en pruebas de razonamiento matemático
  • Variable independiente (Factor): Nivel educativo del padre del alumno

Los datos provienen de las pruebas de rendimiento de alumnos de 4to de secundaria realizadas el 2001 a descargar de la siguiente dirección:

https://sites.google.com/a/pucp.pe/data_est/archivos/educ4sec.rda?attredirects=0&d=1

El libro de códigos del archivo puede verse en el siguiente archivo:

https://sites.google.com/a/pucp.pe/data_est/archivos/DIC_educ4sec.doc?attredirects=0&d=1

load("educ4sec.rda")
educ <- educ4sec

Estadísticos descriptivos por grupo

Veamos los estadísticos descriptivos de la variable dependiente según el grupo:

library(descr)
compmeans(educ$r.mat, educ$educ.pa, plot = FALSE)
## Mean value of "educ$r.mat" according to "educ$educ.pa"
##                       Mean   N Std. Dev.
## Sec. inc. o menos 697.7469 284  41.02232
## Secundaria        708.0465 205  47.50560
## Superior          724.3295 301  52.22912
## Total             710.5479 790  48.55295

Gráfico de medias

Podemos graficar las medias y sus respectivos intervalos de confianza:

library(Rmisc)
tabla.des <- summarySE(educ, measurevar="r.mat", groupvars="educ.pa", na.rm=T)
tabla.des
##             educ.pa   N    r.mat       sd       se       ci
## 1 Sec. inc. o menos 284 697.7469 41.02232 2.434227 4.791489
## 2        Secundaria 205 708.0465 47.50560 3.317935 6.541843
## 3          Superior 301 724.3295 52.22912 3.010437 5.924247
library(ggplot2)
graf.m <- ggplot(tabla.des, aes(x=educ.pa, y=r.mat)) + geom_point() + ylim(650, 750) + 
  geom_errorbar(aes(ymin=r.mat-ci, ymax=r.mat+ci), width=0.2) + 
  xlab("Nivel educativo del padre") + ylab("Puntaje en prueba de matemática") + 
  ggtitle("Intervalo de confianza al 95% de la media del puntaje\n en la prueba de razonamiento matemático,\n según nivel educativo del padre del alumno")

graf.m

Planteamiento de la prueba de hipótesis:

Paso 1: Formular la hipótesis cero y la hipótesis uno

H0 : \(\bar{y}_1 = \bar{y}_2 = \bar{y}_3 = \bar{y}_g\)

H1: \(\bar{y}_1 \neq \bar{y}_2 \neq \bar{y}_3 \neq \bar{y}_g\)

Se compara la media de cada grupo con la media global (g). Las diferencias entre la media de cada grupo y la media global se conoce como efectos principales o efectos de la prueba (EP)

\[EP_{y_i}= \bar{y}_i-\bar{y}_g\]

##                   m.grupo m.global  Ef.Pr
## Sec. inc. o menos  697.75    710.5 -12.75
## Secundaria         708.05    710.5  -2.45
## Superior           724.33    710.5  13.83

Puntuación de desviación y componentes de la varianza

Tomemos el caso del alumno 185, quien obtuvo un puntaje de 765.47 en la prueba de matemáticas y su padre tenía educación superior. En este caso:

  • Su puntuación de desviación respecto de la media global es:

\[D_{y_{185}}= y_{185} - \bar{y}_g\]

\[D_{y_{185}}= 765.47 - 710.5 = 54.97\]

  • Su puntuación de desviación respecto de la media de su grupo es:

\[D_{y_{185}}= y_{185} - \bar{y}_{sup}\]

\[D_{y_{185}}= 765.47 - 724.33 = 41.14\]

Desviación del alumno 185

Desviación del alumno 185

En el caso del alumno 185 podemos decir que su puntuación es igual a:

\[y_{185}= \bar{y}_g + (y_{185} - \bar{y}_g)\]

\[y_{185}= \bar{y}_g + (\bar{y}_{sup} - \bar{y}_g) + (y_{185} - \bar{y}_{sup})\]

\[765.47 = 710.5 + 13.83 + 41.14\]

El puntaje del alumno 185 es igual a la media global + la desviación explicada entre grupos + la desviación no explicada dentro del grupo

Modelo lineal general o de efectos aditivos

La mejor predicción del puntaje de una variable Y será \(\bar{Y}\) más los efectos de una variable independiente X

\[Y = \bar{Y} + bX + e\]

El modelo lineal general descompone cada puntuación de Y en tres partes:

  • La cantidad de Y explicada por la media global \(\bar{Y}\)
  • La cantidad de su desviación explicada por X (el efecto principal)
  • La cantidad de su desviación no explicada por X, es decir el error.

ANOVA y Suma de Cuadrados

ANOVA busca determinar cuánto de Y puede ser explicado por X. Hasta qué punto la varianza explicada por X es mayor a la varianza no explicada por X. Para ello descompone la varianza total de Y en dos partes, aplicando lo que se conoce como el método de la Suma de Cuadrados:

  • Suma total de cuadrados o varianza total:

\[SS_y=\sum_{i=1}^{n}(y_i-\bar{y})^2\]

  • Suma de cuadrados "explicada" o varianza entre grupos:

\[SS_x=\sum_{j=1}^{c}n(\bar{y}_j-\bar{y})^2\]

  • Suma de cuadrados no explicada o varianza dentro de los grupos:

\[SS_{error}=\sum_{j=1}^{c}\sum_{i=1}^{n}(y_{ij}-\bar{y}_j)^2\]

Por tanto:

\[SS_y=SS_x+SS_{error}\]

\[\sum_{i=1}^{n}(y_i-\bar{y})^2=\sum_{j=1}^{c}n(\bar{y}_j-\bar{y})^2+\sum_{j=1}^{c}\sum_{i=1}^{n}(y_{ij}-\bar{y}_j)^2\]

Cuadrados medios y estadística de F

  • Cuadrado medio entre grupos:

\[CM_E= \frac{\sum_{j=1}^{c}n(\bar{y}_j-\bar{y})^2}{K-1}\]

Donde k = al número de grupos; y k-1 son los grados de libertad de la varianza entre grupos

  • Cuadrado medio dentro de los grupos:

\[CM_D = \frac{\sum_{j=1}^{c}\sum_{i=1}^{n}(y_{ij}-\bar{y}_j)^2}{n-K}\]

donde n = número total de casos; y n-k son los grados de libertad de la varianza dentro de los grupos.

  • Estadística de F: razón entre la varianza explicada y la varianza no explicada

\[F=\frac{CM_E}{CM_D}=\frac{{\sum_{j=1}^{c}n(\bar{y}_j-\bar{y})^2}/{(K-1)}}{{\sum_{j=1}^{c}\sum_{i=1}^{n}(y_{ij}-\bar{y}_j)^2}/(n-1)}\]

Pasos de ANOVA

Paso 2: Seleccionar la distribución de muestro de la hipótesis cero

Se utiliza la distribución de F.

Paso 3: Seleccionar un nivel de significancia

Puede ser \(\alpha = 0.05\) ó \(\alpha = 0.01\), por ejemplo

Paso 4: Se calcula el estadístico de la prueba

Es este caso será un valor de F, que deberá compararse con un valor crítico de F dados los grados de libertad. F tiene dos tipos de grados de libertad, los del numerador (k-1) y los del denominador (n-k):

\[F=\frac{CM_E}{CM_D}=\frac{{\sum_{j=1}^{c}n(\bar{y}_j-\bar{y})^2}/{(K-1)}}{{\sum_{j=1}^{c}\sum_{i=1}^{n}(y_{ij}-\bar{y}_j)^2}/(n-1)}\]

Paso 5: Se decide si se acepta o rechaza HO

Se acepta H0 cuando:

  • el estadístico de F es menor que el valor crítico de F ó
  • cuando el p-value o significancia del estadístico de la prueba es mayor que el nivel de significancia de la prueba.

Caso contrario, se rechaza H0

ANOVA con R

Se utiliza la función "summary(aov(x~f))" para generar una tabla de ANOVA que muestra los diferentes componentes de la varianza de Y

summary(aov(educ$r.mat~educ$educ.pa))
##               Df  Sum Sq Mean Sq F value   Pr(>F)    
## educ$educ.pa   2  104990   52495   23.54 1.18e-10 ***
## Residuals    787 1754989    2230                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso los resultados de la prueba nos lleva a rechazar H0 y afirmar que sí existen diferencias estadísticamente significativas en los promedios de las pruenas de matemáticas entre los grupos analizados. El efecto principal de X en Y es estadísticamente significativo.

Nótese que la varianza total es:

\(S^2=(SS_x+SS_{error})/(n-1) = (104990+1754989) / 787 = 2363.38\)

Prueba de diferencias entre grupos específicos

ANOVA sólo nos dice que el efecto principal es estadísticamente significativo (las medias de los grupos no son iguales). Pero no nos dice entre qué grupos se encuentran las diferencias.

Para identificar las diferencias específicas entre grupos se utiliza el método de las diferencias altamente significativas inventado por John Tukey (TukeyHSD):

m <- aov(educ$r.mat~educ$educ.pa)
TukeyHSD(m)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = educ$r.mat ~ educ$educ.pa)
## 
## $`educ$educ.pa`
##                                  diff        lwr      upr     p adj
## Secundaria-Sec. inc. o menos 10.29964  0.1372835 20.46199 0.0461402
## Superior-Sec. inc. o menos   26.58263 17.4096332 35.75563 0.0000000
## Superior-Secundaria          16.28300  6.2416745 26.32432 0.0004435

Graficar las diferencias significativas

plot(TukeyHSD(m))

Un mejor gráfico en ggplot

p <- TukeyHSD(m)
p1 <- p[[1]]
p1 <- as.data.frame(p1)
p1$nivel <- row.names(p1)
graf.dif <- ggplot(p1, aes(x=nivel, y=diff)) + geom_point() + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=0.2)  +
  ylim(-10, 40) + geom_hline(yintercept=0, col="red", linetype = "longdash") + coord_flip() 

graf.dif

Otro ejemplo: Rendimiento según nivel educativo de la madre

compmeans(educ$r.mat, educ$educ.mam, plot = FALSE)
## Mean value of "educ$r.mat" according to "educ$educ.mam"
##                       Mean   N Std. Dev.
## Sec. inc. o menos 698.6035 379  38.35287
## Secundaria        712.4066 186  52.87332
## Superior          729.1312 225  54.04877
## Total             710.5479 790  48.55295

Gráfico de medias

Prueba de ANOVA

m2 <- aov(educ$r.mat~educ$educ.mam)
summary(m2)
##                Df  Sum Sq Mean Sq F value   Pr(>F)    
## educ$educ.mam   2  132415   66208   30.16 2.39e-13 ***
## Residuals     787 1727565    2195                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de HSD de Tukey

TukeyHSD(m2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = educ$r.mat ~ educ$educ.mam)
## 
## $`educ$educ.mam`
##                                  diff       lwr      upr     p adj
## Secundaria-Sec. inc. o menos 13.80308  3.953808 23.65234 0.0029944
## Superior-Sec. inc. o menos   30.52767 21.268687 39.78666 0.0000000
## Superior-Secundaria          16.72460  5.822024 27.62717 0.0009773

Gráfico de diferencias significativas de Tukey

Ejemplo: Horas dedicadas a las tareas según nivel educativo de la madre

compmeans(educ$h.tareas.mat, educ$educ.mam, plot = FALSE)
## Mean value of "educ$h.tareas.mat" according to "educ$educ.mam"
##                       Mean   N Std. Dev.
## Sec. inc. o menos 11.03694 379  4.727350
## Secundaria        10.77957 186  4.938521
## Superior          11.10222 225  4.922003
## Total             10.99494 790  4.828882

summary(aov(educ$h.tareas.mat~educ$educ.mam))
##                Df Sum Sq Mean Sq F value Pr(>F)
## educ$educ.mam   2     12   5.943   0.254  0.775
## Residuals     787  18386  23.362