Ejercicio de Multicolinealidad, Tarea Individual.

Carga de datos:

Utilizando los datos del dataframe hprice1: disponible en el paquete wooldridge use el siguiente código para generar el dataframe:

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

1. Modelo Estimado

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

#Sigma matriz
library(stargazer)
modelo_estimado<-lm(formula = price~lotsize+sqrft+bdrms,data = hprice1)
stargazer(modelo_estimado,type = "text",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:

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

Indice de Condicion de forma manual

Basado en la matriz \(X^{t}X\) :

Matriz X

#Matriz X
options(scipen = 999999)
Matriz_X <- model.matrix(modelo_estimado)
stargazer(head(Matriz_X, n = 6), type = "text",title = "Matriz X")

## 
## Matriz X
## =================================
##   (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  
## ---------------------------------

Matriz XX

#Matriz XX
Matriz_XX <- t(Matriz_X)%*%Matriz_X
stargazer(Matriz_XX, type = "text",title = "Matriz XX")

## 
## Matriz XX
## ==============================================================
##             (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  
## --------------------------------------------------------------

Cálculo de la matriz de normalización:

#Normalización de la matriz XX
library(stargazer)
options(scipen = 99999)
Sn <- solve(diag(sqrt(diag(Matriz_XX))))
stargazer(Sn, type = "text",title = "Matriz de normalización")

## 
## Matriz de normalización
## ==========================
## 0.107    0      0      0  
## 0     0.00001   0      0  
## 0        0    0.0001   0  
## 0        0      0    0.029
## --------------------------

\(X^{t}X\) normalizada

# Sigma matriz normalizada
library(stargazer)
options(scipen = 99999)
normalizada_XX<-(Sn%*%Matriz_XX)%*%Sn
stargazer(normalizada_XX,type = "text",title = "Matriz Normalizada")

## 
## Matriz Normalizada
## =======================
## 1     0.666 0.962 0.974
## 0.666   1   0.678 0.671
## 0.962 0.678   1   0.970
## 0.974 0.671 0.970   1  
## -----------------------

Autovalores de \(X^{t}X\) Normalizada:

#Autovalores de la matriz XX normalizada
lambdas<-eigen(normalizada_XX,symmetric = TRUE)
stargazer(lambdas$values,type = "text",title = "Autovalores de matriz XX normalizada")

## 
## Autovalores de matriz XX normalizada
## =======================
## 3.482 0.455 0.039 0.025
## -----------------------

Cálculo de \(k(x)=√\frac{\gamma max}{\gamma min}\)

#Cálculo de K (número de condición)
K<-sqrt(max(lambdas$values)/min(lambdas$values))
stargazer(K,title = "K (número de condición)",type = "text")

## 
## K (número de condición)
## ======
## 11.868
## ------

Como k(x)=<20, la multicolinealidad se considera que no es un problema, ya que es leve.

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

library(mctest)
matriz_X<-model.matrix(modelo_estimado)
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

Cálculo del Indice de Condición usando la librería “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

El índice de condición es de 11.8678, por lo que existe evidencia que los regresores presentan leve multicolinealidad.

Prueba de Farrar-Glaubar

Cálculo “manual”:

Calculo de |R|

#Cálculo de R
Zn<-scale(Matriz_X[,-1])
stargazer(head(Zn,n=6),type = "text",title = "Cálculo de R")

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

#Calcular la matriz R 
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "text",title = "Calcular la matriz R",digits = 4)

## 
## Calcular la matriz R
## =============================
##         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 χ2FG

m<-ncol(matriz_X[,-1])
n<-nrow(matriz_X[,-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

Como χ2FG ≥ V.C se rechaza H0, por lo tanto hay evidencia de colinealidad en los regresores.

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

Grafica

library(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(Matriz_X[,-1])
print(FG_test)

## $chisq
## [1] 31.38122
## 
## $p.value
## [1] 0.0000007065806
## 
## $df
## [1] 3

Grafica

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

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)

B) Factores inflacionarios de la varianza, presente sus resultados de forma tabular y de forma gráfica

Referencia entre \(R\frac{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:

R_inversa <-solve(R)
print(R_inversa)

##             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(R_inversa)
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)
print(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_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)

Interpretación: Los Valores Inflacionarios de la Varianza de los regresores son menores que 2, existe evidencia de multicolinealidad leve.