Ejercicio de multicolinealidad

data frame

library(wooldridge)
data(hprice1)
head(force(hprice1),n=5) #mostras 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

estime el siguiente modelo : price = ˆα + ˆα1(lotsize) + ˆα2(sqrft) + ˆα3(bdrms) + ε

library(stargazer)
modelo_estimado<-lm(price~lotsize+sqrft+bdrms,data = hprice1)
stargazer(modelo_estimado,type = "html",little = "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
R²	0.672
Adjusted R²	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

modelo estimado

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

Indíce 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.

índice 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 $κ (x)$ , es igual a la raíz cuadrada de la razón entre la raíz característica más grande $λ_{m a x}$ , y la raíz característica más pequeña $λ_{m i n}$ , de la matríz $X^{t} X$ , normalizada; es decir: $κ (x) = \sqrt{\frac{λ_{m a x}}{λ_{m i n}}}$

interpretación

Sí $κ (x)$ , es inferior o igual a 20, la multicolinealidad es leve, ya que no se considera un problema.

Para 20 < $κ (x)$ < 30, la multicolinealidad se considera moderada.

En el caso de que $κ (x)$ ≥ 30, la multicolinealidad es severa.

Cálculo “Manual”:

Se basa la matríz $X^{t} X$ :

library(stargazer)
X_mat<-model.matrix(modelo_estimado)
stargazer(head(X_mat,n=6),type = "html")


	(Intercept)	assess	bdrms	lotsize	colonial	llotsize

1	1	349.100	4	6,126	1	8.720
2	1	351.500	3	9,903	1	9.201
3	1	217.700	3	5,200	0	8.556
4	1	231.800	3	4,600	1	8.434
5	1	319.100	4	6,095	1	8.715
6	1	414.500	5	8,566	1	9.056

XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "html")


	(Intercept)	assess	bdrms	lotsize	colonial	llotsize

(Intercept)	88	27,784.800	314	793,748	61	783.649
assess	27,784.800	9,563,053.000	102,507.500	278,300,049.000	19,578.900	250,005.600
bdrms	314	102,507.500	1,182	2,933,767	228	2,802.953
lotsize	793,748	278,300,049.000	2,933,767	16,165,159,010	555,967	7,457,452.000
colonial	61	19,578.900	228	555,967	61	544.060
llotsize	783.649	250,005.600	2,802.953	7,457,452.000	544.060	7,004.230

La matriz $X^{t} X$ debe ser normalizada, para evitar sesgo de la escala de las variables.

Cálculo de la matriz de normalización: $S_{n} = {\begin{matrix} \frac{1}{\sqrt{{X^{t} X}_{i j}}} & ; i = j \\ 0 & ; i \neq j \end{matrix}}$

library(stargazer)
options(scipen = 999)
Sn<- solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "html")


0.107	0	0	0	0	0
0	0.0003	0	0	0	0
0	0	0.029	0	0	0
0	0	0	0.00001	0	0
0	0	0	0	0.128	0
0	0	0	0	0	0.012

$X^{t} X$ normalizada

library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "html",digits = 4)


1	0.9578	0.9736	0.6655	0.8326	0.9982
0.9578	1	0.9642	0.7078	0.8106	0.9660
0.9736	0.9642	1	0.6712	0.8491	0.9742
0.6655	0.7078	0.6712	1	0.5599	0.7008
0.8326	0.8106	0.8491	0.5599	1	0.8323
0.9982	0.9660	0.9742	0.7008	0.8323	1

Autovalores de $X^{t} X$ normalizada:

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


5.193	0.492	0.237	0.049	0.028	0.0005

Cálculo de $κ (x) = \sqrt{\frac{λ_{m a x}}{λ_{m i n}}}$ :

K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)

## [1] 106.4903

Debido a la interpretación antes vista, cocluimos de que la multicolinealidad no se considera un problema, por lo siguiente: $κ (x)$ < 20.

Cálculo del Indíce de Condición usando librería “mctest”.

library(mctest)
X_mat<-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.1420         0
## Farrar Chi-Square:       164.9525         1
## Red Indicator:             0.3832         0
## Sum of Lambda Inverse:    12.3289         0
## Theil's Method:           -0.8940         0
## Condition Number:        106.4903         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

library(olsrr)
ols_eigen_cindex(model = modelo_estimado)

##     Eigenvalue Condition Index    intercept       assess       bdrms
## 1 5.1933823834        1.000000 3.563203e-05 0.0013514449 0.001229095
## 2 0.4923113019        3.247919 4.447199e-05 0.0001049883 0.001469387
## 3 0.2365049515        4.686030 2.860775e-04 0.0112620908 0.004101462
## 4 0.0490596083       10.288762 4.832476e-03 0.4576465408 0.014764598
## 5 0.0282837924       13.550531 1.198647e-03 0.2138228412 0.907363692
## 6 0.0004579625      106.490334 9.936027e-01 0.3158120940 0.071071766
##       lotsize    colonial     llotsize
## 1 0.003714445 0.007882944 3.054565e-05
## 2 0.296269130 0.061464755 1.170881e-05
## 3 0.032237901 0.854490636 2.042611e-04
## 4 0.010083947 0.002127499 2.652757e-03
## 5 0.009254094 0.069640172 1.673352e-03
## 6 0.648440482 0.004393993 9.954274e-01

