Practica pruebas de multicolinealidad
MELVIN ORLANDO LUNA MARTINEZ
25 de mayo de 2020
Importancion de datos.
options(scipen = 99999)
library(haven)
hprice1 <- read_dta("C:/Users/melvi/Desktop/Econometria/Datos/hprice1.dta")
myformula<-as.formula("price~assess+bdrms+lotsize+colonial+llotsize")
modelo_de_regresion<-lm(formula = price~assess+bdrms+lotsize+colonial+llotsize,data = hprice1)
summary(modelo_de_regresion)##
## Call:
## lm(formula = price ~ assess + bdrms + lotsize + colonial + llotsize,
## data = hprice1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -104.837 -21.870 0.209 19.143 202.718
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.090453 146.133300 0.466 0.642
## assess 0.940298 0.071579 13.137 <0.0000000000000002 ***
## bdrms 8.619563 6.790853 1.269 0.208
## lotsize 0.001087 0.000818 1.329 0.188
## colonial 10.031313 10.580456 0.948 0.346
## llotsize -13.357148 17.813399 -0.750 0.455
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.36 on 82 degrees of freedom
## Multiple R-squared: 0.832, Adjusted R-squared: 0.8218
## F-statistic: 81.22 on 5 and 82 DF, p-value: < 0.00000000000000022Indice de Condicion.
# calculo manual
# Xmat
Xmat<-model.matrix(modelo_de_regresion)
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# 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 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_norm
XXmat_norm<-(Sn%*%XXmat)%*%Sn
print(XXmat_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*
lambas<-eigen(XXmat_norm,symmetric = TRUE)$values
print(lambas)## [1] 5.1933823834 0.4923113019 0.2365049515 0.0490596083 0.0282837924
## [6] 0.0004579625# Indice de Condicion
K<-sqrt(max(lambas)/min(lambas))
print(K)## [1] 106.4903Indice de condicion usando mctess
Para trabajar con la libreria mctest es necesario cargar la correcion de la nueva actualizacion.
library(mctest)
source(file = 'C:/Users/melvi/Desktop/Econometria/correccion_eigprop.R')
my_eigprop(mod = modelo_de_regresion)##
## Call:
## my_eigprop(mod = modelo_de_regresion)
##
## 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 Gloubar.
library(fastGraph)
m<-ncol(Xmat[,-1]) #Cantidad de variables explicativas k-1
n<-nrow(Xmat)
determinante_R<-det(cor(Xmat[,-1]))
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.30704#Nuestro estadistico de prueba es mayor que el VC, por lo tanto se rechaza la H0.
shadeDist(xshade = Chi_FG,ddist = "dchisq",gl,lower.tail = FALSE,sub=paste("VC:",VC,"FG:",Chi_FG))
Farrar Gloubar con mc test
library(mctest)
mctest(modelo_de_regresion)##
## 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 testlibrary(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))
VIF (Factores Inflacionarios de la varianza) MANUAL. [1 forma]
VIF<-diag(solve(cor(Xmat[,-1])))
print(VIF)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674VIF con library car [2 forma]
library(car)
VIF_car<-vif(modelo_de_regresion)
print(VIF_car)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674VIF con mctest [3 forma]
library(mctest)
mc.plot(modelo_de_regresion,vif = 2)