Ejercicio de pruebas de Normalidad

Pruebas de Normalidad

library(wooldridge) #Esta línea carga la biblioteca wooldridge, que contiene conjuntos de datos utilizados en el libro "Introductory Econometrics: A Modern Approach" de Jeffrey M. Wooldridge.
data(hprice1) #Esta línea carga el conjunto de datos hprice1en el entorno de trabajo. Este conjunto de datos contiene información sobre casas vendidas en un suburbio de Chicago en 1987, incluyendo el precio de venta, el número de habitaciones, el tamaño del lote y la edad de la casa.
head(force(hprice1),n=5) #Esta línea imprime las primeras cinco filas del conjunto de datos hprice1. La función forcese utiliza para forzar la carga del conjunto de datos si aún no se ha cargado en el entorno de trabajo. La función headse utiliza para seleccionar las primeras cinco filas del conjunto de datos y mostrarlas en la consola.

##   price assess bdrms lotsize sqrft colonial   lprice  lassess llotsize   lsqrft
## 1   300  349.1     4    6126  2438        1 5.703783 5.855359 8.720297 7.798934
## 2   370  351.5     3    9903  2076        1 5.913503 5.862210 9.200593 7.638198
## 3   191  217.7     3    5200  1374        0 5.252274 5.383118 8.556414 7.225482
## 4   195  231.8     3    4600  1448        1 5.273000 5.445875 8.433811 7.277938
## 5   373  319.1     4    6095  2514        1 5.921578 5.765504 8.715224 7.829630

#mostrar las primeras 5 observaciones

1. estime el siguiente modelo,

price = ˆα + ˆα1(lotsize) + ˆα2(sqrft) + ˆα3(bdrms) + E

library(stargazer) #La primera línea carga la librería "stargazer" en la sesión de R:
modelo_estimado<-lm(price ~ lotsize + sqrft + bdrms,data = hprice1) #se está ajustando un modelo que predice la variable priceen función de las variables lotsize, sqrfty bdrms. Los datos que se están utilizando para ajustar el modelo provienen del conjunto de datos hprice1.
stargazer(modelo_estimado,title = "Modelo para Ejemplo", type = "text")

## 
## Modelo para Ejemplo
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                price           
## -----------------------------------------------
## lotsize                      0.002***          
##                               (0.001)          
##                                                
## sqrft                        0.123***          
##                               (0.013)          
##                                                
## bdrms                         13.853           
##                               (9.010)          
##                                                
## Constant                      -21.770          
##                              (29.475)          
##                                                
## -----------------------------------------------
## Observations                    88             
## R2                             0.672           
## Adjusted R2                    0.661           
## Residual Std. Error      59.833 (df = 84)      
## F Statistic           57.460*** (df = 3; 84)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus) #Esta línea carga la biblioteca fitdistrplus, que proporciona funciones para ajustar distribuciones de probabilidad a datos.
fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm") #En esta línea, se ajusta una distribución normal a los residuos del modelo estimado. La función fitdistde fitdistrplusrealizar el ajuste y devuelve un objeto fitdistque contiene los resultados.
plot(fit_normal) #Esta línea traza el ajuste de la distribución normal realizada en la línea anterior. La función plotse usa para visualizar el objeto fitdisty sus resultados.

summary(fit_normal) #Esta línea imprime un resumen de los resultados del ajuste de distribución en la consola. El resumen incluye información como los medios, la desviación estándar y la probabilidad de que los residuos sigan una distribución normal.

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##          estimate Std. Error
## mean 9.992007e-16   6.231624
## sd   5.845781e+01   4.406424
## Loglikelihood:  -482.8775   AIC:  969.7549   BIC:  974.7096 
## Correlation matrix:
##      mean sd
## mean    1  0
## sd      0  1

2. Verique el supuesto de normalidad, a través de:

1. La prueba JB

library(tseries) # Carga la librería tseries, que proporciona funciones para el análisis de series temporales.
salida_jb<-jarque.bera.test(modelo_estimado$residuals) #realiza la prueba de Jarque-Bera sobre los residuos del modelo estimado y almacena los resultados en la variable salida_jb. La prueba de Jarque-Bera es una prueba de normalidad que utiliza la asimetría y la curvatura de una distribución para determinar si la distribución es normal o no.
salida_jb #simplemente imprime los resultados de la prueba.

## 
##  Jarque Bera Test
## 
## data:  modelo_estimado$residuals
## X-squared = 32.278, df = 2, p-value = 9.794e-08

