Pruebas de Normalidad para los residuos
Econometria UES
Importación de Datos
library(readr)
ejemplo_regresion_1_ <- read_csv("C:/Users/Genovez/Downloads/ejemplo_regresion (1).csv")
head(ejemplo_regresion_1_)## # A tibble: 6 × 3
## X1 X2 Y
## <dbl> <dbl> <dbl>
## 1 3.92 7298 0.75
## 2 3.61 6855 0.71
## 3 3.32 6636 0.66
## 4 3.07 6506 0.61
## 5 3.06 6450 0.7
## 6 3.11 6402 0.72
Modelo Estimado
library(stargazer)
modelo_estimado<-lm("Y~.",data=ejemplo_regresion_1_)
stargazer(modelo_estimado,title = "Modelo para Ejemplo", type = "html")| Dependent variable: | |
| NA | |
| X1 | 0.237*** |
| (0.056) | |
| X2 | -0.0002*** |
| (0.00003) | |
| Constant | 1.564*** |
| (0.079) | |
| Observations | 25 |
| R2 | 0.865 |
| Adjusted R2 | 0.853 |
| Residual Std. Error | 0.053 (df = 22) |
| F Statistic | 70.661*** (df = 2; 22) |
| Note: | p<0.1; p<0.05; p<0.01 |
Ajuste de los residuos a la Distribución Normal
Verificando el ajuste de los residuos a la distribución normal, se usará la librería fitdistrplus
library(fitdistrplus)
fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm")
plot(fit_normal) ## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 1.804926e-18 0.010000382
## sd 5.000191e-02 0.007058615
## Loglikelihood: 39.41889 AIC: -74.83778 BIC: -72.40002
## Correlation matrix:
## mean sd
## mean 1.000000e+00 -9.796171e-14
## sd -9.796171e-14 1.000000e+00
1. Prueba de Normalidad de Jarque Bera
Usando tseries
##
## Jarque Bera Test
##
## data: modelo_estimado$residuals
## X-squared = 0.93032, df = 2, p-value = 0.628
library(fastGraph)
alpha_sig<-0.05
JB<-salida_JB$statistic
gl<-salida_JB$parameter
VC<-qchisq(1-alpha_sig,gl,lower.tail = TRUE)
shadeDist(JB,ddist = "dchisq",
parm1 = gl,
lower.tail = FALSE,xmin=0,
sub=paste("VC:",round(VC,2)," ","JB:",round(JB,2)))2. Prueba de Kolmogorov Smirnov -Lilliefors
Cálculo Manual
library(dplyr) # Carga la librería dplyr para manipulación de datos
library(gt) # Carga la librería gt para crear tablas de datos
library(gtExtras) # Carga la librería gtExtras para agregar funcionalidades a las tablas creadas con gt
residuos<-modelo_estimado$residuals # Crea un vector con los residuos del modelo estimado
residuos %>% # Utiliza el operador %>% para encadenar las operaciones siguientes al vector residuos
as_tibble() %>% # Convierte el vector residuos en una tibble (tabla) de una columna
mutate(posicion=row_number()) %>% # Agrega una columna llamada "posicion" con el número de fila
arrange(value) %>% # Ordena la tabla por los valores de residuos en orden ascendente
mutate(dist1=row_number()/n()) %>% # Agrega una columna "dist1" con los percentiles según su posición en la tabla (usando la función row_number() y n() para obtener el número de filas)
mutate(dist2=(row_number()-1)/n()) %>% # Agrega una columna "dist2" con los percentiles según su posición en la tabla, pero ajustando en una unidad para evitar problemas con los extremos de la distribución
mutate(zi=as.vector(scale(value,center=TRUE))) %>% # Agrega una columna "zi" con los valores de residuos escalados para tener media cero y varianza uno
mutate(pi=pnorm(zi,lower.tail = TRUE)) %>% # Agrega una columna "pi" con los valores de la función de distribución acumulada (CDF) de una distribución normal estándar evaluada en los valores de zi
mutate(dif1=abs(dist1-pi)) %>% # Agrega una columna "dif1" con las diferencias absolutas entre los percentiles según la posición y los valores de pi
mutate(dif2=abs(dist2-pi)) %>% # Agrega una columna "dif2" con las diferencias absolutas entre los percentiles ajustados según la posición y los valores de pi
rename(residuales=value) -> tabla_KS # Renombra la columna "value" como "residuales" y asigna la tabla resultante a la variable tabla_KS
#Formato
tabla_KS %>% # Utiliza el operador %>% para encadenar las operaciones siguientes a la tabla tabla_KS
gt() %>% # Crea una tabla con la función gt()
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( # Cambia el estilo de algunas celdas de la tabla
style = list(
cell_fill(color = "#A569BD"), # Cambia el color de fondo de las celdas a un tono de morado
cell_text(style = "italic") # Cambia el estilo de texto de las celdas a itálico
),
locations = cells_body( # Aplica el estilo a las celdas del cuerpo de la tabla que cumplan las siguientes condiciones:
columns = dif1, # Que pertenezcan a la columna "dif1"
rows = dif1==max(dif1) # Que pertenezcan a la fila donde el valor de "dif1" es máximo
)) %>%
tab_style( # Cambia el estilo de algunas celdas de la tabla
style = list(
cell_fill(color = "#3498DB"), # Cambia el color de fondo de las celdas a un tono de azul
cell_text(style = "italic") # Cambia el estilo de texto de las celdas a itálico
),
locations = cells_body( # Aplica el estilo a las celdas del cuerpo de la tabla que cumplan las siguientes condiciones:
columns = dif2, # Que pertenezcan a la columna "dif2"
rows = dif2==max(dif2) # Que pertenezcan a la fila donde el valor de "dif2" es máximo
))| Tabla para calcular el Estadistico KS | |||||||
| residuales | posicion | dist1 | dist2 | zi | pi | dif1 | dif2 |
|---|---|---|---|---|---|---|---|
| -0.085089833 | 19 | 0.04 | 0.00 | -1.66734963 | 0.04772245 | 0.007722452 | 0.047722452 |
| -0.069637086 | 22 | 0.08 | 0.04 | -1.36455045 | 0.08619719 | 0.006197193 | 0.046197193 |
| -0.062183196 | 4 | 0.12 | 0.08 | -1.21849022 | 0.11151887 | 0.008481131 | 0.031518869 |
| -0.057379932 | 25 | 0.16 | 0.12 | -1.12436944 | 0.13042817 | 0.029571834 | 0.010428166 |
| -0.057277771 | 12 | 0.20 | 0.16 | -1.12236757 | 0.13085309 | 0.069146907 | 0.029146907 |
| -0.048434733 | 11 | 0.24 | 0.20 | -0.94908676 | 0.17128824 | 0.068711757 | 0.028711757 |
| -0.039102258 | 3 | 0.28 | 0.24 | -0.76621535 | 0.22177409 | 0.058225910 | 0.018225910 |
| -0.033971640 | 16 | 0.32 | 0.28 | -0.66568001 | 0.25280783 | 0.067192174 | 0.027192174 |
| -0.024886579 | 17 | 0.36 | 0.32 | -0.48765671 | 0.31289651 | 0.047103492 | 0.007103492 |
| -0.022166535 | 13 | 0.40 | 0.36 | -0.43435700 | 0.33201461 | 0.067985390 | 0.027985390 |
| -0.018597878 | 8 | 0.44 | 0.40 | -0.36442855 | 0.35776901 | 0.082230992 | 0.042230992 |
| -0.010121525 | 21 | 0.48 | 0.44 | -0.19833299 | 0.42139227 | 0.058607729 | 0.018607729 |
| -0.003341163 | 2 | 0.52 | 0.48 | -0.06547065 | 0.47389964 | 0.046100362 | 0.006100362 |
| 0.012424659 | 6 | 0.56 | 0.52 | 0.24346329 | 0.59617674 | 0.036176739 | 0.076176739 |
| 0.016240338 | 5 | 0.60 | 0.56 | 0.31823217 | 0.62484558 | 0.024845584 | 0.064845584 |
| 0.017281530 | 20 | 0.64 | 0.60 | 0.33863450 | 0.63255746 | 0.007442545 | 0.032557455 |
| 0.026439478 | 10 | 0.68 | 0.64 | 0.51808602 | 0.69780087 | 0.017800874 | 0.057800874 |
| 0.026966239 | 18 | 0.72 | 0.68 | 0.52840800 | 0.70139191 | 0.018608095 | 0.021391905 |
| 0.030236216 | 7 | 0.76 | 0.72 | 0.59248375 | 0.72323665 | 0.036763350 | 0.003236650 |
| 0.035565142 | 15 | 0.80 | 0.76 | 0.69690497 | 0.75706887 | 0.042931135 | 0.002931135 |
| 0.040753758 | 14 | 0.84 | 0.80 | 0.79857675 | 0.78773206 | 0.052267936 | 0.012267936 |
| 0.072518915 | 23 | 0.88 | 0.84 | 1.42102038 | 0.92234458 | 0.042344582 | 0.082344582 |
| 0.073469743 | 1 | 0.92 | 0.88 | 1.43965200 | 0.92501706 | 0.005017061 | 0.045017061 |
| 0.074601871 | 24 | 0.96 | 0.92 | 1.46183625 | 0.92810696 | 0.031893044 | 0.008106956 |
| 0.105692240 | 9 | 1.00 | 0.96 | 2.07105726 | 0.98082328 | 0.019176722 | 0.020823278 |
| Fuente: Elaboración propia | |||||||
Usando nortest
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_estimado$residuals
## D = 0.082345, p-value = 0.9328
En este caso dado que 0.9327912 >0.05 No se rechaza la Hipótesis Nula: ϵ∼N(0,σ2) , por lo que los residuos siguen una distribución normal.
3. Prueba de Shapiro - Wilk
Cálculo Manual
library(dplyr)
library(gt)
residuos<-modelo_estimado$residuals
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
m<-sum(tabla_SW$mi^2)
n<-nrow(ejemplo_regresion_1_)
theta<-1/sqrt(n)
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)
tabla_SW$ai[1]<- -tabla_SW$ai[n]
tabla_SW$ai[2]<- -tabla_SW$ai[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)
tabla_SW$ai[3:(n-2)]<-tabla_SW$mi[3:(n-2)]/sqrt(omega)
tabla_SW %>%
mutate(ai_ui=ai*residuales,ui2=residuales^2) ->tabla_SW
tabla_SW %>%
gt() %>% tab_header("Tabla para calcular el Estadistico W") %>% # Agrega un encabezado a la tabla
tab_source_note(source_note = "Fuente: Elaboración propia")| Tabla para calcular el Estadistico W | |||||
| residuales | pi | mi | ai | ai_ui | ui2 |
|---|---|---|---|---|---|
| -0.085089833 | 0.02475248 | -1.96421684 | -0.44179982 | 0.0375926725 | 7.240280e-03 |
| -0.069637086 | 0.06435644 | -1.51919718 | -0.31084122 | 0.0216460772 | 4.849324e-03 |
| -0.062183196 | 0.10396040 | -1.25930333 | -0.25447584 | 0.0158241211 | 3.866750e-03 |
| -0.057379932 | 0.14356436 | -1.06444156 | -0.21509882 | 0.0123423556 | 3.292457e-03 |
| -0.057277771 | 0.18316832 | -0.90335666 | -0.18254732 | 0.0104559034 | 3.280743e-03 |
| -0.048434733 | 0.22277228 | -0.76286392 | -0.15415701 | 0.0074665537 | 2.345923e-03 |
| -0.039102258 | 0.26237624 | -0.63603674 | -0.12852819 | 0.0050257427 | 1.528987e-03 |
| -0.033971640 | 0.30198020 | -0.51871372 | -0.10481995 | 0.0035609055 | 1.154072e-03 |
| -0.024886579 | 0.34158416 | -0.40814351 | -0.08247629 | 0.0020525526 | 6.193418e-04 |
| -0.022166535 | 0.38118812 | -0.30236184 | -0.06110028 | 0.0013543815 | 4.913553e-04 |
| -0.018597878 | 0.42079208 | -0.19986756 | -0.04038857 | 0.0007511418 | 3.458811e-04 |
| -0.010121525 | 0.46039604 | -0.09943603 | -0.02009370 | 0.0002033789 | 1.024453e-04 |
| -0.003341163 | 0.50000000 | 0.00000000 | 0.00000000 | 0.0000000000 | 1.116337e-05 |
| 0.012424659 | 0.53960396 | 0.09943603 | 0.02009370 | 0.0002496574 | 1.543722e-04 |
| 0.016240338 | 0.57920792 | 0.19986756 | 0.04038857 | 0.0006559241 | 2.637486e-04 |
| 0.017281530 | 0.61881188 | 0.30236184 | 0.06110028 | 0.0010559063 | 2.986513e-04 |
| 0.026439478 | 0.65841584 | 0.40814351 | 0.08247629 | 0.0021806299 | 6.990460e-04 |
| 0.026966239 | 0.69801980 | 0.51871372 | 0.10481995 | 0.0028265997 | 7.271780e-04 |
| 0.030236216 | 0.73762376 | 0.63603674 | 0.12852819 | 0.0038862062 | 9.142287e-04 |
| 0.035565142 | 0.77722772 | 0.76286392 | 0.15415701 | 0.0054826161 | 1.264879e-03 |
| 0.040753758 | 0.81683168 | 0.90335666 | 0.18254732 | 0.0074394891 | 1.660869e-03 |
| 0.072518915 | 0.85643564 | 1.06444156 | 0.21509882 | 0.0155987331 | 5.258993e-03 |
| 0.073469743 | 0.89603960 | 1.25930333 | 0.25447584 | 0.0186962745 | 5.397803e-03 |
| 0.074601871 | 0.93564356 | 1.51919718 | 0.31084122 | 0.0231893370 | 5.565439e-03 |
| 0.105692240 | 0.97524752 | 1.96421684 | 0.44179982 | 0.0466948120 | 1.117085e-02 |
| Fuente: Elaboración propia | |||||
Cálculo del Wny 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)
Wn<-(log(1-W)-mu)/sigma
print(Wn)## [1] -0.3726406
## [1] 0.645292
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)
##
## Shapiro-Wilk normality test
##
## data: modelo_estimado$residuals
## W = 0.97001, p-value = 0.6453
Mismos resultados que en el cálculo manual.
Importante, a partir de esta salida se puede calcular el Wn si se llegará a necesitar:
## [1] -0.3726406