Estadística básica agronómica
Raúl Ortiz
21 de noviembre de 2017
Puedes seguir el tutorial por vídeo en YouTube
1 Introducción al diseño de experimentos.
2 Elementos de inferencia estadística.
3 Pruebas paramétricas y procedimientos análogos no-paramétricos.
4 Cómo comprobar la normalidad de los datos.
4.1 Intro.
Con una cantidad de datos de la muestra suficientemente grande (n>30) la violación del supuesto de la normalidad no debería causar mayores problemas. Esto implica que podemos ignorar la distribución de los datos y usar pruebas paramétricas si estamos tratando con muestras de gran tamaño.
El teorema del límite central nos dice que sin importar el tipo de distribución que tengan las cosas, la distribución de la muestra tiende a ser normal si ésta es lo suficientemente grande (n>30).
Sin embargo, para ser constantes, la normalidad puede ser comprobada mediante:
- Inspección visual.
- Gráficos de normalidad (histogramas). Nos dan una imagen sobre si la distribución de los datos tienen forma de campana.
- Gráficos Q-Q (cuartil-cuartil). También conocidos como gráficos de papel probabilístico, dibujan la correlación entre una muestra y la distribución normal.
- Gráficos de densidad.
- Comprobación de la curtosis y la skewness (asimetría).
- Mediante pruebas de significación (Shapiro-Wilk o Kolmogorov-Smirnoff). Las pruebas de significación. La inspeción visual es siempre recomendable, pero no es fiable. Lo mejor es usar un test de significación para comparar la distribución de la muestra con una normal, para determinar así si los datos muestran una desviación seria de la normalidad. Las pruebas comprueban la distribución de la muestra con una normal para determinar si los datos tienen o no, una desviación seria de la normal.La Hipótesis nula (o de partida) de esta pruebas es que “la distribución de los datos de la muestra es normal”. Si la prueba resulta significativa, la distribución es no-normal.
Es importante incidir que la normalidad está influida por el tamaño de la muestra. Muestras pequeña normalmente pasan la prueba. Sin embargo es importante combinar las inspección visual con las pruebas de significación para tomar la decisión correcta.
Las pruebas de normalidad y el resto de supuestos deberían comprobarse antes de continuar con la prueba principal. Sin embargo, en investigación médica por ejemplo, la normalidad de los datos es más la excepción que la regla. En esas situaciones, se recomienda usar métodos no-paramétricos, como la prueba Wilcoxon para dos muestras, en detrimento de los paramétricos.
4.2 Ejemplo para la comprobación de la normalidad partiendo de vectores.
Generamos algunos vectores.
set.seed(123)
Datos.Normales <- rnorm(n = 100, mean = 0, sd = 1)
Datos.Cola.Izquierda <- rbeta(100,5,2)
Datos.Cola.Derecha <- rbeta(100,2,5)4.2.0.1 Método visual, con vectores.
par(mfrow = c(2,3))
hist(Datos.Normales, breaks = 10, prob = TRUE) #Una indicación del número de columnas a dibujar la podemos obtener haciendo la raíz cuadrada del número de observaciones. En este caso, la raíz cuadrada de 100 es 10.
lines(density(Datos.Normales), lwd = 4, col = "chocolate3")
hist(Datos.Cola.Izquierda, breaks = 10, prob = TRUE)
lines(density(Datos.Cola.Izquierda), lwd = 4, col = "chocolate3")
hist(Datos.Cola.Derecha, breaks = 10, prob = TRUE)
lines(density(Datos.Cola.Derecha), lwd = 4, col = "chocolate3")
qqnorm(Datos.Normales)
qqline(Datos.Normales)
qqnorm(Datos.Cola.Izquierda)
qqline(Datos.Cola.Izquierda)
qqnorm(Datos.Cola.Derecha)
qqline(Datos.Cola.Derecha)
4.2.0.2 Comprobación de la curtosis y la skewness.
Curtosis o achatamiento. En teoría nos da una idea de lo picuda que es la curva de la distribución. El valor de la kurtosis se hace mayor cuanto más picuda es la curva. Las distribuciones se describen como Mesocurticas, Leptocurticas o Platicurtica si sus kurtosis son cercanas a 3 (normal), superiores o inferiores, respectivamente.
Ejemplos de curvas con distinta curtosis.
if (!require(ggplot2)) {install.packages("ggplot2")}
library(ggplot2)ggplot(data.frame(x = c(-4, 4)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = 1), aes(colour = "Mesocurtica"), size = 2) +
stat_function(fun = dnorm, args = list(mean = 0, sd = 0.5), aes(colour = "Leptocurtica"), size = 2) +
stat_function(fun = dnorm, args = list(mean = 0, sd = 2), aes(colour = "Platicurtica"), size = 2) +
scale_colour_manual("Curtosis", values = c("black", "blue", "red"))
Para que en la descripción de los datos aparezca la curtosis y la skewness.
if (!require(psych)) {install.packages("psych")}
library(psych)describe(Datos.Normales)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 100 0.09 0.91 0.06 0.08 0.89 -2.31 2.19 4.5 0.06 -0.22
## se
## X1 0.09#Trimmed mean: La media recortada, utilizada en datos preferentemente simétricos, con muchas observaciones anómalas. Se obtiene eliminando un determinado porcentaje de los datos menores y mayores. Skewness o simetría. En teoría os da una idea de lo simétrica que es nuestra distribución.
par(mfrow = c(1,2))
plot(density(Datos.Cola.Izquierda), main = "Skew negativo", lwd = 10)
plot(density(Datos.Cola.Derecha), main = "Skew positivo", lwd = 10)
par(mfrow = c(1,1))describe(Datos.Cola.Izquierda) #Skewness negativo, cola hacia la izquierda.## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 100 0.74 0.14 0.76 0.75 0.15 0.29 0.97 0.68 -0.73 0.06
## se
## X1 0.01skew(Datos.Cola.Izquierda)## [1] -0.7317864.2.0.3 Prueba de normalidad. Shapiro-Wilk’s test con vectores.
Como el resto de tests que analizan la normalidad, se considera como hipótesis nula que los datos sí proceden de una distribución normal, la hipótesis alternativa que no lo hacen.
El p-value de estos test indica la probabilidad de obtener una distribución como la observada si los datos realmente proceden de una población con una distribución normal.
Este test se emplea para contrastar normalidad cuando el tamaño de la muestra es menor de 50 (¿30?). Para muestras grandes es equivalente al test de Kolmogorov-Smirnov.
shapiro.test(Datos.Normales)##
## Shapiro-Wilk normality test
##
## data: Datos.Normales
## W = 0.99388, p-value = 0.9349shapiro.test(Datos.Cola.Izquierda)##
## Shapiro-Wilk normality test
##
## data: Datos.Cola.Izquierda
## W = 0.95288, p-value = 0.001291Dado que el p-value de Datos.Normales es mayor que el nivel de sinificación alfa = 0.05, podemos deducir que los datos no difieren significativamente de la distribución normal. En otras palabras, podemos asumir la normalidad. Sin embargo, el p-value de Datos.Cola.Iquierda es menor que 0.05, así que rechazamos la hipótsis nula y determinamos que los datos de esta muestra no siguen una ditribución normal.
4.3 Ejemplo para la comprobación de la normalidad partiendo de un data frame.
En el supuesto de que estuviésemos trabajando con un data frame, podríamos hacer los siguiente. Vamos a trabajar con el conjunto de datos (data set) PlantGrowth que está dentro del paquete “agricolae”.
if (!require(agricolae)) {install.packages("agricolae")}
library(agricolae)Vemos la estructura del data set PlantGrowth.
str(PlantGrowth) #Vemos que tiene dos variables "weight" (peso seco planta) y "group" (tratamientos). Mäs información sobre el data set ?PlantGrowth## 'data.frame': 30 obs. of 2 variables:
## $ weight: num 4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
## $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...summary(PlantGrowth) #Los tres grupos tienen el mismo número de datos## weight group
## Min. :3.590 ctrl:10
## 1st Qu.:4.550 trt1:10
## Median :5.155 trt2:10
## Mean :5.073
## 3rd Qu.:5.530
## Max. :6.310describeBy(PlantGrowth$weight, PlantGrowth$group) #Estadísticos diferenciados por grupo##
## Descriptive statistics by group
## group: ctrl
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 10 5.03 0.58 5.15 5 0.72 4.17 6.11 1.94 0.23 -1.12
## se
## X1 0.18
## --------------------------------------------------------
## group: trt1
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 10 4.66 0.79 4.55 4.62 0.53 3.59 6.03 2.44 0.47 -1.1
## se
## X1 0.25
## --------------------------------------------------------
## group: trt2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 10 5.53 0.44 5.44 5.5 0.36 4.92 6.31 1.39 0.48 -1.16
## se
## X1 0.14Como estamos viendo los requisitos previos que se deben cumplir para hacer un T-test de dos muestras, vamos a seleccionar sólo dos grupos, el de control “ctrl” y el del Tratamiento 1 “trt1”.
#PlantGrowth2 <- subset(PlantGrowth, group != "trt2") ; PlantGrowth2
PlantGrowth2 <- subset(PlantGrowth, group == "ctrl" | group == "trt1") ; PlantGrowth2## weight group
## 1 4.17 ctrl
## 2 5.58 ctrl
## 3 5.18 ctrl
## 4 6.11 ctrl
## 5 4.50 ctrl
## 6 4.61 ctrl
## 7 5.17 ctrl
## 8 4.53 ctrl
## 9 5.33 ctrl
## 10 5.14 ctrl
## 11 4.81 trt1
## 12 4.17 trt1
## 13 4.41 trt1
## 14 3.59 trt1
## 15 5.87 trt1
## 16 3.83 trt1
## 17 6.03 trt1
## 18 4.89 trt1
## 19 4.32 trt1
## 20 4.69 trt14.3.0.1 Método visual, con data frame.
Hacemos el análisis visual de cada variable agrupada, usando los paquetes “dplyr” y “ggplot2”.
Calculamos el intercepto y la pendiente para cada grupo.
if (!require(dplyr)) {install.packages("dplyr")}
library(dplyr)intptemedia <- PlantGrowth2 %>% group_by(group) %>%
summarize(q25 = quantile(weight,0.25),
q75 = quantile(weight,0.75),
norm25 = qnorm( 0.25),
norm75 = qnorm( 0.75),
pendiente = (q25 - q75) / (norm25 - norm75),
int = q25 - pendiente * norm25,
media = mean(weight)) %>%
select(group, pendiente, int, media)
intptemedia## # A tibble: 2 x 4
## group pendiente int media
## <fctr> <dbl> <dbl> <dbl>
## 1 ctrl 0.5504161 4.92125 5.032
## 2 trt1 0.4911120 4.53875 4.661Hacemos la gráfica qq para cada grupo.
if (!require(ggplot2)) {install.packages("ggplot2")}
library(ggplot2)Gráficos Q-Q.
Graficos_qq <- ggplot(PlantGrowth2, aes(sample = weight, color = group)) +
stat_qq(distribution = qnorm, size = 5) +
geom_abline(data = intptemedia, aes(intercept = int, slope = pendiente), col = "black", size = 3) +
facet_wrap(~group, nrow = 1) + ylab("Weight")
Graficos_qq
Histogramas.
Graficos_hist <- ggplot(PlantGrowth2, aes(x = weight, fill = group)) +
geom_histogram(aes(y = ..density..), binwidth = .5, colour = "black", alpha = .5) +
geom_vline(data = intptemedia, aes(xintercept = media, colour = group), # Ignore NA values for mean
linetype = "dashed", size = 3) +
geom_density(alpha = .2) +
facet_wrap(~group, nrow = 1)
Graficos_hist
Organizamos los gráficos para que queden unos encima de otros.
if (!require(cowplot)) {install.packages("cowplot")}
library(cowplot)plot_grid(Graficos_qq, Graficos_hist,labels = c("Q-Q", "Histogramas"), label_size = 10, ncol = 1, align = 'v')
4.3.0.2 Prueba de normalidad. Shapiro-Wilk’s test con data frame.
Ejecutamos el test de Shapiro-Wilk para los datos de “weight” según el “group” “ctrl”.
shapiro.test(PlantGrowth2$weight[PlantGrowth2$group == "ctrl"])##
## Shapiro-Wilk normality test
##
## data: PlantGrowth2$weight[PlantGrowth2$group == "ctrl"]
## W = 0.95668, p-value = 0.7475Ejecutamos el test de Shapiro-Wilk para el conjunto de datos de cada variable agrupada.
aggregate(formula = weight ~ group,
data = PlantGrowth2,
FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})## group weight.W weight.V2
## 1 ctrl 0.9566815 0.7474734
## 2 trt1 0.9304107 0.45194404.3.0.3 Prueba de normalidad. Kolmogorov-Smirnov y modificación de Lillefors, con data frame.
Este test se emplea para contrastar normalidad cuando el tamaño de la muestra es mayor de 50 (¿30?). Para muestras grandes es equivalente al test de Shapiro-Wilk.
A pesar de que continuamente se alude al test Kolmogorov-Smirnov como un test válido para contrastar la normalidad, esto no es del todo cierto. El Kolmogorov-Smirnov asume que se conocen la media y varianza poblacional, lo que en la mayoría de los casos no es posible. Esto hace que el test sea muy conservador y poco potente. Para solventar este problema se desarrolló una modificación del Kolmogorov-Smirnov conocida como test Lilliefors. El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para testar la normalidad.
Ejecutamos el test de Kolmogorov-Smirnov con la modificación de Lillefors para el conjunto total de datos.
if (!require(nortest)) {install.packages("nortest")}
library(nortest)Ejecutamos el test de Shapiro-Wilk para los datos de “weight” según el “group” “ctrl”.
lillie.test(PlantGrowth2$weight[PlantGrowth2$group == "ctrl"])##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: PlantGrowth2$weight[PlantGrowth2$group == "ctrl"]
## D = 0.17347, p-value = 0.535aggregate(formula = weight ~ group,
data = PlantGrowth2,
FUN = function(x) {y <- lillie.test(x); c(y$statistic, y$p.value)})## group weight.D weight.V2
## 1 ctrl 0.1734716 0.5350235
## 2 trt1 0.1864702 0.41713324.4 Consecuencias de la falta de normalidad.
Es importante resaltar que cuando se habla de que normalidad influye principalmente en los test de hipótesis paramétricos (t-test, anova,…), esa distribución normal no se refiere a los datos de nuestro esayo, los resíduos de nuestro modelo de ajuste.
En el caso de que las muestras no se distribuyan de forma normal pero se tenga certeza de que las poblaciones de origen sí lo hacen, entonces los resultados obtenidos por los contrastes paramétricos sí son válidos. El teorema del límite central, permite reducir los requerimientos de normalidad cuando las muestras suficientemente grandes.
Si no se cumple esta premisa hay que recurir a al test no-paramétrico de Wilcoxon.