Aplicación de pruebas de normalidad y pruebas de bondad y ajuste

Se realizó la medicion del peso de 206 personas. Se tiene el interés de saber si esta variable proviene de una población con distribución normal con media 10 y desviación estandar igual a 3. Realizar la evaluación utilizando distintas pruebas de bondad y ajuste y la evaluación a través de gráficos.

Validación de supuestos

La prueba de bondad de ajuste requiere dos supuestos. El primero es que la variable sea continua, lo cual en este caso se cumple y segundo que la distribución de los datos sea normal. \

Se realizan los siguientes pasos:

Realizar pruebas para definir si la variable tiene distribución normal
Realizar la estandarización de la variable y dibujar la ecdf
Dibujar la cdf de la distribución teórica
Realizar la prueba de bondad y ajuste de KS para definir si la variable procede de una distribución normal con los parámetros definidos

Pruebas y gráficos de normalidad

Utilizando las herramientas del paquete Graphics

# Histograma
hist(pesos, main="Distribución de los pesos")

# Gráfico de caja 
boxplot(pesos, main="Distribución de los pesos")

# Estandarización del peso 
z.pesos<-(pesos-mean(pesos))/sd(pesos)
# Gráfico QQ
qqnorm(z.pesos)
abline(0,1, col="Red")

# Gráfica de la ecdf 
plot(ecdf(pesos), col ="#db6a28",main="Distribución empírica")
grid(col="Darkgray")

# Gráfica de la normal con parámetros n(10,3)
curve(pnorm(x,10,3),0,20, col = "Blue", main="Distribución teórica")

# Ahora las dos curvas en una gráfica, valores sin estandarizar 
plot(ecdf(pesos), col ="#db6a28",main="Distribuciones")
curve(pnorm(x,10,3),0,20, col = "Blue", add=T)
grid(col="Darkgray")
text(20,0.4,"Empírica", col="#db6a28", cex=0.8)
text(20,0.2,"Teórica", col="Blue", cex=0.8)

Ahora usando ggplot

## Estandarización de la variable peso 
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.2.1     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.4
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

media_peso<-mean(pesos)
sd_peso<-sd(pesos)
peso_y<-(pesos-media_peso)/sd_peso
peso_x<-as.data.frame(pesos)

## Datos sin estandarizar
ggplot(peso_x,aes(x=pesos))+
stat_ecdf(geom="step")+
stat_function(fun=pnorm,color="blue",args = list(10,3))+
labs(title="Empírica_normal vs teórica", y="densidad",x="x")

## Datos estandarizados
peso_estan<-(pesos-media_peso)/sd_peso
peso_estan<-as.data.frame(peso_estan)

ggplot(peso_estan,aes(x=pesos))+
stat_ecdf(geom="step")+
stat_function(fun=pnorm,color="blue",args = list(10,3))+
labs(title="Empírica_normal vs teórica", y="densidad",x="x")

Ahora se realizan las pruebas de normalidad. Se aplican varias pruebas de dos librerías diferentes. Primero usando la libreria nortest

library(nortest)
## Usando la libreria nortest 
usa_nortest<-function(x){
  anderson<-ad.test(x)
  cramer<-cvm.test(x)
  lillies<-lillie.test(x)
  shapiro.francia<-sf.test(x)
  ## Salida a pantalla 
 print(anderson)
 print(cramer)
 print(lillies)
 print(shapiro.francia)
}
usa_nortest(pesos)

## 
##  Anderson-Darling normality test
## 
## data:  x
## A = 0.82594, p-value = 0.03251
## 
## 
##  Cramer-von Mises normality test
## 
## data:  x
## W = 0.10697, p-value = 0.08954
## 
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x
## D = 0.045054, p-value = 0.39
## 
## 
##  Shapiro-Francia normality test
## 
## data:  x
## W = 0.96819, p-value = 0.0002705

Ahora usando el paquete normtest

library(normtest)
usa_normtest<-function(x){
  jarque_a<-ajb.norm.test(x)
  frosini<-frosini.norm.test(x)
  geary<-geary.norm.test(x)
  heg_1<-hegazy1.norm.test(x)
  heg_2<-hegazy2.norm.test(x)
  jarque_b<-jb.norm.test(x)
  kurt<-kurtosis.norm.test(x)
  skek<-skewness.norm.test(x)
  spieg<-spiegelhalter.norm.test(x)
  wb<-wb.norm.test(x)
  ## Salida de datos 
  print(jarque_a)
  print(jarque_b)
  print(frosini)
  print(geary)
  print(heg_1)
  print(heg_2)
  print(kurt)
  print(skek)
  print(wb)
}
usa_normtest(pesos)

## 
##  Adjusted Jarque-Bera test for normality
## 
## data:  x
## AJB = 27.7, p-value = 5e-04
## 
## 
##  Jarque-Bera test for normality
## 
## data:  x
## JB = 25.493, p-value = 5e-04
## 
## 
##  Frosini test for normality
## 
## data:  x
## B = 0.29117, p-value = 0.0425
## 
## 
##  Geary test for normality
## 
## data:  x
## d = 0.80206, p-value = 0.4185
## 
## 
##  Hegazy-Green test for normality
## 
## data:  x
## T = 0.097174, p-value = 0.011
## 
## 
##  Hegazy-Green test for normality
## 
## data:  x
## T = 0.033492, p-value < 2.2e-16
## 
## 
##  Kurtosis test for normality
## 
## data:  x
## T = 4.3557, p-value = 0.0025
## 
## 
##  Skewness test for normality
## 
## data:  x
## T = 0.53201, p-value = 0.002
## 
## 
##  Weisberg-Bingham test for normality
## 
## data:  x
## WB = 0.96819, p-value < 2.2e-16

Prueba de bondad de ajuste

El estadístico D de KS es el valor de la distancia máxima vertical entre la función empírica y teórica. La expresión que representa este valor es:

\(D = max(|S_n(x)-F_n(x)|)\)

Gráficamente es dificil definir cual es el punto porque puede que además no este representado en los datos muestrales.

## Aproximación exacta 
ks.test(pesos, "pnorm", 10,3, alternative="two.sided", exact = T)

## Warning in ks.test(pesos, "pnorm", 10, 3, alternative = "two.sided", exact = T):
## ties should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pesos
## D = 0.31143, p-value = 3.331e-15
## alternative hypothesis: two-sided

## Aproximación asintótica 
ks.test(pesos, "pnorm", 10,3, alternative="two.sided", exact = F)

## Warning in ks.test(pesos, "pnorm", 10, 3, alternative = "two.sided", exact = F):
## ties should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pesos
## D = 0.31143, p-value < 2.2e-16
## alternative hypothesis: two-sided

Se rechaza la hipótesis nula de que los datos provienen de la distribución normal N(10,3).

Taller KS