multicolinealidad- BM22010
Gabriela Belloso
2024-05-12
library(wooldridge)
data(hprice1)
head(force(hprice1),n=5) #mostrar las primeras 5 observaciones## 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.8296301. estime el siguiente modelo:
\[ price = \hat{\alpha} + \hat{\alpha}_1 \times \text{lotsize} + \hat{\alpha}_2 \times \text{sqrft} + \hat{\alpha}_3 \times \text{bdrms} + \epsilon\]
# Creación del modelo:
modelo <- lm(formula= price ~ lotsize + sqrft + bdrms , data= hprice1)
library(stargazer)
stargazer(modelo,title = "Modelo", type = "html")| 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.01 |
2. Verifique si hay evidencia de la independencia de los regresores (no colinealidad), a través de:
a) Indice de condición y prueba de FG, presente sus resultados de manera tabular en ambos casos y para la prueba de FG presente también sus resultados de forma gráfica usando la librería fastGraph
ÍNDICE DE CONDICIÓN
cálculo manual del índice de condición
\[ \text{Calcular la matriz: } \mathbf{X}^{\text{t}} \mathbf{X} \]
# Cargar la librería stargazer
library(stargazer)
# Crear la matriz de diseño a partir del modelo
X_mat <- model.matrix(modelo)
# Imprimir las primeras 6 filas de la matriz de diseño en formato HTML
stargazer(head(X_mat, n = 6), type = "html") | (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 |
# Calcular la matriz de producto de la matriz de diseño transpuesta por sí misma
XX_matrix <- t(X_mat) %*% X_mat
# Imprimir la matriz de producto cruzado en formato HTML
stargazer(XX_matrix, type = "html") | (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 |
\[ \text{Normalizar la matriz: } \mathbf{X}^{\text{t}} \mathbf{X} \]
library(stargazer)
options(scipen = 999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "html")| 0.107 | 0 | 0 | 0 |
| 0 | 0.00001 | 0 | 0 |
| 0 | 0 | 0.0001 | 0 |
| 0 | 0 | 0 | 0.029 |
\[ \text{matriz normalizada: } \mathbf{X}^{\text{t}} \mathbf{X} \]
library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "html",digits = 4)| 1 | 0.6655 | 0.9617 | 0.9736 |
| 0.6655 | 1 | 0.6776 | 0.6712 |
| 0.9617 | 0.6776 | 1 | 0.9696 |
| 0.9736 | 0.6712 | 0.9696 | 1 |
\[ \text{Autovalores normalizados de: } \mathbf{X}^{\text{t}} \mathbf{X} \]
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "html")| 3.482 | 0.455 | 0.039 | 0.025 |
\[ \text{Cálculo de: } K = \sqrt{\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}}\]
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)## [1] 11.86778Como \(\kappa(x) \leq 20\) se considera que la multicolinealidad es leve.
Índice de condición Usando librería “mctest”
library(mctest)
X_mat<-model.matrix(modelo)
mctest(mod = modelo)##
## 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 librería “olsrr”
library(olsrr)
ols_eigen_cindex(model = modelo)## 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 manual
\[ \text{Cálculo de } |\mathbf{R}| \]
library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "html")| 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 |
\[ \text{Calcular la matriz } R \]
library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "html",digits = 4)| lotsize | sqrft | bdrms | |
| lotsize | 1 | 0.1838 | 0.1363 |
| sqrft | 0.1838 | 1 | 0.5315 |
| bdrms | 0.1363 | 0.5315 | 1 |
#También se puede calcular R a través de cor(X_mat[,-1])
cor(X_mat[,-1])## lotsize sqrft bdrms
## lotsize 1.0000000 0.1838422 0.1363256
## sqrft 0.1838422 1.0000000 0.5314736
## bdrms 0.1363256 0.5314736 1.0000000\[ \text{Calcular } |\mathbf{R}| \]
determinante_R<-det(R)
print(determinante_R)## [1] 0.6917931Interpretación: Se sospecha que no exixte colinealidad
\[ \text{Estadístico } \chi^2_{FG} \]
m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)## [1] 31.38122\[ \text{Valor Crítico} \]
gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)## [1] 7.814728# También se puede hacer de la ssiguiente manera
vc <- qchisq(p = 0.95, df = gl, lower.tail = TRUE)
print(vc)## [1] 7.814728\[ \text{Regla de desición} \]
Como \(\chi^2_{FG} \geq V.C.\), se rechaza \(H_0\), por lo tanto hay evidencia de colinealidad en los regresores.
Cálculo de FG usando “mctest”
library(mctest)
mctest::omcdiag(mod = modelo)##
## Call:
## mctest::omcdiag(mod = modelo)
##
##
## 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)
alpha_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 la “psych”
library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test)## $chisq
## [1] 31.38122
##
## $p.value
## [1] 0.0000007065806
##
## $df
## [1] 3chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
gl<-m*(m-1)/2
vc<-qchisq(p = 0.95,df = gl)
alpha_sig <- 0.05
library(fastGraph)
shadeDist(chi_FG, ddist = "dchisq", parm1 = gl, lower.tail = FALSE, xmin = 0, xlab = "Valor de chi-cuadrado",
main = "Prueba de Bartlett")
abline(v = vc, col = "red", lty = 2)
axis(1, at = vc, labels = paste("vc:", round(vc, 2)), col.axis = "black", las = 1)
if (chi_FG > vc) {
text(x = vc + 0, y = 0.22, labels = "Rechazar H0", col = "blue", cex = 0.8)
} else {
text(x = vc + 0, y = 0.22, labels = "No rechazar H0", col = "blue", cex = 0.8)
}
text(vc, 0, expression(alpha == 0.05), pos = 4, col = "black", cex = 0.8)
# valor critico de la prueba
FG_test$p.value## [1] 0.0000007065806b) Factores inflacionarios de la varianza, presente sus resultados de forma tabular y de forma gráfica.
cálculo de los VIF
\[ \text{Referencia entre } R^2_j \]
library(dplyr)
R.cuadrado.regresores<-c(0,0.5,.8,.9)
as.data.frame(R.cuadrado.regresores) %>% mutate(VIF=1/(1-R.cuadrado.regresores))## R.cuadrado.regresores VIF
## 1 0.0 1
## 2 0.5 2
## 3 0.8 5
## 4 0.9 10Cálculo manual
\[ \text{Matriz de correlación } \]
# Ya que se calculó anteiormente, se muestra la matriz de correlación del modelo.
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\[ \text{Inversa de la matriz de correlación } R^{-1} \]
inversa_R<-solve(R)
print(inversa_R)## lotsize sqrft bdrms
## lotsize 1.03721145 -0.1610145 -0.05582352
## sqrft -0.16101454 1.4186543 -0.73202696
## bdrms -0.05582352 -0.7320270 1.39666321\[ \text{VIF's para el modelo estimado:} \]
VIFs<-diag(inversa_R)
print(VIFs)## lotsize sqrft bdrms
## 1.037211 1.418654 1.396663Cálculo de los VIF’s usando “performance”
library(performance)
VIFs<-multicollinearity(x = modelo,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)
print(VIFs_car)## lotsize sqrft bdrms
## 1.037211 1.418654 1.396663Cálculo de los VIF’s usando “mctest”
library(mctest)
mc.plot(mod = modelo,vif = 2)