Ejercicios De Multicolinealidad

Ezequiel Benjamín López Coto - LC22057

Datos del modelo.

#Carga de los datos.
library(wooldridge)
library(haven)
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.829630

1. Modelo Estimado.

#Modelo estimado.
library(stargazer)
Estimación_del_modelo<-lm(price~lotsize+sqrft+bdrms,data = hprice1)
stargazer(Estimación_del_modelo,type = "html",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.01

2. Verifique si hay evidencia de la independencia de los regresores (no colinealidad), a través de:

Indice de Condición.

El número de condición mide la sensibilidad de las estimaciones mínimo cuadráticas ante pequeños cambios en los datos.

El número de condición es igual a:

\[\kappa(x) = \sqrt{\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}}\]

Dónde:

λmin = la raíz característica más pequeña.

λmax = la raíz característica más grande.

Interpretación:

Si \(k(x) \leq 20\) Si \(20 < k(x) < 30\) Si \(k(x) \geq 30\)
La multicolinealidad es leve y no se considera un problema. La multicolinealidad se considera moderada. La multicolinealidad es severa.

Cálculo del Indice de Condición de forma “manual”

Se basa la matriz Xt X.

library(stargazer)
X_mat<-model.matrix(Estimación_del_modelo)
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 
XX_matrix<-t(X_mat)%*%X_mat
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

La matriz Xt X debe ser normalizada, para evitar sesgo de la escala de las variables.

Cálculo de la matriz de normalización: \[ S_{n} = \begin{cases} 1 & \text{si } i = j \\ \frac{X^tX_{ij}}{\sqrt{0}} & \text{si } i \neq j \end{cases} \]

library(stargazer)
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

Xt X normalizada:

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

Autovalores de Xt X normalizada:

#Autovalores
library(stargazer)
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "html")
3.482 0.455 0.039 0.025

Cálculo de \[\kappa(x) = \sqrt{\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}}\]

K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)
## [1] 11.86778

En base a la interpretación antes vista, podemos concluir que se considera que la multicolinealidad es severa, ya que k(x) > 30.

Cálculo del Indice de Condición usando librería “mctest”

#Utilizando a librería "mctest"
library(mctest)
X_mat<-model.matrix(Estimación_del_modelo)
mctest(mod = Estimación_del_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 test

Cálculo del Indice de Condición usando librería “olsrr”

#Utilizando a librería "olsrr”
library(olsrr)
ols_eigen_cindex(model = Estimación_del_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.975025872

Prueba de Farrar-Glaubar.

(Implementación de la prueba de Bartlett) Esta prueba identifica si a nivel poblacional, los regresores del modelo presentan independencia estadistica (son ortogonales), a través de la matriz de correlación muestral R, y se verifica si a nivel poblacional dicha matriz de correlación corresponde a una matriz identidad, las hipótesis de la prueba son las siguientes:

\[ \begin{cases} H_0: & R = I \\ H_1: & R \neq I \end{cases} \]

Si no se rechaza H0 Si se rechaza H0
No hay evidencia de multicolinealidad Hay evidencia de multicolinealidad

El contraste se realiza a través de: \(\chi^2_{FG} = -(n-1-2m+56) \ln(|R|)\)

con, \(m=k−1, gl=m(m−1)/2\)

Regla de decisión:

Rechazar \(H_0\) si \(\chi^2_{FG} \geq V.C.\), o si \(p \leq \alpha\).

Cálculo de la Prueba de Farrar-Glaubar de forma “manual”

Cálculo de |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

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

Calcular |R|.

determinante_R<-det(R)
print(determinante_R)
## [1] 0.6917931

Aplicando la prueba de Farrer Glaubar (Bartlett).

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

Valor Crítico.

gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)
## [1] 7.814728

Gráfico de la prueba Farrer-Glaubar.

library(fastGraph)
chi_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
shadeDist(chi_FG, ddist = "dchisq", parm1 = gl, lower.tail = FALSE, xmin = 0, xlab = "VC: 7.81   FG: 31.38",
          main = "Gráfico de la prueba de FG")
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)

Cálculo de FG usando la librería “mctest”

#Utilizando la libreria "mctest"
library(mctest)
mctest::omcdiag(mod = Estimación_del_modelo)
## 
## Call:
## mctest::omcdiag(mod = Estimación_del_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 test

Cálculo de FG usando la librería “psych”

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test) 
## $chisq
## [1] 31.38122
## 
## $p.value
## [1] 7.065806e-07
## 
## $df
## [1] 3

Gráfico de la prueba Farrer-Glaubar.

library(fastGraph)
alpha_sig<-0.05
VC<-qchisq( p = alpha_sig,df= FG_test$df,lower.tail = FALSE)
shadeDist(chi_FG,ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE, xmin = 0,
          sub=paste("VC:",round(VC,2),"","FG:",round(FG_test$chisq,2)),
          main = "Gráfico de la prueba de FG")

Interpretación: En base a la regla de decisión se puede concluir que: se rechaza la H0, esto debido a que \(\chi^2_{FG} \geq V.C.\), además se puede decir que hay evidencia de colinealidad en los regresores.

Factores Inflacionarios de la Varianza (FIV).

Los denominados, variance inflation factor (VIF), por sus siglas en inglés, determinan el tamaño relativo de la varianza del j-ésimo parámetro estimado, respecto a la varianza esperada del estimador en ausencia de colinealidad.

  • Estimador de la varianza para el j-ésimo parámetro:

\[\text{Var}(\beta_j) = \frac{\sigma^2_\epsilon}{(n-1)\cdot\text{Var}(X_j)}\cdot \frac{1}{1 - R^2_j}\]

  • Estimador mínimo de la varianza para el j-ésimo parámetro (en ausencia de colinealidad \(R^2_{j}\) = 0

\[\text{Var}(\beta_j)_{\text{min}} =\frac{\sigma^2_\epsilon}{(n-1)\cdot\text{Var}(X_j)}\]

Fórmula para el cálculo de los VIF.

\[\text{VIF}_j = \frac{\text{Var}(\beta_j)}{\text{Var}(\beta_j)_{\text{mín}}} = \frac{1}{1 - R^2_j}\]

Donde \(R^2_{j}\), es el coeficiente de determinación de la regresión de \(X_{j}\) contra el resto de regresores.

Matricialmente los VIF, se obtienen de la diagonal principal de la inversa de la matriz de Correlación:

\[\text{VIF}_{j = 1,2,...,k-1} = \text{diag}(R^{-1})\]

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  10

Cálculo manual:

Matriz de Correlación de los regresores del modelo (Como se obtuvo con anterioridad):

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^{-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

VIF’s para el modelo estimado:

VIFs<-diag(inversa_R)
print(VIFs)
##  lotsize    sqrft    bdrms 
## 1.037211 1.418654 1.396663

Cálculo de los VIF’s usando “car”

library(car)
VIFs_car<-vif(Estimación_del_modelo)
print(VIFs_car)
##  lotsize    sqrft    bdrms 
## 1.037211 1.418654 1.396663

Cálculo de los VIF’s usando “mctest”

library(mctest)
mc.plot(mod = Estimación_del_modelo,vif = 2)