Ejercicio multicolinealidad
FERNANDO JAVIER ASCENCIO ORELLANA
2023-05-12
Cargar los datos
library(wooldridge)
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.829630Estimar el modelo
library(stargazer)
modelo_estimado <- lm(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.01Indice de condición
Cálculo manual
library(stargazer)
x_mat <- model.matrix(modelo_estimado)
stargazer(head(x_mat,n=6),type="text")##
## =================================
## (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_mat <- t(x_mat)%*%x_mat
stargazer(xx_mat,type = "text")##
## ==============================================================
## (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
## --------------------------------------------------------------library(stargazer)
options(scipen = 999999)
sn <- solve(diag(sqrt(diag(xx_mat))))
stargazer(sn,type = "text") ##
## ==========================
## 0.107 0 0 0
## 0 0.00001 0 0
## 0 0 0.0001 0
## 0 0 0 0.029
## --------------------------\(X^tX\) normalizada:
library(stargazer)
xx_norm <- (sn%*%xx_mat)%*%sn
stargazer(xx_norm,type = "text", 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 \(X^tX\) normalizada:
library(stargazer)
lambdas <- eigen(xx_norm,symmetric = TRUE)
stargazer(lambdas$values, type = "text")##
## =======================
## 3.482 0.455 0.039 0.025
## -----------------------Cálculo de \(\kappa(x)=\sqrt{\lambda_{max}}/{\lambda_{min}}\)
K <- sqrt(max(lambdas$values)/min(lambdas$values))
print(K)## [1] 11.86778Como \(\kappa(x)< 30\) se considera que la multicolinealidad es leve.
Farrar-Glaubar
Cálculo manual
Cálculo de \(|R|\)
library(stargazer)
zn <- scale(x_mat[,-1])
stargazer(head(zn,n=6), type = "text")##
## =======================
## 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
## -----------------------Cáclulo la matriz R
library(stargazer)
n <- nrow(zn)
R <- (t(zn)%*%zn)*(1/(n-1))
stargazer(R,type = "text", 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.6917931Prueba de Farrer-Glaubar (Bartlett)
Estadístico \(\chi_{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.38122Valor crítico
gl <- m*(m-1)/2
vc <- qchisq(p=0.95,df=gl)
print(vc)## [1] 7.814728Regla de decisión
library(fastGraph)
shadeDist(qchisq(0.1, chi_FG, lower.tail = FALSE),
ddist = "dchisq",
parm1 = chi_FG,
parm2 = vc,
lower.tail = FALSE,
col = c("black", "red"),
sub=paste("Valor crítico:",round(vc,2)," ","FB:",round(chi_FG,2)))
como \(\chi_{FG}^2 \geq V.C\). Se rechaza la \(H_0\); es decir, hay evidencia de colinealidad en los regresores.
Factores inflacionarios de la Varianza (FIV)
Relación entre \(R_{j}^2\)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionR_cuadrado <- c(0,0.5,.8,.9)
as.data.frame(R_cuadrado) %>% mutate(VIF=1/(1-R_cuadrado))## R_cuadrado VIF
## 1 0.0 1
## 2 0.5 2
## 3 0.8 5
## 4 0.9 10Cálculo manual
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.0000000Inversa 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.39666321VIF´s para el modelo estimado
VIFs <- diag(inversa.R)
print(VIFs)## lotsize sqrft bdrms
## 1.037211 1.418654 1.396663plot(VIFs)
library(mctest)
mc.plot(mod=modelo_estimado,vif = 2)