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.
library(readxl)
ventas_empresa <- read_excel("D:/harold/Descargas/ventas_empresa.xlsx")
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
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.
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.
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.
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.
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\)].
Más bien, verificar que la relación entre las variables
explicativas no resulte problemática para su
interpretación.
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.
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.
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.
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.
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)))
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.
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.
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.
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.
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.
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.
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.