Guía de Laboratorio: Diagnóstico de Supuestos en Modelos de Regresión Lineal Múltiple

# 1. Cargar el dataset (ya viene incluido en R)
data(mtcars)

# 2. Ver las primeras filas del dataset
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Instrucciones: Carga el dataset mtcars en R. El modelo a analizar busca explicar el consumo de combustible (mpg) en

función de: disp (desplazamiento), hp (potencia), wt (peso) y qsec (tiempo en cuarto de milla).

 #Ajustar el modelo de regresión lineal
library(stargazer)

## Warning: package 'stargazer' was built under R version 4.5.2

## 
## 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)

#Ver el resumen del modelo
stargazer(modelo_estimado,type = "text",title = "modelo_estimado")

## 
## modelo
## ===============================================
##                         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

2. Diagnóstico de MulticolinealidadPrueba de Farrar-Glaubar

# - Cálculo “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 = 999)
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
## -----------------------------

XtX normalizada:

library(stargazer)
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 XtX Normalizada

library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")

## 
## =============================
## 4.721 0.217 0.050 0.010 0.001
## -----------------------------

Cálculo de κ(x)=√λmaxλmin3

K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)

## [1] 57.48052

Como κ(x)>30 se considera que la multicolinealidad es severa. # Cálculo del Indice de Condición usando librería “mctest”

library(mctest)

## Warning: package 'mctest' was built under R version 4.5.2

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

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

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 de |R| - Cálculo “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))
#También se puede calcular R a través de cor(X_mat[,-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) Estadistico χ2FG

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

Como χ2FG≥V.C.se rechaza H0, 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

- 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

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)
as.data.frame(R.cuadrado.regresores) %>% mutate(VIF=1/(1-R.cuadrado.regresores))

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

Cálculo manual:

print(R)

##            disp         hp         wt       qsec
## disp  1.0000000  0.7909486  0.8879799 -0.4336979
## hp    0.7909486  1.0000000  0.6587479 -0.7082234
## wt    0.8879799  0.6587479  1.0000000 -0.1747159
## qsec -0.4336979 -0.7082234 -0.1747159  1.0000000

#Inversa de la matriz de correlación R−1
inversa_R<-solve(R)
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

## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##  qsec 3.13 [2.10,  5.15]     1.77      0.32     [0.19, 0.48]
## 
## Moderate Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##  disp 7.99 [4.92, 13.44]     2.83      0.13     [0.07, 0.20]
##    hp 5.17 [3.28,  8.62]     2.27      0.19     [0.12, 0.30]
##    wt 6.92 [4.30, 11.61]     2.63      0.14     [0.09, 0.23]

plot(VIFs)

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

library(car)

## Warning: package 'car' was built under R version 4.5.3

## Cargando paquete requerido: carData

## Warning: package 'carData' was built under R version 4.5.2

## 
## 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     qsec 
## 7.985439 5.166758 6.916942 3.133119

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

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

Preguntas de reflexión: 1. Si el valor de VIF es mayor a 5 o 10, ¿qué implica esto para la interpretación de los coeficientes? 2. Observando los resultados: ¿Existe alguna variable en este modelo que parezca estar “robándole” significancia a las demás? Justifica basándote en su correlación técnica (ej. desplazamiento vs. potencia). 3. ¿Qué indica el Indice de condición en este modelo?

1. Cuando el VIF es mayor a 5 o 10, implica la presencia de multicolinealidad en el modelo, es decir, que algunas variables explicativas están fuertemente relacionadas entre sí. Esto provoca que los coeficientes estimados sean inestables y poco confiables, dificultando su interpretación individual, ya que el efecto de una variable puede estar mezclado con el de otra.

2. En el modelo analizado, sí existe evidencia de que algunas variables están “robándose” significancia, especialmente variables como desplazamiento (disp) y potencia (hp), que suelen estar altamente correlacionadas. Debido a esta relación, una de ellas puede absorber el efecto explicativo de la otra, haciendo que una variable parezca no significativa aunque en realidad sí influya sobre la variable dependiente.

3. El índice de condición obtenido, al ser mayor a 30, indica que existe multicolinealidad severa en el modelo. Esto significa que hay una fuerte dependencia lineal entre las variables explicativas, lo que afecta la estabilidad de los coeficientes y la calidad de las inferencias realizadas a partir del modelo.

3. Diagnóstico de Normalidad de los Residuos

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)

## 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

## Warning: package 'survival' was built under R version 4.5.2

fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm")
plot(fit_normal)

summary(fit_normal)

## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##                        estimate Std. Error
## mean -0.00000000000000003122502  0.4257752
## sd    2.40854809054366647558254  0.3010683
## Loglikelihood:  -73.5348   AIC:  151.0696   BIC:  154.0011 
## Correlation matrix:
##                     mean                  sd
## mean  1.0000000000000000 -0.0000000004554131
## sd   -0.0000000004554131  1.0000000000000000