library(fastGraph)
alpha_sig<-0.05 #define el nivel de significancia que se ajusta en la prueba. En este caso, se utiliza un nivel de significancia de 0.05, lo que significa que se rechazará la hipótesis nula si la probabilidad de obtener un estadístico de prueba tan extremo o más extremo que el valor observado es menor que 0.05.
jb<-salida_jb$statistic #almacena el valor del estadístico de prueba de Jarque-Bera en la variable jb.
gl<-salida_jb$parameter #almacena el número de grados de libertad de la distribución chi-cuadrado asociado a la prueba de Jarque-Bera en la variable gl.
VC<-qchisq(1- alpha_sig,gl,lower.tail = TRUE) #calcula el valor crítico para la distribución chi-cuadrado correspondiente a un nivel de significancia de 0.05 y el número de grados de libertad obtenidos en la prueba de Jarque-Bera. Este valor crítico se almacena en la variable VC.
shadeDist(jb,ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE, xmin = 0,
          sub=paste("VC",round(VC,2)," ","jb:",round(jb,2))) #utiliza la función shadeDistde la librería "fastGraph" para graficar la distribución chi-cuadrado con el número de grados de libertad obtenidos en la prueba de Jarque-Bera y sombrear el área de la distribución correspondiente al valor del estadístico de prueba observado . Además, la línea también agrega un subtítulo al gráfico que indica el valor crítico calculado y el valor del estadístico de prueba observado. En resumen, este gráfico ayuda a interpretar visualmente el resultado de la prueba de Jarque-Bera y determinar si los residuos del modelo siguen una distribución normal.

Usando notmtest (descatalogada)

library(nortest) #Esta línea carga la biblioteca nortest, que proporciona funciones para realizar pruebas de normalidad.
jarque.bera.test(modelo_estimado$residuals) #Esta línea realiza una prueba de Jarque-Bera en los residuos del modelo estimado. La prueba de Jarque-Bera es una prueba de bondad de ajuste que se utiliza para determinar si una muestra de datos tiene una distribución normal. La función jarque.bera.testde nortestrealizar la prueba y devuelve un objeto htestque contiene los resultados de la prueba, incluyendo el estadístico de prueba y el valor p.

## 
##  Jarque Bera Test
## 
## data:  modelo_estimado$residuals
## X-squared = 32.278, df = 2, p-value = 9.794e-08

##b) La prueba KS

library(dplyr) # Cargamos la librería dplyr para trabajar con los datos de manera más eficiente
residuos<-modelo_estimado$residuals # Guardamos los residuos del modelo estimado en una variable llamada "residuos"
residuos %>% # Convertimos los residuos en un tibble
  as_tibble() %>% #Convierte el objeto en un tibble (tabla) para que pueda ser manipulado por la biblioteca "dplyr"
  mutate(posicion=row_number()) %>% # Agregamos una columna llamada "posicion" que indica la posición de cada valor en el tibble
  arrange(value) %>% # Ordenamos el tibble en orden ascendente de valores
  mutate(dist1=row_number()/n()) %>% # Agregamos una columna llamada "dist1" que indica la distancia relativa entre la posición de cada valor y la distribución teórica
  mutate(dist2=(row_number()-1)/n()) %>% # Agregamos una columna llamada "dist2" que indica la distancia relativa entre la posición de cada valor menos uno y la distribución teórica
  mutate(zi=as.vector(scale(value,center=TRUE))) %>% # Agregamos una columna llamada "zi" que indica los residuos estandarizados
  mutate(pi=pnorm(zi,lower.tail = TRUE)) %>% # Agregamos una columna llamada "pi" que indica la probabilidad acumulada de los residuos estandarizados en la distribución normal
  mutate(dif1=abs(dist1-pi)) %>% # Agregamos una columna llamada "dif1" que indica la diferencia absoluta entre la distribución teórica y la posición de cada valor
  mutate(dif2=abs(dist2-pi)) %>% # Agregamos una columna llamada "dif2" que indica la diferencia absoluta entre la distribución teórica y la posición de cada valor menos uno
  rename(residuales=value) -> tabla_KS# Renombramos la columna "value" a "residuales"


#forma
library(gt)
library(gtExtras) #Carga las bibliotecas "gt" y "gtExtras", que se utilizan para crear tablas interactivas en R.

