Pruebas de normalidad

Universidad Nacional de San Cristóbal de Huamanga

Registros sobre velocidad de la luz de Michelson en 1879: Bases de datos (Crawley, 2013, 2014) http://www.bio.ic.ac.uk/research/mjcraw/therbook/index.htm http://www.bio.ic.ac.uk/research/crawley/statistics/

Pruebas de Normalidad

library(nortest)

luz = read.table("./light.txt", header=T)
luz

##    speed
## 1    850
## 2    740
## 3    900
## 4   1070
## 5    930
## 6    850
## 7    950
## 8    980
## 9    980
## 10   880
## 11  1000
## 12   980
## 13   930
## 14   650
## 15   760
## 16   810
## 17  1000
## 18  1000
## 19   960
## 20   960

names(luz) = "velocidad"
luz

##    velocidad
## 1        850
## 2        740
## 3        900
## 4       1070
## 5        930
## 6        850
## 7        950
## 8        980
## 9        980
## 10       880
## 11      1000
## 12       980
## 13       930
## 14       650
## 15       760
## 16       810
## 17      1000
## 18      1000
## 19       960
## 20       960

Contraste de hipótesis

H0: Los datos tienden una distribución Normal

H1: Los datos no tienden una distribución Normal

Comprobar la normalidad de los datos con pruebas estadisticas:

1. Prueba de Anderson-Darling

Es la derivada del Kolmorof

head(cars)

##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

attach(cars)
ad.test(speed)

## 
##  Anderson-Darling normality test
## 
## data:  speed
## A = 0.26143, p-value = 0.6927

ad.test(dist)

## 
##  Anderson-Darling normality test
## 
## data:  dist
## A = 0.74067, p-value = 0.05021

ad.test(luz$velocidad)

## 
##  Anderson-Darling normality test
## 
## data:  luz$velocidad
## A = 0.67243, p-value = 0.0671

2. Prueba de Cramer-von Mises

Es útil para pequeñas muestras y usa los momentos como criterio.

sd <- c(23,76,2,3,23,2,3,2,3,23,23)
cvm.test(sd)

## 
##  Cramer-von Mises normality test
## 
## data:  sd
## W = 0.23633, p-value = 0.001408

cvm.test(dist)

## 
##  Cramer-von Mises normality test
## 
## data:  dist
## W = 0.12632, p-value = 0.04742

cvm.test(luz$velocidad)

## 
##  Cramer-von Mises normality test
## 
## data:  luz$velocidad
## W = 0.11806, p-value = 0.05815

3. Prueba de Kolmogorov-Smirnov

Kolmogorov-Smirnov … no hay que usarla por ser poco confiable (Steinskog et al., 2007; Pedrosa et al., 2015)

ks.test(luz$velocidad,"pnorm",mean(luz$velocidad),sd(luz$velocidad))

## Warning in ks.test.default(luz$velocidad, "pnorm", mean(luz$velocidad), : ties
## should not be present for the Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  luz$velocidad
## D = 0.17931, p-value = 0.541
## alternative hypothesis: two-sided

Crawley (2013) para más información sobre esta prueba

4. Prueba de Lilliefors (Kolmogorov-Smirnov)

Se aplica mas ampliamente cuando la muestra es grande.

Usar K-S con corrección de Lilliefors en lugar de K-S (Razali & Wah, 2011; Pedrosa et al., 2015):

lillie.test(dist)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  dist
## D = 0.12675, p-value = 0.04335

install.packages("nortest")

## Warning: package 'nortest' is in use and will not be installed

library("nortest")
lillie.test(luz$velocidad)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  luz$velocidad
## D = 0.17931, p-value = 0.09138

Otras pruebas de normalidad (ver Yap & Sim, 2011 y librer?a “nortest”):

5. Prueba de Pearson chi-square

basada en una distribución Ji cuadrado y que corresponde a una prueba de bondad de ajuste.

pearson.test(dist)

## 
##  Pearson chi-square normality test
## 
## data:  dist
## P = 11.2, p-value = 0.1301

pearson.test(luz$velocidad)

## 
##  Pearson chi-square normality test
## 
## data:  luz$velocidad
## P = 8.7, p-value = 0.06905

6. Prueba de Shapiro-Wilk

Es más poderosa cuando se compara con otras pruebas de normalidad cuando la muestra es pequeña.

Prueba de normalidad más robusta (Yap & Sim, 2011; Ghasemi & Zahediasl, 2012). Aunque tiene limitaciones (Yap & Sim, 2011; Razali & Wah, 2011)

shapiro.test(luz$velocidad)

## 
##  Shapiro-Wilk normality test
## 
## data:  luz$velocidad
## W = 0.91992, p-value = 0.09876

