Carga de datos

mtcars
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

Estimación del modelo.

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
modelo_estimado <-lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
stargazer(modelo_estimado, title = "modelo estimado", type = "text")
## 
## modelo estimado
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 mpg            
## -----------------------------------------------
## disp                           0.003           
##                               (0.011)          
##                                                
## hp                            -0.019           
##                               (0.016)          
##                                                
## wt                           -4.609***         
##                               (1.266)          
##                                                
## qsec                           0.544           
##                               (0.466)          
##                                                
## Constant                     27.330***         
##                               (8.639)          
##                                                
## -----------------------------------------------
## Observations                    32             
## R2                             0.835           
## Adjusted R2                    0.811           
## Residual Std. Error       2.622 (df = 27)      
## F Statistic           34.195*** (df = 4; 27)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Diagnostico de Multicolinealidad.

# Calculo Manual

library(stargazer)
X_mat <-model.matrix(modelo_estimado)
stargazer(head(X_mat, n = 6), type = "text")
## 
## ===================================================
##                   (Intercept) disp hp   wt    qsec 
## ---------------------------------------------------
## Mazda RX4              1      160  110 2.620 16.460
## Mazda RX4 Wag          1      160  110 2.875 17.020
## Datsun 710             1      108  93  2.320 18.610
## Hornet 4 Drive         1      258  110 3.215 19.440
## Hornet Sportabout      1      360  175 3.440 17.020
## Valiant                1      225  105 3.460 20.220
## ---------------------------------------------------
XX_matrix <-t(X_mat) %*% X_mat
stargazer(XX_matrix, type = "text")
## 
## ==========================================================================
##             (Intercept)     disp           hp           wt        qsec    
## --------------------------------------------------------------------------
## (Intercept)     32        7,383.100       4,694      102.952     571.160  
## disp         7,383.100  2,179,627.000 1,291,364.000 27,091.490 128,801.500
## hp             4,694    1,291,364.000    834,278    16,471.740 81,092.160 
## wt            102.952    27,091.490    16,471.740    360.901    1,828.095 
## qsec          571.160    128,801.500   81,092.160   1,828.095  10,293.480 
## --------------------------------------------------------------------------
library(stargazer)
options(scipen = 9999)
Sn <-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn, type = "text")
## 
## =============================
## 0.177   0     0     0     0  
## 0     0.001   0     0     0  
## 0       0   0.001   0     0  
## 0       0     0   0.053   0  
## 0       0     0     0   0.010
## -----------------------------

XᵗX Normalizada.

# Calculando matriz normalizada

XX_norm <- (Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm, type = "text", digits = 4)
## 
## ==================================
## 1      0.8840 0.9085 0.9580 0.9952
## 0.8840   1    0.9576 0.9659 0.8599
## 0.9085 0.9576   1    0.9493 0.8751
## 0.9580 0.9659 0.9493   1    0.9485
## 0.9952 0.8599 0.8751 0.9485   1   
## ----------------------------------

Autovalores de XᵗX

library(stargazer)

#Autovalores

lambdas <-eigen(XX_norm, symmetric = TRUE)
stargazer(lambdas$values, type = "text")
## 
## =============================
## 4.721 0.217 0.050 0.010 0.001
## -----------------------------

Calculo de K(X) = √(λmax  λmin)

K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)
## [1] 57.48052

INTERPRETACIÓN: La multicolinealidad es severa por que (X) ES DE 57.48 lo cual es mayor o igual a 20 y se considera un problema para el modelo ya que hay correlación entre las variables.

Cálculo del Indice de Condición usando librería “mctest”

library(mctest)
X_mat<-model.matrix(modelo_estimado)
mctest(mod = modelo_estimado)
## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0247         0
## Farrar Chi-Square:       106.7504         1
## Red Indicator:             0.6542         1
## Sum of Lambda Inverse:    23.2023         1
## Theil's Method:            0.7121         1
## Condition Number:         57.4805         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

INTERPRETACIÓN: El modelo evidencia multicolinealidad severa, lo que afecta la estabilidad e interpretación de los estimadores por lo que se considera un problema para seguir analizando el modelo.

Cálculo del Indice de Condición usando librería “olsrr”