tabla_KS %>%
  gt() %>% #Crea una tabla interactiva de la tabla "tabla_KS" utilizando la función "gt" de la biblioteca "gt". El operador "%>%" se utiliza para pasar el objeto "tabla_KS" a la siguiente función.
  tab_header("Tabla para calcular el Estadistico KS") %>%  #Agrega un encabezado a la tabla.
  tab_source_note(source_note = "Fuente: Elaboración propia") %>%  #Agrega una nota de fuente a la tabla.
  tab_style(style = list(
      cell_fill(color = "#82E0AA"),  
      cell_text(style = "italic")), 
    locations = cells_body( 
      columns = dif1,  
      rows = dif1==max(dif1) #Aplica un estilo a las celdas que tienen la columna "dif1" y la fila correspondiente al valor máximo de "dif1". El estilo aplicado es un color de fondo verde claro ("#82E0AA") y texto en cursiva. 
    )) %>%
   tab_style(  
    style = list(
      cell_fill(color = "#CB4335"),  
      cell_text(style = "italic") 
      ),
    locations = cells_body(  
      columns = dif2, 
      rows = dif2==max(dif2)))#Aplica otro estilo a las celdas que tienen la columna "dif2" y la fila correspondiente al valor máximo de "dif2". El estilo aplicado es un color de fondo rojo oscuro ("#CB4335") y texto en cursiva.

###Cálculo del estadístico D:

D<-max(max(tabla_KS$dif1),max(tabla_KS$dif2))
print(D)

## [1] 0.0754392

Valor crítico de la tabla de Lilliefors

n=88 α=0.05

n>50 (0,875897/√n)=0,0933709

Conclusión:

En este caso dado que 0.0754392 < 0,0933709 No se rechaza la Hipótesis Nula: ϵ∼N(0,σ2) , por lo que los residuos siguen una distribución normal.

library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals) #En esta línea, se realiza la prueba de Lilliefors para comprobar si los residuos del modelo estimado siguen una distribución normal. La función lillie.testde la biblioteca nortestrealiza la prueba y devuelve un objeto htestque contiene los resultados de la prueba, como el estadístico de prueba y el valor p.
print(prueba_KS) #Esta línea imprime los resultados de la prueba de Lilliefors en la consola. En particular, la salida obtuvo el valor del estadístico de prueba y el valor p.

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_estimado$residuals
## D = 0.075439, p-value = 0.2496

p.value<-prueba_KS$p.value #Esta línea extrae el valor p del objeto htestgenerado por la prueba de Lilliefors y lo almacena en la variable p.value. El valor p es una medida de la evidencia en contra de la hipótesis nula, que en este caso sería que los residuos siguen una distribución normal. Si el valor p es menor que el nivel de significancia elegido, se rechazaría la hipótesis nula y se concluiría que los residuos no siguen una distribución normal.

library(dplyr)
library(gt) # La biblioteca dplyrproporciona funciones para manipular y transformar datos, mientras que gtes una biblioteca para generar tablas elegantes y bien formateadas.
residuos<-modelo_estimado$residuals #Esta línea extrae los residuos del modelo estimado y los almacena en la variable residuos.

residuos %>%  
  as_tibble() %>%
  rename(residuales=value) %>%
  arrange(residuales) %>%
  mutate(pi=(row_number()-0.375)/(n()+0.25)) %>%
  mutate(mi=qnorm(pi,lower.tail = TRUE)) %>% 
  mutate(ai=0)->tabla_SW #En esta sección de código se realiza un cálculo de los valores necesarios para la prueba de normalidad de Shapiro-Wilk. Primero, los residuos se definen en un tibble utilizando as_tibble(). Luego, se renombra la columna de residuos como "residuales" utilizando rename(). Después, los residuos se ordenan en orden creciente arrange(). Luego, se calculan los valores pi, mi, y ainecesarios para la prueba de Shapiro-Wilk de normalidad. Estos cálculos se realizan utilizando la función mutate(). Finalmente, la tabla quedó asignada a la variable tabla_SW.

m<-sum(tabla_SW$mi^2) #Esta línea calcula la suma de los cuadrados de los valores mide la tabla tabla_SWy lo almacena en la variable m.
n<-nrow(hprice1) #Esta línea calcula el número de filas de la base de datos hprice1(que contiene los datos del modelo) y lo almacena en la variable n.
theta<-1/sqrt(n) #Esta línea calcula el valor de thetaque se reafirma en los cálculos posteriores y lo almacena en la variable theta.


tabla_SW$ai[n]<- -2.706056*theta^5+4.434685*theta^4-2.071190*theta^3-0.147981*theta^2+0.2211570*theta+tabla_SW$mi[n]/sqrt(m)
tabla_SW$ai[n-1]<- -3.582633*theta^5+5.682633*theta^4-1.752461*theta^3-0.293762*theta^2+0.042981*theta+tabla_SW$mi[n-1]/sqrt(m) #Estas dos líneas calculan y asignan valores a los elementos aien las posiciones ny n-1de la tabla tabla_SW. Los cálculos se basan en los valores de theta, mi, y my en constantes específicas del método de Shapiro-Wilk.
tabla_SW$ai[1]<- -tabla_SW$ai[n]
tabla_SW$ai[2]<- -tabla_SW$ai[n-1] #Estas dos líneas calculan y asignan valores a los elementos aien las posiciones 1 y 2 de la tabla tabla_SW. Los cálculos se basan en los valores previamente asignados a los elementos aien las posiciones ny n-1.

