Enunciado: sea el conjunto de datos, tomados en 24 meses correspondientes a los gastos de comercialización (C) de una empresa, el nivel de ventas (V) su coste de personal (P) y los costes de materias primas (M); se trata de estimar el nivel de ventas a partir de las variables restantes.

Cargar dataset

library(readxl)
ventas_empresa <- read_excel("D:/harold/Descargas/ventas_empresa.xlsx")

Estimar modelo de regresión lineal

modelo_ventas <- lm( V ~ C + P + M, data = ventas_empresa)
summary(modelo_ventas)
## 
## Call:
## lm(formula = V ~ C + P + M, data = ventas_empresa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.279  -6.966   1.555   6.721  14.719 
## 
## Coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept) 107.4435    18.0575   5.950 0.00000808 ***
## C             0.9226     0.2227   4.142   0.000505 ***
## P             0.9502     0.1558   6.097 0.00000586 ***
## M             1.2978     0.4307   3.013   0.006872 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.506 on 20 degrees of freedom
## Multiple R-squared:  0.9798, Adjusted R-squared:  0.9768 
## F-statistic: 323.6 on 3 and 20 DF,  p-value: < 0.00000000000000022

1. Verificar Supuesto de Normalidad de los residuos

1.1. Prueba de Normalidad Jarque-Bera (JB)

Utilizando la función jarque.bera.test de la librería “tseries”.

library(tseries)
jb_test <- jarque.bera.test(modelo_ventas$residuals)
jb_test
## 
##  Jarque Bera Test
## 
## data:  modelo_ventas$residuals
## X-squared = 1.4004, df = 2, p-value = 0.4965

Resultado:
Como \(p=0.4965 > \alpha=0.05\), no se rechaza la hipótesis nula. No hay evidencia suficiente para decir que la distribución se desvía de la normalidad.

1.1.1. Gráfica del estadístico de la prueba JB

El estadístico JB sigue asintóticamente una distribución Chi-cuadrado (\(\chi^2\)) con dos grados de libertad (uno por la asimetría y otro por la curtosis), gráficamente se puede observar de la siguiente manera:

library(fastGraph)
alpha <- 0.05
JB <- jb_test$statistic
gl <- jb_test$parameter
VC <- qchisq(1-alpha, gl, lower.tail = TRUE)
shadeDist(JB, ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE, xmin = 0,
          sub=paste("Valor critico:", round(VC,2),"","JB", round(JB,2)))

El resultado converge con la prueba de JB hecha a través de función jarque.bera.test, de la librerua “tseries”, como el valor de \(\text{JB}=1.4 \leq \text{V.C.=5.99}\), no se rechaza la hipótesis nula. No hay evidencia suficiente para decir que la distribución se desvía de la normalidad.

1.2. Prueba de Kolgomoróv-Smirnov (KS)

Utilizando la función lillie.test para realizar la prueba KS con corrección de Lilliefors, de la librería “nortest”.

library(nortest)
ks_test <- lillie.test(modelo_ventas$residuals)
ks_test
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_ventas$residuals
## D = 0.13659, p-value = 0.2935

Resultado:
Como \(p=0.2935>\alpha=0.05\), no se rechaza la hipotesis nula. No hay evidencia suficiente para decir que la distribución se desvía de la normalidad.

1.3. Prueba de Shapiro-Wilk (SW)

Utilizando la función shapiro.test de la librería precargada “stats”.

sw_test <- shapiro.test(modelo_ventas$residuals)
sw_test
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_ventas$residuals
## W = 0.95315, p-value = 0.3166

Resultado:
Como \(p=0.3166 > \alpha=0.05\), no se rechaza la hipótesis nula. No hay evidencia suficiente para decir que la distribución se desvía de la normalidad.

1.3.1. Gráfica del estadístico de la prueba SW

Para ilustrar el resultado de la prueba de Shapiro-Wilk, se puede graficar el estadístico W, que sigue asintóticamente una distribución normal estándar (Z), de la siguiente manera:
A. Tabla para calcular el estadístico W