1. Prueba de Normalidad de Jarque Bera Usando tseries

library(tseries)

## Warning: package 'tseries' was built under R version 4.5.2

## 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 = 3.316, df = 2, p-value = 0.1905

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)))

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

cat("JB =", salida_JB$statistic, "\n")

## JB = 3.315957

cat("p-value =", salida_JB$p.value, "\n")

## p-value = 0.1905237

2. 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 de TU 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
-3.8664153	21	0.03125	0.00000	-1.58000709	0.05705262	0.025802622	0.057052622
-3.7219399	23	0.06250	0.03125	-1.52096733	0.06413402	0.001634018	0.032884018
-2.9249910	6	0.09375	0.06250	-1.19529489	0.11598592	0.022235918	0.053485918
-2.8335055	22	0.12500	0.09375	-1.15790946	0.12345049	0.001549506	0.029700494
-2.5153146	3	0.15625	0.12500	-1.02788103	0.15200290	0.004247104	0.027002896
-2.1098370	11	0.18750	0.15625	-0.86218297	0.19429342	0.006793419	0.038043419
-1.8775253	14	0.21875	0.18750	-0.76724902	0.22146674	0.002716736	0.033966736
-1.5839621	1	0.25000	0.21875	-0.64728469	0.25872385	0.008723854	0.039973854
-1.5812652	7	0.28125	0.25000	-0.64618256	0.25908057	0.022169435	0.009080565
-1.5256810	32	0.31250	0.28125	-0.62346815	0.26648846	0.046011536	0.014761536
-1.0761488	24	0.34375	0.31250	-0.43976723	0.33005285	0.013697147	0.017552853
-1.0743171	9	0.37500	0.34375	-0.43901872	0.33032399	0.044676013	0.013426013
-0.8170837	29	0.40625	0.37500	-0.33390049	0.36922732	0.037022675	0.005772675
-0.7133657	2	0.43750	0.40625	-0.29151623	0.38532827	0.052171732	0.020921732
-0.4169002	30	0.46875	0.43750	-0.17036590	0.43236119	0.036388805	0.005138805
-0.3833408	10	0.50000	0.46875	-0.15665188	0.43775961	0.062240392	0.030990392
-0.3742943	26	0.53125	0.50000	-0.15295502	0.43921687	0.092033127	0.060783127
-0.3244422	4	0.56250	0.53125	-0.13258300	0.44726159	0.115238410	0.083988410
0.0522694	15	0.59375	0.56250	0.02135984	0.50852070	0.085229304	0.053979304
0.2096827	13	0.62500	0.59375	0.08568663	0.53414223	0.090857767	0.059607767
0.2708056	5	0.65625	0.62500	0.11066445	0.54405878	0.112191219	0.080941219
0.8242600	27	0.68750	0.65625	0.33683309	0.63187864	0.055621362	0.024371362
0.9856164	12	0.71875	0.68750	0.40277125	0.65644174	0.062308257	0.031058257
1.1599816	16	0.75000	0.71875	0.47402542	0.68225911	0.067740891	0.036490891
1.2050398	19	0.78125	0.75000	0.49243841	0.68879527	0.092454727	0.061204727
1.6307713	31	0.81250	0.78125	0.66641320	0.74742649	0.065073512	0.033823512
1.6563961	8	0.84375	0.81250	0.67688474	0.75076046	0.092989544	0.061739544
2.5145182	25	0.87500	0.84375	1.02755556	0.84792053	0.027079468	0.004170532
2.7033586	28	0.90625	0.87500	1.10472502	0.86536062	0.040889380	0.009639380
5.2230317	20	0.93750	0.90625	2.13438716	0.98359445	0.046094449	0.077344449
5.6377518	18	0.96875	0.93750	2.30386212	0.98938481	0.020634808	0.051884808
5.6468467	17	1.00000	0.96875	2.30757875	0.98948871	0.010511289	0.020738711
Fuente: Elaboración propia

Cálculo del estadístico D:

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

## [1] 0.1152384

library(nortest) # calculando valorLilliefors

## Warning: package 'nortest' was built under R version 4.5.2

n <- 32 # OJO CON OBSE...
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.1538328

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

Usando nortest

library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals)
prueba_KS

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

p.value<-prueba_KS$p.value

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