library(olsrr)
## Warning: package 'olsrr' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'olsrr'
## The following object is masked from 'package:datasets':
## 
##     rivers
ols_eigen_cindex(model = modelo_estimado)
##    Eigenvalue Condition Index   intercept        disp          hp           wt
## 1 4.721487187        1.000000 0.000123237 0.001132468 0.001413094 0.0005253393
## 2 0.216562203        4.669260 0.002617424 0.036811051 0.027751289 0.0002096014
## 3 0.050416837        9.677242 0.001656551 0.120881424 0.392366164 0.0377028008
## 4 0.010104757       21.616057 0.025805998 0.777260487 0.059594623 0.7017528428
## 5 0.001429017       57.480524 0.969796790 0.063914571 0.518874831 0.2598094157
##           qsec
## 1 0.0001277169
## 2 0.0046789491
## 3 0.0001952599
## 4 0.0024577686
## 5 0.9925403056

Prueba de Farrar-Glaubar.

# Calculo Manual

library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "text")
## 
## =============================================
##                    disp    hp     wt    qsec 
## ---------------------------------------------
## Mazda RX4         -0.571 -0.535 -0.610 -0.777
## Mazda RX4 Wag     -0.571 -0.535 -0.350 -0.464
## Datsun 710        -0.990 -0.783 -0.917 0.426 
## Hornet 4 Drive    0.220  -0.535 -0.002 0.890 
## Hornet Sportabout 1.043  0.413  0.228  -0.464
## Valiant           -0.046 -0.608 0.248  1.327 
## ---------------------------------------------

Calcular la matriz R

library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))

stargazer(R,type = "text",digits = 4)
## 
## ====================================
##       disp     hp      wt     qsec  
## ------------------------------------
## disp    1    0.7909  0.8880  -0.4337
## hp   0.7909     1    0.6587  -0.7082
## wt   0.8880  0.6587     1    -0.1747
## qsec -0.4337 -0.7082 -0.1747    1   
## ------------------------------------

Calcular |R|

determinante_R<-det(R)
print(determinante_R)
## [1] 0.02466606

Aplicando la prueba de Farrer Glaubar (Bartlett)

m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)
## [1] 106.7504

Valor Critico.

gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)
## [1] 12.59159

Regla de desición:

INTERPRETACIÓN: Como χ2FG≥ es 106.7504 >= lo cual el valor critico es de 12.59159 se rechaza Ho, por lo tanto hay evidencia de colinealidad en los regresores

##Cálculo de FG usando “mctest”

library(mctest)
mctest::omcdiag(mod = modelo_estimado)
## 
## Call:
## mctest::omcdiag(mod = modelo_estimado)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0247         0
## Farrar Chi-Square:       106.7504         1
## Red Indicator:             0.6542         1
## Sum of Lambda Inverse:    23.2023         1
## Theil's Method:            0.7121         1
## Condition Number:         57.4805         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

INTERPRETACIÓN: Otra vez tenemos que el Chi - Square es de 57.4805 superior al valor critico de 7.8 encontrado por lo tanto se considera como bien lo que dice el diagnostico colinealidad en el modelo.

Cálculo de FG usando la “psych”

library(psych)
## Warning: package 'psych' was built under R version 4.5.3
FG_test<-cortest.bartlett(X_mat[,-1])
## R was not square, finding R from data
print(FG_test)
## $chisq
## [1] 106.7504
## 
## $p.value
## [1] 0.000000000000000000009757936
## 
## $df
## [1] 6

Se confirma que el estadistico FG es de 106.7504

Factores Inflacionarios de la Varianza (FIV)

Referencia entre R2j

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
R.cuadrado.regresores<-c(0,0.5,.8,.9)
R<-as.data.frame(R.cuadrado.regresores) %>% mutate(VIF=1/(1-R.cuadrado.regresores))
print(R)
##   R.cuadrado.regresores VIF
## 1                   0.0   1
## 2                   0.5   2
## 3                   0.8   5
## 4                   0.9  10

Calculo Manual

print(R)
##   R.cuadrado.regresores VIF
## 1                   0.0   1
## 2                   0.5   2
## 3                   0.8   5
## 4                   0.9  10

Inversa de la matriz de correlación R−1

modelo_estimado <-mtcars[, c("disp","hp","wt","qsec")]

# Creando la matriz correlación
matriz_cor <-cor(modelo_estimado)

# Calculando la inversa
inversa_R <-solve(matriz_cor)
print(inversa_R)
##           disp        hp        wt      qsec
## disp  7.985439 -1.398162 -5.918456  1.439009
## hp   -1.398162  5.166758 -1.679954  2.759325
## wt   -5.918456 -1.679954  6.916942 -2.548105
## qsec  1.439009  2.759325 -2.548105  3.133119