library(dplyr)
library(gt)
residuos <- modelo_ventas$residual
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(ventas_empresa)
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
-17.27950969 0.02577320 -1.94690278 -0.44751326 7.7328097356 298.5814551698
-15.50384667 0.06701031 -1.49843365 -0.31302439 4.8530821064 240.3692617022
-13.84401165 0.10824742 -1.23590240 -0.25502179 3.5305246473 191.6566585921
-9.91393148 0.14948454 -1.03864671 -0.21431914 2.1247453168 98.2860374647
-8.43534410 0.19072165 -0.87524006 -0.18060106 1.5234321207 71.1550300366
-7.90808316 0.23195876 -0.73241136 -0.15112913 1.1951417479 62.5377792373
-6.65267950 0.27319588 -0.60317579 -0.12446207 0.8280062569 44.2581445625
-3.33614422 0.31443299 -0.48332361 -0.09973122 0.3327177289 11.1298582608
-3.32226362 0.35567010 -0.37005675 -0.07635921 0.2536854377 11.0374355619
-2.43553165 0.39690722 -0.26136061 -0.05393035 0.1313490763 5.9318144083
-0.40251845 0.43814433 -0.15567569 -0.03212284 0.0129300352 0.1620211022
0.01280369 0.47938144 -0.05170609 -0.01066927 -0.0001366061 0.0001639345
3.09667909 0.52061856 0.05170609 0.01066927 0.0330393120 9.5894213810
4.04562361 0.56185567 0.15567569 0.03212284 0.1299569147 16.3670704012
4.98802235 0.60309278 0.26136061 0.05393035 0.2690057954 24.8803669567
5.88324503 0.64432990 0.37005675 0.07635921 0.4492399641 34.6125720380
6.28004614 0.68556701 0.48332361 0.09973122 0.6263166548 39.4389794795
6.69809559 0.72680412 0.60317579 0.12446207 0.8336588371 44.8644845291
6.78865323 0.76804124 0.73241136 0.15112913 1.0259632737 46.0858126562
7.37251119 0.80927835 0.87524006 0.18060106 1.3314833668 54.3539212565
8.70403920 0.85051546 1.03864671 0.21431914 1.8654422374 75.7602984046
9.77138853 0.89175258 1.23590240 0.25502179 2.4919170034 95.4800337177
10.67367778 0.93298969 1.49843365 0.31302439 3.3411214463 113.9273972668
14.71907878 0.97422680 1.94690278 0.44751326 6.5869829488 216.6512801774
Fuente: Elaboración propia

B. Cálculo del estadístico W

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

C. Cálculo del \(W_n\) y su valor \(p\)

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
Wn
## [1] 0.4773519
p.value <- pnorm(Wn, lower.tail = FALSE)
p.value
## [1] 0.3165558

D. Gráfica del estadístico \(Wn\), que sigue asintóticamente una distribución normal estándar (Z)

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

Resultado:
El resultado es convergente, \(p=0.3166 > \alpha=0.05\), no se rechaza la hipótesis nula. No hay evidencia suficiente para decir que la distribución se desvía de la normalidad. El estadístico \(W\) es cercano a uno (1) [\(W=0.9531481\)].

2. Verificar el Supuesto de No Colinealidad entre las variables explicativas

Más bien, verificar que la relación entre las variables explicativas no resulte problemática para su interpretación.

2.1. Indice de Condición (IC)

2.1.1. Indice de Condición con “mctest”

Utilizando la función mctest de la librería “mctest”.

Nota: se debe generar la matriz de diseño (\(\mathbf{X}\)) del modelo de regresión estimado.

library(mctest)
X <- model.matrix(modelo_ventas)
mctest(mod = modelo_ventas)
## 
## 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.0346         0
## Farrar Chi-Square:        71.2080         1
## Red Indicator:             0.8711         1
## Sum of Lambda Inverse:    20.9196         1
## Theil's Method:            0.5430         1
## Condition Number:         71.1635         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

