Pruebas de Multicolinealidad.
CARLOS ELIAS RAMOS GONZALEZ
8/6/2020
Estimación del modelo
options(scipen = 9999)
library(haven)
library(stargazer)
hprice1 <- read_dta("C:/Users/manue/Desktop/Practicas_ECMA/datos_Rdata/Datos Examen/hprice1.dta")
myformula<-as.formula("price~assess+bdrms+lotsize+colonial+llotsize")
modelo_estimado<-lm(myformula,data = hprice1)
stargazer(modelo_estimado,type = "html",title = "modelo estimado")| Dependent variable: | |
| price | |
| assess | 0.940*** |
| (0.072) | |
| bdrms | 8.620 |
| (6.791) | |
| lotsize | 0.001 |
| (0.001) | |
| colonial | 10.031 |
| (10.580) | |
| llotsize | -13.357 |
| (17.813) | |
| Constant | 68.090 |
| (146.133) | |
| Observations | 88 |
| R2 | 0.832 |
| Adjusted R2 | 0.822 |
| Residual Std. Error | 43.364 (df = 82) |
| F Statistic | 81.224*** (df = 5; 82) |
| Note: | p<0.1; p<0.05; p<0.01 |
Indice de condicion( caluclo manual)
Calcular la matriz |X X|
# matriz Xmat
Xmat<-model.matrix(modelo_estimado)
print(Xmat)## (Intercept) assess bdrms lotsize colonial llotsize
## 1 1 349.1 4 6126 1 8.720297
## 2 1 351.5 3 9903 1 9.200593
## 3 1 217.7 3 5200 0 8.556414
## 4 1 231.8 3 4600 1 8.433811
## 5 1 319.1 4 6095 1 8.715224
## 6 1 414.5 5 8566 1 9.055556
## 7 1 367.8 3 9000 1 9.104980
## 8 1 300.2 3 6210 1 8.733916
## 9 1 236.1 3 6000 0 8.699514
## 10 1 256.3 3 2892 0 7.969704
## 11 1 314.0 4 6000 1 8.699514
## 12 1 416.5 5 7047 1 8.860357
## 13 1 434.0 3 12237 1 9.412219
## 14 1 279.3 3 6460 0 8.773385
## 15 1 287.5 3 6519 1 8.782476
## 16 1 232.9 4 3597 1 8.187856
## 17 1 303.8 4 5922 0 8.686430
## 18 1 305.6 3 7123 1 8.871084
## 19 1 266.7 3 5642 1 8.637994
## 20 1 326.0 4 8602 1 9.059750
## 21 1 294.3 3 5494 1 8.611412
## 22 1 318.8 3 7800 1 8.961879
## 23 1 294.2 3 6003 0 8.700015
## 24 1 208.0 4 5218 0 8.559870
## 25 1 239.7 3 9425 1 9.151121
## 26 1 294.1 3 6114 0 8.718336
## 27 1 267.4 3 6710 0 8.811355
## 28 1 359.9 3 8577 1 9.056840
## 29 1 478.1 7 8400 1 9.035987
## 30 1 355.3 4 9773 1 9.187379
## 31 1 217.8 4 4806 1 8.477620
## 32 1 385.0 4 15086 0 9.621523
## 33 1 224.3 3 5763 1 8.659213
## 34 1 251.9 4 6383 1 8.761394
## 35 1 354.9 4 9000 1 9.104980
## 36 1 212.5 4 3500 0 8.160519
## 37 1 452.4 4 10892 1 9.295784
## 38 1 518.1 5 15634 1 9.657204
## 39 1 289.4 4 6400 1 8.764053
## 40 1 268.1 2 8880 0 9.091557
## 41 1 278.5 3 6314 1 8.750525
## 42 1 655.4 5 28231 1 10.248176
## 43 1 273.3 4 7050 1 8.860783
## 44 1 212.1 3 5305 0 8.576406
## 45 1 354.0 5 6637 1 8.800415
## 46 1 252.1 3 7834 1 8.966228
## 47 1 324.0 3 1000 0 6.907755
## 48 1 475.5 4 8112 0 9.001100
## 49 1 256.8 3 5850 1 8.674197
## 50 1 279.2 4 6660 1 8.803875
## 51 1 313.9 3 6637 1 8.800415
## 52 1 279.8 2 15267 0 9.633449
## 53 1 198.7 3 5146 1 8.545975
## 54 1 221.5 3 6017 1 8.702344
## 55 1 268.4 3 8410 1 9.037177
## 56 1 282.3 4 5625 1 8.634976
## 57 1 230.7 4 5600 1 8.630522
## 58 1 287.0 4 6525 1 8.783396
## 59 1 298.7 3 6060 1 8.709465
## 60 1 314.6 4 5539 0 8.619569
## 61 1 291.0 3 7566 0 8.931419
## 62 1 286.4 4 5484 1 8.609590
## 63 1 253.6 6 5348 1 8.584478
## 64 1 482.0 5 15834 1 9.669915
## 65 1 384.3 4 8022 1 8.989944
## 66 1 543.6 4 11966 1 9.389825
## 67 1 336.5 4 8460 1 9.043104
## 68 1 515.1 4 15105 1 9.622781
## 69 1 437.0 4 10859 0 9.292749
## 70 1 263.4 3 6300 1 8.748305
## 71 1 300.4 3 11554 0 9.354787
## 72 1 250.7 3 6000 1 8.699514
## 73 1 708.6 5 31000 0 10.341743
## 74 1 276.3 3 4054 1 8.307459
## 75 1 388.6 2 20700 0 9.937889
## 76 1 252.5 3 5525 0 8.617039
## 77 1 295.2 4 92681 1 11.436919
## 78 1 359.5 3 8178 1 9.009203
## 79 1 276.2 4 5944 1 8.690138
## 80 1 249.8 3 18838 0 9.843632
## 81 1 202.4 4 4315 1 8.369853
## 82 1 254.0 3 5167 1 8.550048
## 83 1 306.8 4 7893 1 8.973732
## 84 1 318.3 3 6056 1 8.708805
## 85 1 259.4 3 5828 0 8.670429
## 86 1 258.1 3 6341 0 8.754792
## 87 1 232.0 2 6362 0 8.758098
## 88 1 252.0 4 4950 1 8.507143
## attr(,"assign")
## [1] 0 1 2 3 4 5# matriz XXmat
XXmat<-t(Xmat)%*%Xmat
print(XXmat)## (Intercept) assess bdrms lotsize colonial
## (Intercept) 88.0000 27784.8 314.000 793748 61.0000
## assess 27784.7998 9563052.8 102507.499 278300049 19578.8999
## bdrms 314.0000 102507.5 1182.000 2933767 228.0000
## lotsize 793748.0000 278300049.1 2933767.000 16165159010 555967.0000
## colonial 61.0000 19578.9 228.000 555967 61.0000
## llotsize 783.6492 250005.6 2802.953 7457452 544.0597
## llotsize
## (Intercept) 783.6492
## assess 250005.6359
## bdrms 2802.9529
## lotsize 7457452.0339
## colonial 544.0597
## llotsize 7004.2300## sn calculo de la matriz de normalizacion
sn<-diag(1/sqrt(diag(XXmat)))
print(sn)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1066004 0.0000000000 0.00000000 0.000000000000 0.0000000 0.00000000
## [2,] 0.0000000 0.0003233715 0.00000000 0.000000000000 0.0000000 0.00000000
## [3,] 0.0000000 0.0000000000 0.02908649 0.000000000000 0.0000000 0.00000000
## [4,] 0.0000000 0.0000000000 0.00000000 0.000007865204 0.0000000 0.00000000
## [5,] 0.0000000 0.0000000000 0.00000000 0.000000000000 0.1280369 0.00000000
## [6,] 0.0000000 0.0000000000 0.00000000 0.000000000000 0.0000000 0.01194868## XXmat Normalizada
XX_norm<-(sn%*%XXmat)%*%sn
print(XX_norm)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 0.9577842 0.9735978 0.6655050 0.8325754 0.9981600
## [2,] 0.9577842 1.0000000 0.9641589 0.7078236 0.8106345 0.9659871
## [3,] 0.9735978 0.9641589 1.0000000 0.6711613 0.8491046 0.9741523
## [4,] 0.6655050 0.7078236 0.6711613 1.0000000 0.5598789 0.7008423
## [5,] 0.8325754 0.8106345 0.8491046 0.5598789 1.0000000 0.8323413
## [6,] 0.9981600 0.9659871 0.9741523 0.7008423 0.8323413 1.0000000## autovalores
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)$values
print(lambdas)## [1] 5.1933823834 0.4923113019 0.2365049515 0.0490596083 0.0282837924
## [6] 0.0004579625# Indice de condicion
K<-sqrt(max(lambdas)/min(lambdas))
print(K)## [1] 106.4903Uso de la libreria “mctest”, para obtener el indice de condición
library(mctest)
library(mctest)
source("C:/Users/manue/Desktop/Practicas_ECMA/datos_Rdata/correccion_eigprop.R")
my_eigprop(mod=modelo_estimado)##
## Call:
## my_eigprop(mod = modelo_estimado)
##
## Eigenvalues CI (Intercept) assess bdrms lotsize colonial llotsize
## 1 5.1934 1.0000 0.0000 0.0014 0.0012 0.0037 0.0079 0.0000
## 2 0.4923 3.2479 0.0000 0.0001 0.0015 0.2963 0.0615 0.0000
## 3 0.2365 4.6860 0.0003 0.0113 0.0041 0.0322 0.8545 0.0002
## 4 0.0491 10.2888 0.0048 0.4576 0.0148 0.0101 0.0021 0.0027
## 5 0.0283 13.5505 0.0012 0.2138 0.9074 0.0093 0.0696 0.0017
## 6 0.0005 106.4903 0.9936 0.3158 0.0711 0.6484 0.0044 0.9954
##
## ===============================
## Row 5==> bdrms, proportion 0.907364 >= 0.50
## Row 6==> lotsize, proportion 0.648440 >= 0.50
## Row 3==> colonial, proportion 0.854491 >= 0.50
## Row 6==> llotsize, proportion 0.995427 >= 0.50Prueba de Farrar-Glaubar ( manual)
library(fastGraph)
m<-ncol(Xmat[,-1]) # cantidad de variables explicativas k-1
n<-nrow(Xmat)
determinante_R<- det(cor(Xmat[,-1])) # determinanre de la matriz de correlacion
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)## [1] 164.9525## Valor Critico
gl<-m*(m-1)/2
VC<-qchisq(0.05,gl,lower.tail = FALSE)
print(VC)## [1] 18.30704shadeDist(xshade = chi_FG,ddist = "dchisq",parm1 = gl,lower.tail = FALSE,sub=paste("VC:",VC,"FG:",chi_FG))
Como χ2FG≥V.C. se rechaza H0, por lo tanto hay evidencia de colinealidad en los regresores.
Prueba de Farrar-Glaubar con mctest
library(mctest)
mctest(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.1420 0
## Farrar Chi-Square: 164.9525 1
## Red Indicator: 0.3832 0
## Sum of Lambda Inverse: 12.3289 0
## Theil's Method: -0.8940 0
## Condition Number: 25.0893 0
##
## 1 --> COLLINEARITY is detected by the test
## 0 --> COLLINEARITY is not detected by the testPrueba de Farrar-Glaubar Uso de la libreria psych
library(psych)
library(fastGraph)
FG_test<-cortest.bartlett(Xmat[,-1])
VC_1<-qchisq(0.05,FG_test$df,lower.tail = FALSE)
print(FG_test)$chisq [1] 164.9525
$p.value [1] 0.000000000000000000000000000003072151
$df [1] 10
shadeDist(xshade = FG_test$chisq,ddist = "dchisq",parm1 = FG_test$df,lower.tail = FALSE,sub=paste("VC:",VC_1,"FG:",FG_test$chisq))
Cálculando los VIF para el modelo estimado(manual)
Matriz de Correlación de los regresores del modelo (Como se obtuvo con anterioridad):
VIF<-diag(solve(cor(Xmat[,-1])))
print(VIF)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674Uso de la librería “car” y “mctest”
Obtención de los VIF’s, a través de la librería “car”
library(car)
VIFs_car<-vif(modelo_estimado)
print(VIFs_car)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674Obtención de los VIF’s, a través de la librería “mctest”
library(mctest)
mc.plot(modelo_estimado,vif = 2,)