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

1. Estime el modelo

library(stargazer)
options(scipen = 999999)

modelo_hprice1<-lm(formula = price~bdrms+lotsize+sqrft,data = hprice1)

stargazer(modelo_hprice1,title = "Modelo hprice1",type = "text",digits = 8)
## 
## Modelo hprice1
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                price           
## -----------------------------------------------
## bdrms                       13.85252000        
##                            (9.01014500)        
##                                                
## lotsize                    0.00206771***       
##                            (0.00064213)        
##                                                
## sqrft                      0.12277820***       
##                            (0.01323741)        
##                                                
## Constant                   -21.77031000        
##                            (29.47504000)       
##                                                
## -----------------------------------------------
## Observations                    88             
## R2                          0.67236220         
## Adjusted R2                 0.66066090         
## Residual Std. Error    59.83348000 (df = 84)   
## F Statistic         57.46023000*** (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.

#Calcular matriz X'X
library(stargazer)
X_mat<-model.matrix(modelo_hprice1)
stargazer(head(X_mat,n=6),type= "text")
## 
## =================================
##   (Intercept) bdrms lotsize sqrft
## ---------------------------------
## 1      1        4    6,126  2,438
## 2      1        3    9,903  2,076
## 3      1        3    5,200  1,374
## 4      1        3    4,600  1,448
## 5      1        4    6,095  2,514
## 6      1        5    8,566  2,754
## ---------------------------------
XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "text")
## 
## ==============================================================
##             (Intercept)   bdrms      lotsize         sqrft    
## --------------------------------------------------------------
## (Intercept)     88         314       793,748        177,205   
## bdrms           314       1,182     2,933,767       654,755   
## lotsize       793,748   2,933,767 16,165,159,010 1,692,290,257
## sqrft         177,205    654,755  1,692,290,257   385,820,561 
## --------------------------------------------------------------
#Cálculo de la matriz de normalizacion
library(stargazer)
options(scipen = 999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "text")
## 
## ==========================
## 0.107   0      0      0   
## 0     0.029    0      0   
## 0       0   0.00001   0   
## 0       0      0    0.0001
## --------------------------
#X'X Normalizada
library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "text",digits = 4)
## 
## ===========================
## 1      0.9736 0.6655 0.9617
## 0.9736   1    0.6712 0.9696
## 0.6655 0.6712   1    0.6776
## 0.9617 0.9696 0.6776   1   
## ---------------------------

Cálculo del indice de condición

#Autovalores de X'X Normalizada
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")
## 
## =======================
## 3.482 0.455 0.039 0.025
## -----------------------
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)
## [1] 11.86778

INTERPRETACIÓN: El indice de condición es inferior a 20, por tanto,se puede evidenciar que el modelo adolece de multicolinealidad leve, lo cual no representa un problema.

Uso de librería “mctest”

library(mctest)
eigprop(mod = modelo_hprice1)
## 
## Call:
## eigprop(mod = modelo_hprice1)
## 
##   Eigenvalues      CI (Intercept)  bdrms lotsize  sqrft
## 1      3.4816  1.0000      0.0037 0.0029  0.0278 0.0042
## 2      0.4552  2.7656      0.0068 0.0051  0.9671 0.0061
## 3      0.0385  9.5082      0.4726 0.0169  0.0051 0.8161
## 4      0.0247 11.8678      0.5170 0.9750  0.0000 0.1737
## 
## ===============================
## Row 4==> bdrms, proportion 0.975026 >= 0.50 
## Row 2==> lotsize, proportion 0.967080 >= 0.50 
## Row 3==> sqrft, proportion 0.816079 >= 0.50

Cálculo de R

#Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "text")
## 
## =======================
##   bdrms  lotsize sqrft 
## -----------------------
## 1 0.513  -0.284  0.735 
## 2 -0.675  0.087  0.108 
## 3 -0.675 -0.375  -1.108
## 4 -0.675 -0.434  -0.980
## 5 0.513  -0.287  0.867 
## 6 1.702  -0.045  1.283 
## -----------------------

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 = "text",digits = 4)
## 
## =============================
##         bdrms  lotsize sqrft 
## -----------------------------
## bdrms     1    0.1363  0.5315
## lotsize 0.1363    1    0.1838
## sqrft   0.5315 0.1838    1   
## -----------------------------

Calcular R

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

Aplicando la prueba de Farrer Glaubar (Bartlett)

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

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

Uso de librería psych

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

Graficando los resultados

library(fastGraph)
shadeDist(xshade = FG_test$chisq, ddist = "dchisq", parm1 = FG_test$df,
          lower.tail = FALSE, sub = paste("VC:", VC, "FG:", FG_test$chisq))

b) Factores inflacionarios de la varianza

#Matriz de Correlación de los regresores del modelo 
print(R)
##             bdrms   lotsize     sqrft
## bdrms   1.0000000 0.1363256 0.5314736
## lotsize 0.1363256 1.0000000 0.1838422
## sqrft   0.5314736 0.1838422 1.0000000

Inversa de la matriz de correlación

inversa_R<-solve(R)
print(inversa_R)
##               bdrms     lotsize      sqrft
## bdrms    1.39666321 -0.05582352 -0.7320270
## lotsize -0.05582352  1.03721145 -0.1610145
## sqrft   -0.73202696 -0.16101454  1.4186543

VIF’s para el modelo estimado

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

Obtención de los VIF’s, a través de la librería “car”

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

Obtención de los VIF’s, a través de la librería “mctest”

library(mctest)
mc.plot(mod = modelo_hprice1,vif = 2)