Resultado:
Como se observa, Condition Numer es de 71.16 (\(\kappa(x)=71.16\)). Dado que el valor del Indice de Condición es mayor a 30 (\(\kappa(x)>30\)), se puede afirmar que hay evidencia de colinealidad entre las variables explicativas del modelo. Esto implica que las variables explicativas están altamente correlacionadas entre sí, lo que puede dificultar la interpretación de los coeficientes de regresión y afectar la estabilidad del modelo.

2.2. Indice de Condición con “olsrr”

Utilizando la función ols_eigen_cindex de la librería “olsrr”.

library(olsrr)
ols_eigen_cindex(modelo_ventas)
##     Eigenvalue Condition Index    intercept             C            P
## 1 3.9869237681         1.00000 0.0007213226 0.00009632167 0.0002714936
## 2 0.0095007154        20.48523 0.8775733849 0.00494748750 0.0877003286
## 3 0.0027882470        37.81406 0.1182991896 0.15940219484 0.8478049007
## 4 0.0007872695        71.16349 0.0034061029 0.83555399599 0.0642232771
##               M
## 1 0.00008216652
## 2 0.00754342318
## 3 0.06362314121
## 4 0.92875126909

Resultado:
El resultado es convergente al obtenido a través de la librería “mctest”, ya que el Indice de Condición es de 71.16, lo que indica, nuevamente, una colinealidad severa entre las variables explicativas del modelo.

2.3. Prueba de Farrar-Glaubar

2.3.1. Prueba de Farrar-Glaubar con “mctest”

utilizando la función omcdiag de la librería “mctest”.

library(mctest)
mctest::omcdiag(modelo_ventas)
## 
## Call:
## mctest::omcdiag(mod = modelo_ventas)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0346         0
## Farrar Chi-Square:        71.2080         1
## Red Indicator:             0.8711         1
## Sum of Lambda Inverse:    20.9196         1
## Theil's Method:            0.5430         1
## Condition Number:         71.1635         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

Valor crítico:
El valor crítico de la prueba de Farrar-Glaubar para un nivel de significancia del 5%, con 3 variables explicativas:

m <- ncol(X[,-1]) # Número de variables explicativas
n <- nrow(X[,-1]) # Número de observaciones
gl <- m*(m-1)/2 # Grados de libertad
VC <- qchisq(0.95, df = gl)
VC
## [1] 7.814728

Resultado:
Como esta prueba nos entrega solo el coeficiente de \(\chi^2\) de la prueba de Farrar-Glaubar, \(\chi^2_{FG}=71.2080>\text{V.C.}=7.814728\), se rechaza la hipótesis nula, la matriz de correlación no corresponde a una matriz identidad.

2.3.2. Prueba de Farrar-Glaubar con “psych”

Utilizando la función cortest.bartlett de la librería “psych”.

library(psych)
fg_test <- cortest.bartlett(X[,-1])
fg_test
## $chisq
## [1] 71.20805
## 
## $p.value
## [1] 0.000000000000002352605
## 
## $df
## [1] 3

Resultado:

El resultado es convergente al obtenido a través de la librería “mctest”, como se observa; además, se nos brinda el \(p\) valor de la prueba: como \(p=0.000000000000002352605<\alpha=0.05\), se rechaza la hipótesis nula, la matriz de correlación no corresponde a una matriz identidad.

2.3.3. Gráfica del estadístico de la prueba de Farrar-Glaubar

El estadístico de la prueba de Farrar-Glaubar sigue asintóticamente una distribución Chi-cuadrado (χ2\(\chi^2\)) con \(gl=m(m-1)/2\) grados de libertad, gráficamente se puede observar de la siguiente manera:

library(fastGraph)
FG <- fg_test$chisq
gl <- fg_test$df

shadeDist(FG, ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE, xmin = 0,
          main = "Resultados de la Prueba de Farrar-Glaubar (Chi-cuadrado)",
          sub=paste("Valor critico:", round(VC,2),"","FG", round(FG,2)))

2.4. Factor de Inflación de la Varianza (VIF)

2.4.1. VIF con “perfomance”

