Practica de colinealidad

Se presenta a continuación el desarrollo de la practica, pruebas de colinealidad a partir de la siguiente información.

Importación de datos

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.829630

Desarrollo de la practica.

1. Estimación del modelo.

Se estima el modelo a partir de los datos y se presenta en formato Apa utilizando la libreria Stargazer

\(price = ˆα + ˆα1(lotsize) + ˆα2(sqrft) + ˆα3(bdrms) + 𝑒\)

options(scipen = 999999)
library(wooldridge)
library(stargazer)
data(hprice1)
modelo_estimado<-lm(price~lotsize+sqrft+bdrms,data=hprice1)
stargazer(modelo_estimado, title = "Modelo estimado", type = "html")
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. Verificación de evidencia de no colinealidad.

Se verifica la evidencia de colinealidad menciante el indice de condición, prueba FG y Factores inflacionarios de la varianza.

a) Indice de condición y prueba de FG.

Indice de condición.

En primer lugar, el número que brinda el indice de condición mide la sensibilidad de las estimaciones mínimo cuadráticas ante pequeños cambios en los datos. El indice se representa de la siguiente manera:

\[k(x)= √ 𝛌max/𝛌min\]

Se interpreta de la siguiente manera:

  • Sí \(k(x) ≤ 20\), la multicolinealidad es leve.

  • Para \(20< k(x) <30\), la multicolinealidad se considera moderada

  • En el caso de que \(K(x)≥30\) la multicolinealidad es severa.

Calculo manual.

Calculo de la sigma matriz \(X^tX\)

library(stargazer)
X_mat<-model.matrix(modelo_estimado)
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

Posteriomente se normaliza la matriz, con la finalidad de evitar sesgos de escala.

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

Matriz 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

Cálculo de autovalores

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

Finalmente se realiza el cálculo del indice.\(k(x)= √ 𝛌max/𝛌min\)

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

Interpretación: \[11.86778 < 20\] De acuerdo a la interpretación inicial, el indice obtenido es menor a 30 y menor a 20, por lo cual indica que la multicolieanidad es leve y no se considera un problema.

Utilizando libreria “mctest”

library(mctest)
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 test

Interpretación:Los resultados aplicando la libreria mctest, indica que existe colinealidad, el numero de condición es de 11.8678. De acuerdo con la interpretación del indice K, 11.8670 es menos que 20 lo cual confirma que existe colinealidad a nivel leve y no se considera un problema. Este resultado coincide con el obtenido mediante el calculo manual.

Utilizando libreria “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.975025872

Interpretacion: De igual forma, el resultado de la segunda columna y ultima fila, es de 11.867781 que coincide con las pruebas anteriores.

Prueba de Farrar - Glaubar

En general, la prueba FG identifica si a nivel poblacional, existe presencia de independencia estadistica en los regresores del modelo. Mediante la matriz de correlación muestral R, ademas se verifica si a nivel poblacional dicha matriz de correlacion corresponde a una matriz identidad, para ello las hipotesis de la prueba son las siguientes:

\[H0:R∼I\] \[H1:R≁I\]

  • Si no se rechaza H0, no hay evidencia de multicolinealidad, caso contrario
  • Si se rechaza H0, hay evidencia de multicolinealidad.

Regla de desición: Rechazar {H_0} si \(X ^2_{FG} ≥ V.C\) , o si \(p ≤ α\)

Cálculo manual.

Calculo 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))
#También se puede calcular R a través de cor(X_mat[,-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)

Estadistico \(X_{FG}^2\)

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 critico:

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

Regla de desición: Rechazar {H_0} si \(X ^2_{FG} ≥ V.C\) , o si \(p ≤ α\)

Interpretación: El valor crítico obtenido de forma manual es de 7.814728 y el estadistico \(X ^2_{FG}\) es de 31.38122. En conclusión, el estadistico es mayor que el valor critico, existe evidencia de colinealidad en los regresores.

\(X ^2_{FG} > V.C\) \(31.38122 > 7.814728\)

Cálculo de FG mediante libreria “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 test

Interpretación: El valor del estadistico es de 31.3812 y el indice de condición es de 11.8678, ambos valores coinciden con las salidas anteriores. Además, segun este resultado la prueba si detecta colinealidad.

Cálculo de FG mediante libreria “psych”

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
## R was not square, finding R from data
print(FG_test)
## $chisq
## [1] 31.38122
## 
## $p.value
## [1] 0.0000007065806
## 
## $df
## [1] 3

Interpretación: Se aplica: Rechazar \({H_0}\) si \(p ≤ α\) El valor del estadistico coincide, sin embargo la salida proporciona el valor P- Value, el cual es de 0.0000007065806. Al aplicar la regla desición se concluye que el valor P-value es mucho menor que el valor de significancia .05 y por lo tanto se rechaza \({H_0}\)

Gráfica con libreria FastGraph.

# Calcular el valor crítico
vc <- qchisq(p = 0.05, df = FG_test$df, lower.tail = FALSE)

# Graficar la distribución chi-cuadrado
shadeDist(xshade = FG_test$chisq,
          ddist = "dchisq",
          parm1 = FG_test$df,
          main = "Prueba FG",
          sub = paste("vc:", vc, ", FG:", FG_test$chisq),
          lower.tail = FALSE)

# Agregar texto para indicar la región de rechazo
text(x = vc, y = 0.02, labels = "Zona de Rechazo", pos = 4, col = "red")
# Agregar línea vertical en el límite de la región de rechazo
abline(v = vc, col = "red", lty = 2)

b) Factores inflacionarios de la varianza.

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:

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 “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)

Interpretación: La variable Lotsize es la que presente mayor colinealidad, sin embargo para el valor tomado se considera que la colinealidad es minima y no genera problemas. Por lo tanto, las pruebas muestran evidencia de correlación estadistica.Sin embargo, al análizar el índice de condición y los VIF´s se da a comprender que la correlación no es tan significativa.

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

Cálculo con “car”

library(car)
VIFs_car<-vif(modelo_estimado)

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 = modelo_estimado,vif = 2)