VIF’s para el modelo estimado

VIFs<-diag(inversa_R)
print(VIFs)
##     disp       hp       wt     qsec 
## 7.985439 5.166758 6.916942 3.133119

Cálculo de los VIF’s usando “performance”

library(performance)
## Warning: package 'performance' was built under R version 4.5.3
VIFs<-multicollinearity(x = modelo_estimado,verbose = FALSE)
VIFs
## NULL
library(performance)
library(see)
## Warning: package 'see' was built under R version 4.5.3
modelo_estimado <- lm( qsec ~ disp + hp + wt, data = mtcars)
colinealidad <- check_collinearity(modelo_estimado)
plot(colinealidad)

Cálculo de los VIF’s usando “car”

library(car)
## Cargando paquete requerido: carData
## 
## Adjuntando el paquete: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:psych':
## 
##     logit
VIFs_car<-vif(modelo_estimado)
print(VIFs_car)
##     disp       hp       wt 
## 7.324517 2.736633 4.844618

Cálculo de los VIF’s usando “mctest”

library(mctest)
mc.plot(mod = modelo_estimado,vif = 2)

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 4.5.3
## Cargando paquete requerido: MASS
## 
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:olsrr':
## 
##     cement
## Cargando paquete requerido: survival
fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm")
plot(fit_normal)

Prueba de Normalidad de Jarque Bera

Usando tseries

library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
salida_JB<-jarque.bera.test(modelo_estimado$residuals)
salida_JB
## 
##  Jarque Bera Test
## 
## data:  modelo_estimado$residuals
## X-squared = 11.533, df = 2, p-value = 0.003131
library(fastGraph)
## Warning: package 'fastGraph' was built under R version 4.5.3
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)))

cat("JB =", salida_JB$statistic, "\n")
## JB = 11.53256
cat("p-value =", salida_JB$p.value, "\n")
## p-value = 0.003131386

INTERPRETACIÓN: Con un p-value 0.1905, mayor nivel de significancia del 5%, no se rechaza la hipótesis nula de normalidad. Por lo tanto los datos analizados siguen una distribución normal.

Prueba de Kolmogorov Smirnov -Lilliefors

Cálculo Manual

library(dplyr)
library(gt)
## Warning: package 'gt' was built under R version 4.5.3
library(gtExtras)
## Warning: package 'gtExtras' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'gtExtras'
## The following object is masked from 'package:MASS':
## 
##     select
# Residuos del Modelo

residuos<-residuals(modelo_estimado)

# Tabla KS

residuos %>%
  as_tibble() %>%
  mutate(posicion = row_number()) %>%
  arrange(value) %>%
  mutate(dist1 = row_number()/n()) %>%
  mutate(dist2 = (row_number()-1)/n()) %>%
  mutate(zi = as.vector(scale(value, center = TRUE))) %>%
  mutate(pi = pnorm(zi, lower.tail = TRUE)) %>%
  mutate(dif1 = abs(dist1 - pi)) %>%
  mutate(dif2 = abs(dist2 - pi)) %>%
  rename(residuales = value) -> tabla_KS

# Formato

tabla_KS %>%
  gt() %>%
  tab_header(title = "Tabla para calcular el Estadístico KS") %>%
  tab_source_note(source_note = "Fuente: Elaboración propia") %>%
  
# Resaltar máximo dif1
  
  tab_style(
    style = list(
      cell_fill(color = "#A569BD"),
      cell_text(style = "italic")
    ),
    locations = cells_body(
      columns = dif1,
      rows = dif1 == max(tabla_KS$dif1)
    )
  ) %>%

# Resaltar máximo dif2
  tab_style(
    style = list(
      cell_fill(color = "#3498DB"),
      cell_text(style = "italic")
    ),
    locations = cells_body(
      columns = dif2,
      rows = dif2 == max(tabla_KS$dif2)
    )
  )