library(performance)
vif_performance <- multicollinearity(x = modelo_ventas, verbose = FALSE)
vif_performance
## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##     P 3.84 [2.41,  6.73]     1.96      0.26     [0.15, 0.42]
## 
## Moderate Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##     C 7.63 [4.49, 13.61]     2.76      0.13     [0.07, 0.22]
##     M 9.45 [5.49, 16.91]     3.07      0.11     [0.06, 0.18]
plot(vif_performance)

Resultado:
Como se observa en la gráfica, el VIF de las variables C y M superan el umbral de 5, lo que indica que estas dos variable comparten demasiada información, la correlación es alta. El VIF de la variable P se encuentre dentro del límite aceptable, es moderado.

2.4.2. VIF con “car”

Usando la función vif de la librería “car”.

library(car)
vif_car <- vif(modelo_ventas)
vif_car
##        C        P        M 
## 7.631451 3.838911 9.449210

Resultado:
Los resultados son convergentes a los obtenidos con la librería “performance”, el VIF de las variables C y M superan el umbral de 5, lo que indica que estas dos variable comparten demasiada información, la correlación es alta. El VIF de la variable P se encuentre dentro del límite aceptable, es moderado.

2.4.3. VIF con “mctest”

Usando la función mc.plot de la librería “mctest”.

library(mctest)
mc.plot(modelo_ventas, type = "vif", vif = 5)

Resultado:
Los resultados son convergentes a los obtenidos con la librería “performance”, el VIF de las variables C y M superan el umbral de 5, lo que indica que estas dos variable comparten demasiada información, la correlación es alta. El VIF de la variable P se encuentre dentro del límite aceptable, es moderado.

3. Supuesto de Homocedasticidad de los Residuos

3.1. Prueba de White con “lmtest”

Utilizando la función bptest de la librería “lmtest”, con los términos cruzados incluidos.

library(lmtest)
white_test <- bptest(modelo_ventas, ~ (C + P + M)^2 + I(C^2) + I(P^2) + I(M^2), data = ventas_empresa)
white_test
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ventas
## BP = 7.1227, df = 9, p-value = 0.6244

Resultado:
Como \(p=0.6244>\alpha=0.05\), no se rechaza la hipótesis nula. La varianza de los errores es constante. Hay evidencia de que la varianza de los residuos es homocedástica.

3.2. Gráfica del estadístico de la prueba de White

El estadístico de la prueba de White sigue asintóticamente una distribución Chi-cuadrado (χ2\(\chi^2\)) con \(gl=k(k+1)/2\) grados de libertad, gráficamente se puede observar de la siguiente manera:

library(fastGraph)
white_stat <- white_test$statistic
gl_white <- white_test$parameter

shadeDist(
  xstat = white_stat, 
  ddist = "dchisq",
  parm1 = gl_white,
  lower.tail = FALSE,
  main = "Resultados de la Prueba de White (Chi-cuadrado)",
  sub = paste("Valor critico:", round(VC, 2), "","White", round(white_stat, 2)),
)

Este resultado refuerza el pasado, como \(LM_W=7.12<\text{V.C.}=7.81\), no se rechaza la hipótesis nula. La varianza de los errores es constante. Hay evidencia de que la varianza de los residuos es homocedástica.

4. Supuesto de No-Autocorrelación en los Residuos

4.1. Prueba del Múltiplicador de Lagrange [Breusch-Godfrey(BG)]

4.1.1. Prueba BG de Primer Orden con “lmtest”

library(lmtest)
bgtest(modelo_ventas, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_ventas
## LM test = 2.5963, df = 1, p-value = 0.1071

Resultado:
Como \(p=0.1071>\alpha=0.05\), no se rechaza la hipótesis nula, no hay evidencia de autocorrelación de primer orden en los residuos del modelo.

4.1.2. Prueba BG de Segundo Orden con “lmtest”

library(lmtest)
bgtest(modelo_ventas, order = 2)
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  modelo_ventas
## LM test = 3.8409, df = 2, p-value = 0.1465

Resultado:
Como \(p=0.1456>\alpha=0.05\), no se rechaza la hipótesis nula, no hay evidencia de autocorrelación de segundo orden los residuos del modelo.