Prueba de Farrar - Glaubar.

(Implementación de la prueba de Bartlett) Está prueba identifica si a nivel poblacional, los regresores del modelo presentan independencia estadística (son ortogonales), a través de la matriz de correlación muestral R y se verifica si a nivel poblacional dicha matriz identidad, la hipótesis de la prueba son las siguientes: ${\begin{matrix} H_{0} : & R \sim I \\ H_{1} : & R ≁ I \end{matrix}$

Si no se rechaza $H_{0}$ , no hay evidencia de multicolinealidad, caso contrario;

Si se rechaza $H_{0}$ hay evidencia de multicolinealidad El contraste se realiza a través de: $χ_{F G}^{2} = - (n - 1 - \frac{2 m + 5}{6}) \ln (| R |)$ con $g l = m (m - 1) / 2$

Rechazar; $χ_{F G}^{2}$ ≥ V.C., ó si p≤α

Cálculo “Manual”.

Cáñculo de $| R |$

library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "html")


	assess	bdrms	lotsize	colonial	llotsize

1	0.350	0.513	-0.284	0.662	-0.340
2	0.375	-0.675	0.087	0.662	0.543
3	-1.029	-0.675	-0.375	-1.495	-0.641
4	-0.881	-0.675	-0.434	0.662	-0.866
5	0.035	0.513	-0.287	0.662	-0.349
6	1.036	1.702	-0.045	0.662	0.277

Calcular la matriz R.

library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "html",digits = 4)


	assess	bdrms	lotsize	colonial	llotsize

assess	1	0.4825	0.3281	0.0829	0.5717
bdrms	0.4825	1	0.1363	0.3046	0.1695
lotsize	0.3281	0.1363	1	0.0140	0.8079
colonial	0.0829	0.3046	0.0140	1	0.0386
llotsize	0.5717	0.1695	0.8079	0.0386	1

Calcular

R |

determinante_R<-det(R)
print(determinante_R)

## [1] 0.1419755

Aplicando la prueba de Farrer Glaubar (Bartlett).

Estadístico $χ_{F G}^{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] 164.9525

Valor Crítico.

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

## [1] 18.30704

Regla de desición.

En base a la Regla de decisión concluimos con lo siguiente: $χ_{F G}^{2} \geq V . C .$ ≥V.C. se rechaza $H_{0}$ , 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.1420         0
## Farrar Chi-Square:       164.9525         1
## Red Indicator:             0.3832         0
## Sum of Lambda Inverse:    12.3289         0
## Theil's Method:           -0.8940         0
## Condition Number:        106.4903         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

Gráfica.

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

## $chisq
## [1] 164.9525
## 
## $p.value
## [1] 3.072151e-30
## 
## $df
## [1] 10

Gráfica.

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.

Factores Inflacionarios de la Varianza (FIV).

Los denominados, variantes inflation factor (VIF), por sus siglas en inglés, determinante 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: $V a r (β_{j}) = \frac{σ_{ε}^{2}}{(n - 1) . V a r (X_{j})} . \frac{1}{1 - R_{j}^{2}}$

Estimador mínimo de la varianza para el j-ésimo parámetro (en ausencia de colinealidad $R_{j}^{2} = 0$ $V a r (β_{j})_{m í n} = \frac{σ_{ε}^{2}}{(n - 1) . V a r (X_{j})}$

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

$V I F_{j} = \frac{V a r (β_{j})}{V a r (β_{j})_{m í n}} = \frac{1}{1 - R_{j}^{2}}$

Donde $R_{j}^{2}$ , es el coeficiente de determinación de la regresión de $X_{j}$

Marticialmente los VIF, se obtienen de la diagonal principal de la inversa de la matriz de Correlación: $V I F_{j = 1, 2, . . ., k - 1} = d i a g (R^{- 1})$

Referencia entre.

$R_{j}^{2}$

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)

##              assess     bdrms    lotsize   colonial  llotsize
## assess   1.00000000 0.4824739 0.32814633 0.08293582 0.5716654
## bdrms    0.48247394 1.0000000 0.13632563 0.30457549 0.1694902
## lotsize  0.32814633 0.1363256 1.00000000 0.01401865 0.8078552
## colonial 0.08293582 0.3045755 0.01401865 1.00000000 0.0386421
## llotsize 0.57166539 0.1694902 0.80785523 0.03864210 1.0000000

Inversa de la matriz de correlación $R^{- 1}$

inversa_R<-solve(R)
print(inversa_R)

##              assess      bdrms    lotsize   colonial   llotsize
## assess    2.1535576 -0.9010888  0.8347216  0.1520968 -1.7586001
## bdrms    -0.9010888  1.5104833 -0.3640744 -0.4021983  0.5687703
## lotsize   0.8347216 -0.3640744  3.2049651  0.1130042 -3.0089889
## colonial  0.1520968 -0.4021983  0.1130042  1.1142184 -0.1531266
## llotsize -1.7586001  0.5687703 -3.0089889 -0.1531266  4.3456744

VIF’s para el modelo estimado:

VIFs<-diag(inversa_R)
print(VIFs)

##   assess    bdrms  lotsize colonial llotsize 
## 2.153558 1.510483 3.204965 1.114218 4.345674

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)

## Variable `Component` is not in your data frame :/

Cálculo de los VIF´s usando “car”.

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:psych':
## 
##     logit

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)