3. Prueba de Shapiro - Wilk Cálculo 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
-3.8664153	0.01937984	-2.06672907	-0.407971522	1.577387330	14.94916721
-3.7219399	0.05038760	-1.64110706	-0.296267293	1.102689045	13.85283628
-2.9249910	0.08139535	-1.39574699	-0.248692215	0.727422479	8.55557209
-2.8335055	0.11240310	-1.21384688	-0.216281511	0.612834858	8.02875360
-2.5153146	0.14341085	-1.06511983	-0.189781537	0.477360277	6.32680768
-2.1098370	0.17441860	-0.93684711	-0.166926086	0.352186833	4.45141219
-1.8775253	0.20542636	-0.82239396	-0.146532987	0.275119397	3.52510140
-1.5839621	0.23643411	-0.71782010	-0.127900165	0.202589020	2.50893608
-1.5812652	0.26744186	-0.62056827	-0.110571972	0.174843606	2.50039947
-1.5256810	0.29844961	-0.52886482	-0.094232382	0.143768554	2.32770247
-1.0761488	0.32945736	-0.44141204	-0.078650169	0.084639285	1.15809624
-1.0743171	0.36046512	-0.35721583	-0.063648210	0.068378362	1.15415730
-0.8170837	0.39147287	-0.27548228	-0.049085041	0.040106585	0.66762569
-0.7133657	0.42248062	-0.19555148	-0.034843085	0.024855860	0.50889056
-0.4169002	0.45348837	-0.11685275	-0.020820656	0.008680136	0.17380580
-0.3833408	0.48449612	-0.03887224	-0.006926201	0.002655095	0.14695018
-0.3742943	0.51550388	0.03887224	0.006926201	-0.002592437	0.14009622
-0.3244422	0.54651163	0.11685275	0.020820656	-0.006755099	0.10526272
0.0522694	0.57751938	0.19555148	0.034843085	0.001821227	0.00273209
0.2096827	0.60852713	0.27548228	0.049085041	0.010292282	0.04396682
0.2708056	0.63953488	0.35721583	0.063648210	0.017236289	0.07333566
0.8242600	0.67054264	0.44141204	0.078650169	0.064828188	0.67940453
0.9856164	0.70155039	0.52886482	0.094232382	0.092876982	0.97143970
1.1599816	0.73255814	0.62056827	0.110571972	0.128261451	1.34555727
1.2050398	0.76356589	0.71782010	0.127900165	0.154124786	1.45212086
1.6307713	0.79457364	0.82239396	0.146532987	0.238961790	2.65941502
1.6563961	0.82558140	0.93684711	0.166926086	0.276495710	2.74364789
2.5145182	0.85658915	1.06511983	0.189781537	0.477209125	6.32280165
2.7033586	0.88759690	1.21384688	0.216281511	0.584686474	7.30814753
5.2230317	0.91860465	1.39574699	0.248692215	1.298927330	27.28006053
5.6377518	0.94961240	1.64110706	0.296267293	1.670281449	31.78424483
5.6468467	0.98062016	2.06672907	0.407971522	2.303752631	31.88687740
Fuente: Elaboración propia

Cálculo del Estadistico W

W<-(sum(tabla_SW$ai_ui)^2)/sum(tabla_SW$ui2)
print(W)

## [1] 0.9366138

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)
Wn<-(log(1-W)-mu)/sigma
print(Wn)

## [1] 1.554447

p.value<-pnorm(Wn,lower.tail = FALSE)
print(p.value)

## [1] 0.06003893

library(fastGraph)
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE)

En este caso, dado que 0.645292 > 0.05, no se rechaza la Hipótesis Nula: ϵ ∼ N(0, σ²), 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)
print(salida_SW)

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

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

## [1] 1.554447

En general, todas las pruebas de normalidad aplicadas (Lilliefors, Jarque-Bera y Shapiro-Wilk) son coincidentes, ya que en todos los casos el p-value es mayor a 0.05, por lo que no se rechaza la Hipótesis Nula: ϵ ∼ N(0, σ²), concluyéndose que los residuos siguen una distribución normal.

4. Análisis Crítico Final (Reporte)

Redacta un párrafo (máximo 150 palabras) respondiendo lo siguiente:

 ¿Es este modelo confiable para realizar inferencia estadística?

 Si detectaste problemas de multicolinealidad o falta de normalidad, ¿qué medidas correctivas propondrías?

(Pista: ¿Sería correcto eliminar una variable, transformar los datos o buscar una muestra más grande?)

El modelo presenta resultados mixtos para la inferencia estadística. Por un lado, las pruebas de normalidad (Lilliefors, Jarque-Bera y Shapiro-Wilk) indican que los residuos siguen una distribución normal, lo cual es favorable. Sin embargo, se detecta multicolinealidad severa (VIF elevados e índice de condición mayor a 30), lo que afecta la estabilidad e interpretación de los coeficientes. Por ello, el modelo no es completamente confiable para inferencias individuales sobre las variables explicativas. Como medidas correctivas, se recomienda eliminar o combinar variables altamente correlacionadas (como desplazamiento y potencia), aplicar transformaciones o utilizar métodos alternativos como regresión ridge osea una forma de ajustar el modelo para que no se “confunda”. También podría considerarse ampliar la muestra para mejorar la robustez del modelo.