omega<-(m-2*tabla_SW$mi[n]^2-2*tabla_SW$mi[n-1]^2)/(1-2*tabla_SW$ai[n]^2-2*tabla_SW$ai[n-1]^2) #Esta línea calcula el valor de omegautilizar los valores previamente calculados de m, mi, y ai.
tabla_SW$ai[3:(n-2)]<-tabla_SW$mi[3:(n-2)]/sqrt(omega) #Esta línea calcula y asigna valores a los elementos aien las posiciones 3 a n-2de la tabla tabla_SW. Los cálculos se basan en los valores previamente calculados de miy omega.

tabla_SW %>% 
  mutate(ai_ui=ai*residuales,ui2=residuales^2) ->tabla_SW
tabla_SW %>% #se crea una nueva columna en la tabla tabla_SWllamada "ai_ui", que corresponde al producto entre la columna "ai" y "residuales". También se crea una columna llamada "ui2", que corresponde al cuadrado de los residuos.
  gt() %>% tab_header("Tabla para calcular el Estadistico W") %>% #se genera una tabla utilizando la función gt(), que muestra las columnas "residuales", "pi", "mi", "ai_ui" y "ui2". Esta tabla está encabezada por el texto "Tabla para calcular el Estadistico W" usando la función tab_header(). 
  tab_source_note(source_note = "Fuente: Elaboración propia") # se agrega una nota de fuente con el texto "Fuente: Elaboración propia" utilizando la función tab_source_note().

W<-(sum(tabla_SW$ai_ui)^2)/sum(tabla_SW$ui2) # se calcula el estadístico de Shapiro-Wilk utilizando la fórmula W = (sum(ai_ui)^2)/sum(ui2), que se almacena en la variable W.
print(W) #se muestra el valor del estadistico utilizando la funcion print().

## [1] 0.9413208

Cálculo del Wn y su p value

mu<-0.0038915*log(n)^3-0.083751*log(n)^2-0.31082*log(n)-1.5861
sigma<-exp(0.0030302*log(n)^2-0.082676*log(n)-0.4803) #se calculan los valores de mu y sigma necesarios para normalizar el estadístico. Los valores de mu y sigma se calculan utilizando las ecuaciones obtenidas por Royston (1982) para muestras de tamaño n >= 20.
Wn<-(log(1-W)-mu)/sigma #Luego, se calcula el valor estandarizado del estadístico de Shapiro-Wilk utilizando la fórmula Wn = (log(1-W)-mu)/sigma, que se almacena en la variable Wn.
print(Wn)

## [1] 3.241867

p.value.sw<-pnorm(Wn,lower.tail = FALSE) #se calcula el valor-p asociado al estadístico normalizado utilizando la función de distribución acumulativa de la normal estándar pnorm(), con el argumento lower.tail=FALSEpara obtener el valor-p de la cola superior. El valor-p se almacena en la variable p.value.sw
print(p.value.sw) #El valor-p se almacena en la variable p.value.swy se muestra utilizando la función print().

## [1] 0.0005937472

library(fastGraph) #llama a la función shadeDist()del paquete fastGraph, que se utiliza para sombrear una distribución de probabilidad a partir de un valor crítico.
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE) #El argumento lower.tail = FALSEindica que se quiere sombrear el area de la cola superior. Si se quisiera sombrear la cola inferior, se usaría lower.tail = TRUE.

En este caso dado que 0.645292 > 0.05 No se rechaza la Hipótesis Nula: ϵ∼N(0,σ2), por lo que los residuos siguen una distribución normal.

Usando la librería stats (precargada al iniciar R)

salida_SW<-shapiro.test(modelo_estimado$residuals) #llama a la función shapiro.test()de R, que se utiliza para realizar la prueba de normalidad de Shapiro-Wilk en una muestra de datos. En este caso, se está aplicando esta prueba a los residuos del modelo estimado.
print(salida_SW) #La salida que se muestra al imprimir el resultado de la función estadística de prueba de Shapiro-Wilk ( W) y su valor-p ( p-value), así como un mensaje que indica si se rechaza o no la hipótesis nula de normalidad.

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_estimado$residuals
## W = 0.94132, p-value = 0.0005937

Si el valor-p es menor que el nivel de significancia (generalmente 0.05), se rechaza la hipótesis nula de normalidad y se concluye que los residuos no siguen una distribución normal. Si el valor-p es mayor que el nivel de significancia, se acepta la hipótesis nula de normalidad y se concluye que los residuos se distribuyen normalmente.

salida_Wn<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(salida_Wn)

## [1] 3.241867