“El valor de P es una medida de credibilidad de la hipótesis nula” (Crawley, 2014).

La hipótesis nula de la prueba de Shapiro-Wilk (así como de la mayoría de las pruebas de normalidad) es que la distribución de nuestros datos sea normal.

Para nuestro ejemplo de la velocidad de la luz (registrada por Michelson en 1879), se tiene un porcentaje de credibilidad del 9.9% de que la velocidad de la luz tenga una distribución normal. Puesto que es mayor al 5%, concluimos que tiene una distribución normal.

La manera en como se reportaría, sería algo como lo siguiente: Las medidas registradas sobre la velocidad de la luz por Michelson, en 1879, presentaron una distribución normal (W=0.92; P=0.099) [NOTA: la “P” siempre irá en mayúsculas e itálicas].

7. Prueba de Shapiro-Francia

simplificacion de la prueba shapiro-W

sf.test(dist)

## 
##  Shapiro-Francia normality test
## 
## data:  dist
## W = 0.95206, p-value = 0.04179

sf.test(luz$velocidad)

## 
##  Shapiro-Francia normality test
## 
## data:  luz$velocidad
## W = 0.91666, p-value = 0.07992

8. Prueba de Frosini

library(normtest)
frosini.norm.test(speed)

## 
##  Frosini test for normality
## 
## data:  speed
## B = 0.14841, p-value = 0.734

9. Prueba de Geary

Usa los valores acumulados muestrales, sus medias y desviaciones estándar.

geary.norm.test(speed)

## 
##  Geary test for normality
## 
## data:  speed
## d = 0.82835, p-value = 0.181

10. Prueba de Hegazy-Green

hegazy1.norm.test(dist, nrepl=20000) ###nrepl: considera el número de replicas en simulación de Monte Carlo

## 
##  Hegazy-Green test for normality
## 
## data:  dist
## T = 0.16856, p-value = 0.03545

11. Prueba de Jarque-Bera

Utiliza un estadístico en la prueba que involucra la curtosis y la asimetría. –Usada por economistas.

jb.norm.test(dist, nrepl=2000)

## 
##  Jarque-Bera test for normality
## 
## data:  dist
## JB = 5.2305, p-value = 0.053

12. Prueba de Kurtosis

kurtosis.norm.test(speed, nrepl=2000)

## 
##  Kurtosis test for normality
## 
## data:  speed
## T = 2.4229, p-value = 0.325

13. Prueba de Skewness

skewness.norm.test(speed, nrepl=2000)

## 
##  Skewness test for normality
## 
## data:  speed
## T = -0.11395, p-value = 0.711

14. Prueba de Spiegelhalter

spiegelhalter.norm.test(speed, nrepl=2000)

## 
##  Spiegelhalter test for normality
## 
## data:  speed
## T = 1.2295, p-value = 0.7475

15. Puerba de Weisberg-Bingham

wb.norm.test(speed, nrepl=2000)

## 
##  Weisberg-Bingham test for normality
## 
## data:  speed
## WB = 0.98416, p-value = 0.652

16. Prueba de Agostino

library(moments)
agostino.test(dist)

## 
##  D'Agostino skewness test
## 
## data:  dist
## skew = 0.78248, z = 2.31281, p-value = 0.02073
## alternative hypothesis: data have a skewness

Método gráficos

Comprobar normalidad con gráficas:

Histograma:

hist(luz$velocidad)

Gráfico cuantil-cuantil (Crawley, 2014):

qqnorm(luz$velocidad)
qqline(luz$velocidad)

Si nuestros datos de interés (puntos) están alineados con la línea graficada, quiere decir que tienen una distribución normal. Si por el contrario, observamos que se acomodan en forma de “S” o de plátano (curva o en “U”), no son normales. Para más información consulta Crawley (2013,2014).

Gráfico para probar normalidad de Crawley (2013):

