ANOVA
Alfredo Cárdenas-Rivera
20/10/2020
CLASE 11
ANÁLISIS DE VARIANZA
1. OBJETIVO
El objetivo de este capítulo es evaluar la igualdad de medias entre 3 o más grupos independientes
================================================================================
2. PUNTOS CLAVE
- La prueba de ANOVA nos permite evaluar la diferencia estadística entre medias de tres a más grupos.
| Parámetro | Muestra | Población |
|---|---|---|
| Proporción | \(\hat{p}\) | \(p\) |
| Media | \(\bar{x}\) | \(\mu\) |
| Desviación estandar | \(s\) | \(\sigma\) |
| Varianza | \(s^2\) | \(\sigma^2\) |
| \(H_0\) es verdadera | \(H_0\) es falsa | |
|---|---|---|
| Rechaza \(H_0\) | ERROR TIPO I (\(\alpha\)) | Decisión correcta |
| Acepta \(H_0\) | Decisión correcta | ERROR TIPO II (\(\beta\)) |
================================================================================
3. TEMARIO
- ANOVA de un factor
- ANOVA de dos factores
================================================================================
4. ANOVA DE UN FACTOR
El método de análisis de varianza de un factor permite evaluar la hipótesis nula donde se cumple que las medias de tres o más poblaciones independientes son iguales. De manera formal lo podemos expresar de la siguiente manera:
\[H_0:\mu_1=\mu_2=\mu_3\]
4.1 Distribución F
Se usa la distribución de F para aproximar los análisis de ANOVA. La distribución de F tiene las siguientes caracterísiticas:
- La distribución F es sesgada a la derecha, por lo tanto NO ES SIMÉTRICA como la distribución Z o la t-student.
- Los valores de F no pueden tomar valores negativos
- La distribución de F varía con los grados de libertad de la muestra.
Para demostrar la distribución F, usaremos R para simular algunas gráficas. Para ello vamos a utilizar algunos paquetes como:
- dplyr
- ggplot2
- tidyr
Asi que primero cargaremos las librerías a nuestro entorno de trabajo.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyr)Ahora que ya tenemos las librerías cargadas, generaremos los valores de x, usando una secuencia de 0 a 7 (recordemos que los valores de la distribución de F no puede ser negativa), usando la función seq(), con intervalos de 0.001
head(data.frame(chisq = seq(from = 0, to = 7, by = 0.01)))## chisq
## 1 0.00
## 2 0.01
## 3 0.02
## 4 0.03
## 5 0.04
## 6 0.05Ahora que tenemos los valores del eje x, generaremos los valores de y (que es la densidad) con la función dchisq(). Esta función tiene dos parámetros que debemos llenar:
- x= donde se define el valor de x antes generado
- df= donde se definen los grados de libertad
Una función util que vamos a utilizar es %>% que está definida en el paquete “magrittr” y que es incorporada en el paquete “dplyr”. Esta función lo que hace es pasar el argumento del lado izquierdo a la función definida en lado derecho, por ejemplo
iris%>%head()## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaOtro ejemplo sencillo
c(5,6,8,7,9)%>%mean()## [1] 7Utilizando esta función, podemos usar los valores definidos en el data.frame para generar los resultados calculados por la función dchisq(). Para crear un data.frame e ir agregando uno a una las columnas que estamos generando utilizaremos la función mutate() que está definida en el paquete “dplyr”. De esta forma, al data.frame generado por la función data.frame() en la cual solo hay una columna con los valores de 0 a 70 separados por 0.001, se iran ageregando las columnas definidas en la función mutate(). Cada columna que se genera con la función dchisq() usa los valores previos del data.frame, para llamar estos valores utilizamos la función %>%, de la siguiente forma
data.frame(chisq = seq(from = 0, to = 70, by = 0.001))%>%
mutate(df_05 = dchisq(x = chisq, df = 5),
df_15 = dchisq(x = chisq, df = 15),
df_30 = dchisq(x = chisq, df = 30))%>%head()## chisq df_05 df_15 df_30
## 1 0.000 0.000000e+00 0.000000e+00 0.000000e+00
## 2 0.001 4.203119e-06 9.330933e-26 3.498844e-58
## 3 0.002 1.188227e-05 8.441175e-24 5.729640e-54
## 4 0.003 2.181822e-05 1.177007e-22 1.671814e-51
## 5 0.004 3.357455e-05 7.632443e-22 9.378060e-50
## 6 0.005 4.689841e-05 3.253580e-21 2.131258e-48Ahora lo que necesitamos es reordenar esta base de datos, para que solo tenga 3 columnas:
- chisq que sea el valor del eje x de la distribución F
- df que contenga el nombre del factor donde se define los grado de libertad
- density que es el resultado del cálculo de la densidad para cada valor de x definido en la columna chisq
Para reordenar las columnas en fila utilizaremos otra función gather() del paquete “tidyr”. En esta función definimos que columna se debe conservar, en este caso la columna sería chisq, para lo cual la definimos como “-chisq” en nuestra función. Usamos el parámetro “key=” para definir el nombre de la columna donde serán contenidos nombre del grupo de datos (factores) y el parámetro “value=” permite definir el nombre de la columna que contendrá los valores.
data.frame(chisq = seq(from = 0, to = 70, by = 0.001)) %>%
mutate(df_05 = dchisq(x = chisq, df = 5),
df_15 = dchisq(x = chisq, df = 15),
df_30 = dchisq(x = chisq, df = 30)) %>%
gather(key = "grados", value = "valor", -chisq) %>%
head()## chisq grados valor
## 1 0.000 df_05 0.000000e+00
## 2 0.001 df_05 4.203119e-06
## 3 0.002 df_05 1.188227e-05
## 4 0.003 df_05 2.181822e-05
## 5 0.004 df_05 3.357455e-05
## 6 0.005 df_05 4.689841e-05Con esta base de datos, tipo data frame podemos generar nuestro gráfico con el paquete “ggplot2”, pero antes verifiquemos la estructura de nuestra base de datos
data.frame(chisq = seq(from = 0, to = 70, by = 0.001)) %>%
mutate(df_05 = dchisq(x = chisq, df = 5),
df_15 = dchisq(x = chisq, df = 15),
df_30 = dchisq(x = chisq, df = 30)) %>%
gather(key = "grados", value = "valor", -chisq) %>%
str()## 'data.frame': 210003 obs. of 3 variables:
## $ chisq : num 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 ...
## $ grados: chr "df_05" "df_05" "df_05" "df_05" ...
## $ valor : num 0.00 4.20e-06 1.19e-05 2.18e-05 3.36e-05 ...Utilizaremos la gráfica de lineas del paquete “ggplot2” que está definida en la función geom_line. En esta función debemos definir el parámetro aes(), donde indicamos asignamos los valores de las variables:
- x= donde se define los valores de x
- y= donde se define los valores de y
- color= aquí indicamos la variable factor que va a agrupar nuestros datos y les va a asignar un color distinto.
data.frame(chisq = seq(from = 0, to = 70, by = 0.001)) %>%
mutate(df_05 = dchisq(x = chisq, df = 5),
df_15 = dchisq(x = chisq, df = 15),
df_30 = dchisq(x = chisq, df = 30)) %>%
gather(key = "grados", value = "valor", -chisq) %>%
ggplot() +
geom_line(aes(x = chisq,
y = valor,
color = grados))
Lo que observamos es que al aumentar la \(n\) de la muestra, la gráfica se va acercando a distribución Z o t-student. Para terminar está gráfica lo que podemos hacer es cambiar título y el nombre de los ejes con la función labs() definida en el paquete “ggplot2”
data.frame(chisq = seq(from = 0, to = 70, by = 0.001)) %>%
mutate(df_05 = dchisq(x = chisq, df = 5),
df_15 = dchisq(x = chisq, df = 15),
df_30 = dchisq(x = chisq, df = 30)) %>%
gather(key = "grados", value = "valor", -chisq) %>%
ggplot() +
geom_line(aes(x = chisq,
y = valor,
color = grados)) +
labs(title = "Distribución F",
x = "Valores de Chi cuadrado",
y = "Densidad")
4.2 Requisitos
- La distribución de las diferentes poblaciones tienen que aproximarse a una distribución normal.
- La varianza de las muestras difieren poco entre ellas (si el tamaño muestral es el mismo, la varianza puede variar hasta 9 veces sin afectar el resultado de ANOVA)
- Los datos son cuantitativos aleatorios independientes.
- Los datos se encuentran categorizados de una sola manera.
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa4.3 Analisis de ANOVA en R
Para el ejemplo trabajaremos con la base de datos de iris, específicamente con la variable “ancho de pétalos”. Para ello, generaremos un data frame con los datos de especie y ancho de pétalo
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosadb <- iris[, c("Petal.Width","Species")]
db## Petal.Width Species
## 1 0.2 setosa
## 2 0.2 setosa
## 3 0.2 setosa
## 4 0.2 setosa
## 5 0.2 setosa
## 6 0.4 setosa
## 7 0.3 setosa
## 8 0.2 setosa
## 9 0.2 setosa
## 10 0.1 setosa
## 11 0.2 setosa
## 12 0.2 setosa
## 13 0.1 setosa
## 14 0.1 setosa
## 15 0.2 setosa
## 16 0.4 setosa
## 17 0.4 setosa
## 18 0.3 setosa
## 19 0.3 setosa
## 20 0.3 setosa
## 21 0.2 setosa
## 22 0.4 setosa
## 23 0.2 setosa
## 24 0.5 setosa
## 25 0.2 setosa
## 26 0.2 setosa
## 27 0.4 setosa
## 28 0.2 setosa
## 29 0.2 setosa
## 30 0.2 setosa
## 31 0.2 setosa
## 32 0.4 setosa
## 33 0.1 setosa
## 34 0.2 setosa
## 35 0.2 setosa
## 36 0.2 setosa
## 37 0.2 setosa
## 38 0.1 setosa
## 39 0.2 setosa
## 40 0.2 setosa
## 41 0.3 setosa
## 42 0.3 setosa
## 43 0.2 setosa
## 44 0.6 setosa
## 45 0.4 setosa
## 46 0.3 setosa
## 47 0.2 setosa
## 48 0.2 setosa
## 49 0.2 setosa
## 50 0.2 setosa
## 51 1.4 versicolor
## 52 1.5 versicolor
## 53 1.5 versicolor
## 54 1.3 versicolor
## 55 1.5 versicolor
## 56 1.3 versicolor
## 57 1.6 versicolor
## 58 1.0 versicolor
## 59 1.3 versicolor
## 60 1.4 versicolor
## 61 1.0 versicolor
## 62 1.5 versicolor
## 63 1.0 versicolor
## 64 1.4 versicolor
## 65 1.3 versicolor
## 66 1.4 versicolor
## 67 1.5 versicolor
## 68 1.0 versicolor
## 69 1.5 versicolor
## 70 1.1 versicolor
## 71 1.8 versicolor
## 72 1.3 versicolor
## 73 1.5 versicolor
## 74 1.2 versicolor
## 75 1.3 versicolor
## 76 1.4 versicolor
## 77 1.4 versicolor
## 78 1.7 versicolor
## 79 1.5 versicolor
## 80 1.0 versicolor
## 81 1.1 versicolor
## 82 1.0 versicolor
## 83 1.2 versicolor
## 84 1.6 versicolor
## 85 1.5 versicolor
## 86 1.6 versicolor
## 87 1.5 versicolor
## 88 1.3 versicolor
## 89 1.3 versicolor
## 90 1.3 versicolor
## 91 1.2 versicolor
## 92 1.4 versicolor
## 93 1.2 versicolor
## 94 1.0 versicolor
## 95 1.3 versicolor
## 96 1.2 versicolor
## 97 1.3 versicolor
## 98 1.3 versicolor
## 99 1.1 versicolor
## 100 1.3 versicolor
## 101 2.5 virginica
## 102 1.9 virginica
## 103 2.1 virginica
## 104 1.8 virginica
## 105 2.2 virginica
## 106 2.1 virginica
## 107 1.7 virginica
## 108 1.8 virginica
## 109 1.8 virginica
## 110 2.5 virginica
## 111 2.0 virginica
## 112 1.9 virginica
## 113 2.1 virginica
## 114 2.0 virginica
## 115 2.4 virginica
## 116 2.3 virginica
## 117 1.8 virginica
## 118 2.2 virginica
## 119 2.3 virginica
## 120 1.5 virginica
## 121 2.3 virginica
## 122 2.0 virginica
## 123 2.0 virginica
## 124 1.8 virginica
## 125 2.1 virginica
## 126 1.8 virginica
## 127 1.8 virginica
## 128 1.8 virginica
## 129 2.1 virginica
## 130 1.6 virginica
## 131 1.9 virginica
## 132 2.0 virginica
## 133 2.2 virginica
## 134 1.5 virginica
## 135 1.4 virginica
## 136 2.3 virginica
## 137 2.4 virginica
## 138 1.8 virginica
## 139 1.8 virginica
## 140 2.1 virginica
## 141 2.4 virginica
## 142 2.3 virginica
## 143 1.9 virginica
## 144 2.3 virginica
## 145 2.5 virginica
## 146 2.3 virginica
## 147 1.9 virginica
## 148 2.0 virginica
## 149 2.3 virginica
## 150 1.8 virginicaUna vez generamos nuestra base de datos, usaremos la función ggboxplot() para generar un boxplot de nuestros datos. Esta función está definida en el paquete “ggpubr”
library(ggpubr)
ggboxplot(db,
x = "Species",
y = "Petal.Width",
color = "Species",
palette = colores1,
order = c("virginica", "versicolor", "setosa"),
ylab = "Ancho de pétalo", xlab = "Especie", title = "Ancho de petalo")
Para el análisis de ANOVA, podemos usar la función aov() definida en R, de la siguiente forma:
res <- aov(formula = Petal.Width~Species, data = db)
res## Call:
## aov(formula = Petal.Width ~ Species, data = db)
##
## Terms:
## Species Residuals
## Sum of Squares 80.41333 6.15660
## Deg. of Freedom 2 147
##
## Residual standard error: 0.20465
## Estimated effects may be unbalancedsummary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 80.41 40.21 960 <2e-16 ***
## Residuals 147 6.16 0.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para hacer una prueba multiple de comparación, podemos usar la prueba poshoc de tukey definida en la función TukeyHSD()
TukeyHSD(res)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Petal.Width ~ Species, data = db)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa 1.08 0.9830903 1.1769097 0
## virginica-setosa 1.78 1.6830903 1.8769097 0
## virginica-versicolor 0.70 0.6030903 0.7969097 0