Tabla para calcular el Estadístico KS
residuales posicion dist1 dist2 zi pi dif1 dif2
-1.81209054 1 0.03125 0.00000 -1.79497163 0.03632911 0.0050791104 0.036329110
-1.63083769 2 0.06250 0.03125 -1.61543108 0.05310869 0.0093913132 0.021858687
-1.60224000 30 0.09375 0.06250 -1.58710355 0.05624460 0.0375054035 0.006255404
-1.55816416 27 0.12500 0.09375 -1.54344410 0.06136153 0.0636384698 0.032388470
-0.84130055 10 0.15625 0.12500 -0.83335274 0.20232291 0.0460729087 0.077322909
-0.77442872 22 0.18750 0.15625 -0.76711265 0.22150727 0.0340072702 0.065257270
-0.65218051 12 0.21875 0.18750 -0.64601933 0.25913342 0.0403834180 0.071633418
-0.31723832 24 0.25000 0.21875 -0.31424135 0.37666887 0.1266688694 0.157918869
-0.31088794 23 0.28125 0.25000 -0.30795097 0.37905982 0.0978098218 0.129059822
-0.30853956 8 0.31250 0.28125 -0.30562477 0.37994518 0.0674451804 0.098695180
-0.29090893 17 0.34375 0.31250 -0.28816070 0.38661186 0.0428618646 0.074111865
-0.24130055 11 0.37500 0.34375 -0.23902097 0.40554466 0.0305446593 0.061794659
-0.22010742 16 0.40625 0.37500 -0.21802806 0.41370363 0.0074536260 0.038703626
-0.19094800 32 0.43750 0.40625 -0.18914411 0.42498994 0.0125100631 0.018739937
-0.14891664 19 0.46875 0.43750 -0.14750982 0.44136482 0.0273851835 0.003864816
-0.08879107 28 0.50000 0.46875 -0.08795225 0.46495731 0.0350426868 0.003792687
0.03969579 25 0.53125 0.50000 0.03932078 0.51568268 0.0155673180 0.015682682
0.04826170 15 0.56250 0.53125 0.04780577 0.51906448 0.0434355210 0.012185521
0.04894094 3 0.59375 0.56250 0.04847860 0.51933259 0.0744174111 0.043167411
0.05281569 13 0.62500 0.59375 0.05231674 0.52086184 0.1041381584 0.072888158
0.09899357 26 0.65625 0.62500 0.09805837 0.53905703 0.1171929704 0.085942970
0.21063854 29 0.68750 0.65625 0.20864863 0.58263873 0.1048612716 0.073611272
0.27340698 18 0.71875 0.68750 0.27082409 0.60673684 0.1120131640 0.080763164
0.34635461 5 0.75000 0.71875 0.34308258 0.63423183 0.1157681679 0.084518168
0.37855154 14 0.78125 0.75000 0.37497535 0.64616060 0.1350894006 0.103839401
0.58000827 7 0.81250 0.78125 0.57452891 0.71719503 0.0953049723 0.064054972
0.93312402 4 0.84375 0.81250 0.92430875 0.82233720 0.0214127952 0.009837205
1.01511828 31 0.87500 0.84375 1.00552840 0.84267876 0.0323212401 0.001071240
1.01593588 6 0.90625 0.87500 1.00633828 0.84287356 0.0633764366 0.032126437
1.17225448 20 0.93750 0.90625 1.16118013 0.87721567 0.0602843265 0.029034327
1.40551509 21 0.96875 0.93750 1.39223712 0.91807470 0.0506753028 0.019425303
3.36926520 9 1.00000 0.96875 3.33743559 0.99957722 0.0004227766 0.030827223
Fuente: Elaboración propia

Calculo de estadístico D

D<-max(max(tabla_KS$dif1),max(tabla_KS$dif2))
print(D)
## [1] 0.1579189
#Calculando valor lilliefors
library(nortest)

n <- 88
alfa <- 0.05

set.seed(123)
iteraciones <- 100000
distribucion_d <- numeric(iteraciones)

for(i in 1:iteraciones) {
  x <- rnorm(n)
  distribucion_d[i] <- lillie.test(x)$statistic
}

valor_critico <- quantile(distribucion_d, 1 - alfa)

cat("Valor crítico Lilliefors:", valor_critico, "\n")
## Valor crítico Lilliefors: 0.09460697
library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals)
prueba_KS
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_estimado$residuals
## D = 0.15792, p-value = 0.04119

INTERPRETACIÓN : con un p-value de 0.3446, mayor al nivel de significancia del 5%, no se rechaza la hipótesis nula de normalidad. Por lo tanto, los residuos del modelo siguen una distribución normal.

Prueba de Shapiro-Wilk SW

# Calculo Manual

library(dplyr)
library(gt)

# Residuales del modelo
residuos <- modelo_estimado$residuals

# Tamaño de muestra
n <- length(residuos)