crawley.plot <- function(y) {
  s <- sd(y)
  plot(c(0,3),
       c(min(y,mean(y)-s*4*qnorm(0.75)),max(y,mean(y)+s*4*qnorm(0.75))),
       xaxt="n",
       xlab="",
       type="n",
       ylab="")
       
       
  # for your data's boxes and whiskers, centred at x = 1
  top <- quantile(y,0.75)
  bottom <- quantile(y,0.25)
  w1u <- quantile(y,0.91)
  w2u <- quantile(y,0.98)
  w1d <- quantile(y,0.09)
  w2d <- quantile(y,0.02)
  rect(0.8,bottom,1.2,top)
  lines(c(0.8,1.2),c(mean(y),mean(y)),lty=3)
  lines(c(0.8,1.2),c(median(y),median(y)))
  lines(c(1,1),c(top,w1u))
  lines(c(0.9,1.1),c(w1u,w1u))
  lines(c(1,1),c(w2u,w1u),lty=3)
  lines(c(0.9,1.1),c(w2u,w2u),lty=3)
  nou <- length(y[y>w2u])
  points(rep(1,nou),jitter(y[y>w2u]))
  lines(c(1,1),c(bottom,w1d))
  lines(c(0.9,1.1),c(w1d,w1d))
  lines(c(1,1),c(w2d,w1d),lty=3)
  lines(c(0.9,1.1),c(w2d,w2d),lty=3)
  nod <- length(y[y<w2d])
  points(rep(1,nod),jitter(y[y<w2d]))
  #for the normal box and whiskers, centred at x = 2
  n75 <- mean(y)+ s * qnorm(0.75)
  n25 <- mean(y)- s * qnorm(0.75)
  n91 <- mean(y)+ s * 2* qnorm(0.75)
  n98 <- mean(y)+ s * 3* qnorm(0.75)
  n9 <- mean(y)- s * 2* qnorm(0.75)
  n2 <- mean(y)- s * 3* qnorm(0.75)
  rect(1.8,n25,2.2,n75)
  lines(c(1.8,2.2),c(mean(y),mean(y)),lty=3)
  lines(c(2,2),c(n75,n91))
  lines(c(1.9,2.1),c(n91,n91))
  lines(c(2,2),c(n98,n91),lty=3)
  lines(c(1.9,2.1),c(n98,n98),lty=3)
  lines(c(2,2),c(n25,n9))
  lines(c(1.9,2.1),c(n9,n9))
  lines(c(2,2),c(n9,n2),lty=3)
  lines(c(1.9,2.1),c(n2,n2),lty=3)
  lines(c(1.2,1.8),c(top,n75),lty=3,col="gray")
  lines(c(1.1,1.9),c(w1u,n91),lty=3,col="gray")
  lines(c(1.1,1.9),c(w2u,n98),lty=3,col="gray")
  lines(c(1.2,1.8),c(bottom,n25),lty=3,col="gray")
  lines(c(1.1,1.9),c(w1d,n9),lty=3,col="gray")
  lines(c(1.1,1.9),c(w2d,n2),lty=3,col="gray")
  # label the two boxes
  axis(1,c(1,2),c("data","normal")) }

crawley.plot(luz$velocidad)

A la derecha del grafico (cajas y bigotes) tenemos una distrubución normal creada con la media y las desviación estantdar de la distribución de nuestros datos, en este caso de la variable velocidad de la luz (las lineas de meidana y media se impalman, los rangos intecuantílicos, osea la seccion de los bigotes y la seccion de las cajas son mas o menos del mismo tamaño o iguales, los datos atípicos son cosas poco comun)

a lado izquierdo tenemos nuestros datos sin ningun tipo de transformación (las lineas punteadas y sombreada no se impalman, quiere decir que la linea punteada es la media y linea sólida es la mediana,en este caso se oberva que el bigote inferior es mucho mas grande la la superior, por último nuestros datos presentan datos atípicos)

Podemos ver que la distribución de la velocidad de la luz NO ES NORMAL sim embargo nuestras pruebas estadísticas del inicio dicen que sí, esto se debe porque solamente tenemos 20 valores en nuestra variable, entonces algunos test no son buenos cuando tenemos menos de 30 datos.

Si utilizamos muestras mayores a 30

peces = read.table("./fishes.txt", header=T)

Podemos observar que tenemos mas de 150 datos

hist(peces$mass)

Graficando el histograma observamos que tiene una distribución normal

shapiro.test(peces$mass)

## 
##  Shapiro-Wilk normality test
## 
## data:  peces$mass
## W = 0.91651, p-value = 1.287e-10

lillie.test(peces$mass)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  peces$mass
## D = 0.10448, p-value = 5.865e-07

ad.test(peces$mass)

## 
##  Anderson-Darling normality test
## 
## data:  peces$mass
## A = 4.6184, p-value = 1.777e-11

Haciendo las pruebas de normalidad observamos que los p son mucho menores que el 0,05, que quiere decir que esta variable no tiene una distribución normal

Doctor en estadística Ian Fellows, por parte de la universidad de California: http://blog.fellstat.com/?p=61

x = rt(1000,df=200)

hist(x)

qqnorm(x)
qqline(x)

shapiro.test(x)

## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.99782, p-value = 0.2162

lillie.test(x)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x
## D = 0.023704, p-value = 0.1889

ad.test(x)

## 
##  Anderson-Darling normality test
## 
## data:  x
## A = 0.62939, p-value = 0.1008

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.