Ejercicio de multicolinealidad
DANIEL HERNANDEZ
2023-05-15
Carga de datos
library(wooldridge)
data("hprice1")
head(force(hprice1),n=5)## price assess bdrms lotsize sqrft colonial lprice lassess llotsize lsqrft
## 1 300 349.1 4 6126 2438 1 5.703783 5.855359 8.720297 7.798934
## 2 370 351.5 3 9903 2076 1 5.913503 5.862210 9.200593 7.638198
## 3 191 217.7 3 5200 1374 0 5.252274 5.383118 8.556414 7.225482
## 4 195 231.8 3 4600 1448 1 5.273000 5.445875 8.433811 7.277938
## 5 373 319.1 4 6095 2514 1 5.921578 5.765504 8.715224 7.829630Modelo Estimado
library(stargazer)
modelo_estimado<-lm(formula = price~lotsize+sqrft+bdrms,data = hprice1)
stargazer(modelo_estimado,type = "text",title = "Modelo Estimado")##
## Modelo Estimado
## ===============================================
## Dependent variable:
## ---------------------------
## price
## -----------------------------------------------
## lotsize 0.002***
## (0.001)
##
## sqrft 0.123***
## (0.013)
##
## bdrms 13.853
## (9.010)
##
## Constant -21.770
## (29.475)
##
## -----------------------------------------------
## Observations 88
## R2 0.672
## Adjusted R2 0.661
## Residual Std. Error 59.833 (df = 84)
## F Statistic 57.460*** (df = 3; 84)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Indice de condición
#Matriz X
mat_X<-model.matrix(modelo_estimado)
stargazer(head(mat_X,n=6),type = "text",title = "Matriz X")##
## Matriz X
## =================================
## (Intercept) lotsize sqrft bdrms
## ---------------------------------
## 1 1 6,126 2,438 4
## 2 1 9,903 2,076 3
## 3 1 5,200 1,374 3
## 4 1 4,600 1,448 3
## 5 1 6,095 2,514 4
## 6 1 8,566 2,754 5
## ---------------------------------#Matriz XX
mat_XX<-t(mat_X)%*%mat_X
stargazer(mat_XX,type = "text",title = "Matriz XX")##
## Matriz XX
## ==============================================================
## (Intercept) lotsize sqrft bdrms
## --------------------------------------------------------------
## (Intercept) 88 793,748 177,205 314
## lotsize 793,748 16,165,159,010 1,692,290,257 2,933,767
## sqrft 177,205 1,692,290,257 385,820,561 654,755
## bdrms 314 2,933,767 654,755 1,182
## --------------------------------------------------------------#Normalización de la matriz XX
options(scipen = 999)
Sn<-solve(diag(sqrt(diag(mat_XX))))
stargazer(Sn,type = "text",title = "Matriz de normalización")##
## Matriz de normalización
## ==========================
## 0.107 0 0 0
## 0 0.00001 0 0
## 0 0 0.0001 0
## 0 0 0 0.029
## --------------------------#Matriz XX normalizada
normalizada_XX<-(Sn%*%mat_XX)%*%Sn
stargazer(normalizada_XX,type = "text",title = "Matriz Normalizada")##
## Matriz Normalizada
## =======================
## 1 0.666 0.962 0.974
## 0.666 1 0.678 0.671
## 0.962 0.678 1 0.970
## 0.974 0.671 0.970 1
## -----------------------#Autovalores de la matriz XX normalizada
lambdas<-eigen(normalizada_XX,symmetric = TRUE)
stargazer(lambdas$values,type = "text",title = "Autovalores de matriz XX normalizada")##
## Autovalores de matriz XX normalizada
## =======================
## 3.482 0.455 0.039 0.025
## -----------------------#Cálculo de K (número de condición)
K<-sqrt(max(lambdas$values)/min(lambdas$values))
stargazer(K,title = "K (número de condición)",type = "text")##
## K (número de condición)
## ======
## 11.868
## ------Como k(x)=<20 se considera multicolinealidad leve, no sería un problema.
Cálculo del Indice de Condición usando la librería “mctest”
library(mctest)
matriz_X<-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.6918 0
## Farrar Chi-Square: 31.3812 1
## Red Indicator: 0.3341 0
## Sum of Lambda Inverse: 3.8525 0
## Theil's Method: -0.7297 0
## Condition Number: 11.8678 0
##
## 1 --> COLLINEARITY is detected by the test
## 0 --> COLLINEARITY is not detected by the testCálculo del Indice de Condición usando la librería “olsrr”
library(olsrr)
ols_eigen_cindex(model = modelo_estimado)## Eigenvalue Condition Index intercept lotsize sqrft bdrms
## 1 3.48158596 1.000000 0.003663034 0.0277802824 0.004156293 0.002939554
## 2 0.45518380 2.765637 0.006800735 0.9670803174 0.006067321 0.005096396
## 3 0.03851083 9.508174 0.472581427 0.0051085488 0.816079307 0.016938178
## 4 0.02471941 11.867781 0.516954804 0.0000308514 0.173697079 0.975025872Prueba de Farrar-Glaubar
#Cálculo de R
Zn<-scale(mat_X[,-1])
stargazer(head(Zn,n=6),type = "text",title = "Cálculo de R")##
## Cálculo de R
## =======================
## lotsize sqrft bdrms
## -----------------------
## 1 -0.284 0.735 0.513
## 2 0.087 0.108 -0.675
## 3 -0.375 -1.108 -0.675
## 4 -0.434 -0.980 -0.675
## 5 -0.287 0.867 0.513
## 6 -0.045 1.283 1.702
## -----------------------#Calcular la matriz R
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "text",title = "Calcular la matriz R",digits = 4)##
## Calcular la matriz R
## =============================
## lotsize sqrft bdrms
## -----------------------------
## lotsize 1 0.1838 0.1363
## sqrft 0.1838 1 0.5315
## bdrms 0.1363 0.5315 1
## -----------------------------#Calcular R
determinante_R<-det(R)
stargazer(determinante_R,type = "text",title = "R")##
## R
## =====
## 0.692
## -----Aplicando la prueba Farrer Glaubar (Barlett)
m<-ncol(mat_X[,-1])
n<-nrow(mat_X[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)## [1] 31.38122Valor Critico
gl<-m*(m-1)/2
VC<-qchisq(p=0.95,df=gl)
print(VC)## [1] 7.814728Regla de decisión: Como chi_FG >= VC, se rechaza la Ho, por lo tanto hay evidencia de colinealidad entre los regresores
Cálculo de FG utilizando “mctest”
library(mctest)
mctest::omcdiag(mod = modelo_estimado)##
## Call:
## mctest::omcdiag(mod = modelo_estimado)
##
##
## Overall Multicollinearity Diagnostics
##
## MC Results detection
## Determinant |X'X|: 0.6918 0
## Farrar Chi-Square: 31.3812 1
## Red Indicator: 0.3341 0
## Sum of Lambda Inverse: 3.8525 0
## Theil's Method: -0.7297 0
## Condition Number: 11.8678 0
##
## 1 --> COLLINEARITY is detected by the test
## 0 --> COLLINEARITY is not detected by the testlibrary(fastGraph)## Warning: package 'fastGraph' was built under R version 4.2.3alpha_sig<-0.05
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
GL<-gl
vc<-qchisq(1-alpha_sig,gl,lower.tail=TRUE)
shadeDist(chi_FG,ddist = "dchisq",
parm1 = gl,
lower.tail = FALSE, xmin = 0,
sub=paste("VC:",round(VC,2),"","chi_FG:",round(chi_FG,2)))
Cálculo de FG usando “psych”
library(psych)
FG_test<-cortest.bartlett(mat_X[,-1])
print(FG_test)## $chisq
## [1] 31.38122
##
## $p.value
## [1] 0.0000007065806
##
## $df
## [1] 3Factores Inflacionarios de la Varianza (FIV)
print(R)## lotsize sqrft bdrms
## lotsize 1.0000000 0.1838422 0.1363256
## sqrft 0.1838422 1.0000000 0.5314736
## bdrms 0.1363256 0.5314736 1.0000000#Inversa de la matriz de correlación R
R_inversa<-solve(R)
stargazer(R_inversa,title = "Inversa de la matriz de correlación R",type = "text")##
## Inversa de la matriz de correlación R
## =============================
## lotsize sqrft bdrms
## -----------------------------
## lotsize 1.037 -0.161 -0.056
## sqrft -0.161 1.419 -0.732
## bdrms -0.056 -0.732 1.397
## -----------------------------#VIF`s para el modelo estimado
VIFs<-diag(R_inversa)
stargazer(VIFs,type = "text",title = "VIF`s para el modelo estimado")##
## VIF`s para el modelo estimado
## ===================
## lotsize sqrft bdrms
## -------------------
## 1.037 1.419 1.397
## -------------------Cálculo de los VIF`s usando “performance”
library(performance)
VIFS<-multicollinearity(x = modelo_estimado,verbose = FALSE)
VIFS## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## lotsize 1.04 [1.00, 11.02] 1.02 0.96 [0.09, 1.00]
## sqrft 1.42 [1.18, 1.98] 1.19 0.70 [0.51, 0.85]
## bdrms 1.40 [1.17, 1.95] 1.18 0.72 [0.51, 0.86]plot(VIFS)
Cálculo de los VIF`s usando “car”
library(car)
VIFS_Car<-vif(modelo_estimado)
stargazer(VIFS_Car,title = "VIF`s usando",type = "text")##
## VIF`s usando
## ===================
## lotsize sqrft bdrms
## -------------------
## 1.037 1.419 1.397
## -------------------Cálculo de los VIF`s usando “mctest”
mc.plot(mod = modelo_estimado,vif = 2)