tabla_SW <- residuos %>%  
  as_tibble() %>%
  rename(residuales = value) %>%
  arrange(residuales) %>%
  mutate(
    pi = (row_number() - 0.375) / (n + 0.25),
    mi = qnorm(pi),
    ai = 0
  )

m <- sum(tabla_SW$mi^2)
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 <- tabla_SW %>%
  mutate(
    ai_ui = ai * residuales,
    ui2 = residuales^2
  )

tabla_SW %>%
  gt() %>%
  tab_header(title = "Tabla para calcular el Estadístico W") %>%
  tab_source_note(source_note = "Fuente: Elaboración propia")
Tabla para calcular el Estadístico W
residuales pi mi ai ai_ui ui2
-1.81209054 0.01937984 -2.06672907 -0.407971522 0.7392813330 3.283672108
-1.63083769 0.05038760 -1.64110706 -0.296267293 0.4831638668 2.659631565
-1.60224000 0.08139535 -1.39574699 -0.248692215 0.3984646131 2.567173005
-1.55816416 0.11240310 -1.21384688 -0.216281511 0.3370020977 2.427875536
-0.84130055 0.14341085 -1.06511983 -0.189781537 0.1596633124 0.707786620
-0.77442872 0.17441860 -0.93684711 -0.166926086 0.1292723547 0.599739840
-0.65218051 0.20542636 -0.82239396 -0.146532987 0.0955659589 0.425339422
-0.31723832 0.23643411 -0.71782010 -0.127900165 0.0405748331 0.100640150
-0.31088794 0.26744186 -0.62056827 -0.110571972 0.0343754927 0.096651311
-0.30853956 0.29844961 -0.52886482 -0.094232382 0.0290744174 0.095196657
-0.29090893 0.32945736 -0.44141204 -0.078650169 0.0228800368 0.084628007
-0.24130055 0.36046512 -0.35721583 -0.063648210 0.0153583482 0.058225957
-0.22010742 0.39147287 -0.27548228 -0.049085041 0.0108039819 0.048447277
-0.19094800 0.42248062 -0.19555148 -0.034843085 0.0066532174 0.036461139
-0.14891664 0.45348837 -0.11685275 -0.020820656 0.0031005421 0.022176165
-0.08879107 0.48449612 -0.03887224 -0.006926201 0.0006149847 0.007883853
0.03969579 0.51550388 0.03887224 0.006926201 0.0002749410 0.001575756
0.04826170 0.54651163 0.11685275 0.020820656 0.0010048402 0.002329191
0.04894094 0.57751938 0.19555148 0.034843085 0.0017052534 0.002395216
0.05281569 0.60852713 0.27548228 0.049085041 0.0025924603 0.002789497
0.09899357 0.63953488 0.35721583 0.063648210 0.0063007636 0.009799727
0.21063854 0.67054264 0.44141204 0.078650169 0.0165667572 0.044368597
0.27340698 0.70155039 0.52886482 0.094232382 0.0257637914 0.074751379
0.34635461 0.73255814 0.62056827 0.110571972 0.0382971123 0.119961515
0.37855154 0.76356589 0.71782010 0.127900165 0.0484168048 0.143301270
0.58000827 0.79457364 0.82239396 0.146532987 0.0849903452 0.336409599
0.93312402 0.82558140 0.93684711 0.166926086 0.1557627394 0.870720428
1.01511828 0.85658915 1.06511983 0.189781537 0.1926507069 1.030465113
1.01593588 0.88759690 1.21384688 0.216281511 0.2197281460 1.032125702
1.17225448 0.91860465 1.39574699 0.248692215 0.2915305623 1.374180561
1.40551509 0.94961240 1.64110706 0.296267293 0.4164081510 1.975472671
3.36926520 0.98062016 2.06672907 0.407971522 1.3745642513 11.351948002
Fuente: Elaboración propia

Calculo del estadístico W

W<-(sum(tabla_SW$ai_ui)^2)/sum(tabla_SW$ui2)
print(W)
## [1] 0.9169524

Calculo 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)
Wn<-(log(1-W)-mu)/sigma
print(Wn)
## [1] 2.115314
p.value<-pnorm(Wn,lower.tail = FALSE)
print(p.value)
## [1] 0.01720159
library(fastGraph)
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE)

INTERPRETACIÓN: El modelo cumple con el supuesto de normalidad, lo que significa que los resultados son confiables.

Usando la libreria Stats

salida_SW<-shapiro.test(modelo_estimado$residuals)
print(salida_SW)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_estimado$residuals
## W = 0.91695, p-value = 0.0172