library(wooldridge)
data(hprice1)
head(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. Estimacion del Modelo.

library(stargazer)
modelo_estimado <- lm(price ~ lotsize + sqrft + bdrms, 
                      data = hprice1)

suppressMessages(
  suppressWarnings(
    stargazer(modelo_estimado,
              type = "text",
              title = "Modelo Estimado",
              digits = 3)
  )
)
## 
## 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:

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

Calculo Manual

library(knitr)
library(kableExtra)

X_mat <- model.matrix(modelo_estimado)

kable(head(X_mat, 6),
      caption = "Matriz de diseño (primeras observaciones)",
      digits = 3) %>%
  kable_styling(
    full_width = TRUE,   
    position = "right",   
    font_size = 16
  )
Matriz de diseño (primeras observaciones)
(Intercept) lotsize sqrft bdrms
1 6126 2438 4
1 9903 2076 3
1 5200 1374 3
1 4600 1448 3
1 6095 2514 4
1 8566 2754 5 
X_mat <- model.matrix(modelo_estimado)

XX_matrix <- t(X_mat) %*% X_mat
XX_matrix_fmt <- format(XX_matrix,
                        big.mark = ".",
                        decimal.mark = ",",
                        scientific = FALSE)
library(knitr)
library(kableExtra)

kable(XX_matrix_fmt,
      caption = "Matriz XtX") %>%
  kable_styling(
    full_width = TRUE,
    position = "right",
    font_size = 14
  )
Matriz XtX
(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 
options(scipen = 999)

library(stargazer)


options(scipen = 999)

X_mat <- model.matrix(modelo_estimado)
XX_matrix <- t(X_mat) %*% X_mat

Sn <- solve(diag(sqrt(diag(XX_matrix))))

suppressMessages(
  suppressWarnings(
    stargazer(Sn,
              type = "text",
              digits = 4,
              summary = FALSE)
  )
)
## 
## ============================
## 0.1066    0      0      0   
## 0      0.00001   0      0   
## 0         0    0.0001   0   
## 0         0      0    0.0291
## ----------------------------
library(stargazer)

options(scipen = 999)

X_mat <- model.matrix(modelo_estimado)
XX_matrix <- t(X_mat) %*% X_mat
Sn <- solve(diag(sqrt(diag(XX_matrix))))
XX_norm <- (Sn %*% XX_matrix) %*% Sn

stargazer(XX_norm,
          type = "html",
          digits = 4,
          summary = FALSE)
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 
library(stargazer)
lambdas <- eigen(XX_norm, symmetric = TRUE)

autovalores_row <- t(lambdas$values)
suppressMessages(
  suppressWarnings(
   stargazer(autovalores_row,
          type = "html",
          digits = 4,
          summary = FALSE)
  )
)
3.4816 0.4552 0.0385 0.0247 
K <- sqrt(max(lambdas$values) / min(lambdas$values))
K
## [1] 11.86778

Como κ(x)=11.86778 es menor a 20, se considera que la multicolinealidad es leve y no representa un problema en el modelo.

#Índice de Condición usando mctest

library(mctest)

X_mat <- model.matrix(modelo_estimado)

resultado_mctest <- mctest(mod = modelo_estimado)

resultado_mctest
## 
## 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

Índice de Condición usando olsrr

library(olsrr)

resultado_olsrr <- ols_eigen_cindex(model = modelo_estimado)

resultado_olsrr
##   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

Prueba de Farrar-Glaubar

Calculo Manual

library(knitr)
library(kableExtra)

Zn <- scale(X_mat[,-1])

kable(head(Zn, 6),
      caption = "Variables Estandarizadas (primeras observaciones)",
      digits = 3,
      row.names = TRUE) %>%
  kable_styling(full_width = TRUE, position = "right")
Variables Estandarizadas (primeras observaciones)
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(knitr)
library(kableExtra)
n <- nrow(Zn)

R <- cor(X_mat[,-1])
kable(R,
      caption = "Matriz de Correlación (R)",
      digits = 4,
      row.names = TRUE) %>%
  kable_styling(full_width = TRUE, position = "right")
Matriz de Correlación (R)
lotsize sqrft bdrms
lotsize 1.0000 0.1838 0.1363
sqrft 0.1838 1.0000 0.5315
bdrms 0.1363 0.5315 1.0000 
library(knitr)
determinante_R <- det(R)
kable(data.frame(Determinante_R = round(determinante_R, 6)),
      caption = "Determinante de la matriz de correlación (|R|)",
      align = "c")
Determinante de la matriz de correlación (|R|)
Determinante_R
0.691793 
library(knitr)

# Número de variables explicativas
m <- ncol(X_mat[,-1])

# Número de observaciones
n <- nrow(X_mat[,-1])

# Estadístico Chi-cuadrado de Farrar-Glauber
chi_FG <- -(n - 1 - (2*m + 5)/6) * log(determinante_R)

# Mostrar resultado
kable(data.frame(Chi_Cuadrado_FG = round(chi_FG, 4)),
      caption = "Estadístico Chi-cuadrado de Farrar-Glauber",
      align = "c")
Estadístico Chi-cuadrado de Farrar-Glauber
Chi_Cuadrado_FG
31.3812 
library(knitr)

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

kable(data.frame(
  Valor_Critico = round(VC, 4)
),
caption = "Valor Crítico de la distribución Chi-cuadrado",
align = "c")
Valor Crítico de la distribución Chi-cuadrado
Valor_Critico
7.8147

Gráfica prueba FG

library(fastGraph)

# Gráfica Chi-cuadrado
curve(dchisq(x, df = gl), from = 0, to = chi_FG + 20,
      main = "Prueba de Farrar-Glauber",
      xlab = "Chi-cuadrado",
      ylab = "Densidad")

# Estadístico
abline(v = chi_FG, col = "red", lwd = 2)

# Valor crítico
abline(v = VC, col = "green", lwd = 2)

Dado que χ²FG ≥ 7.8147, se rechaza la hipótesis nula (H₀). Por lo tanto, existe evidencia estadística de que los regresores del modelo (lotsize, sqrft y bdrms) presentan correlación entre sí, lo que indica la presencia de multicolinealidad.

Cálculo de FG usando “mctest”

library(mctest)

# Aplicación de la prueba
resultado_fg <- mctest::omcdiag(mod = modelo_estimado)

resultado_fg
## 
## 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

Cálculo de FG usando la “psych”

library(psych)

FG_test <- cortest.bartlett(X_mat[,-1])

FG_test
## $chisq
## [1] 31.38122
## 
## $p.value
## [1] 0.0000007065806
## 
## $df
## [1] 3
  1. Factores inflacionarios de la varianza, presente sus resultados de forma tabular y de forma gráfica ## Factores Inflacionarios de la Varianza (FIV)

Calculo Manual.

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
R.cuadrado.regresores <- c(
  summary(lm(lotsize ~ sqrft + bdrms, data = hprice1))$r.squared,
  summary(lm(sqrft ~ lotsize + bdrms, data = hprice1))$r.squared,
  summary(lm(bdrms ~ lotsize + sqrft, data = hprice1))$r.squared
)

as.data.frame(R.cuadrado.regresores) %>%
  mutate(VIF = 1/(1 - R.cuadrado.regresores))
##   R.cuadrado.regresores      VIF
## 1            0.03587644 1.037211
## 2            0.29510664 1.418654
## 3            0.28400778 1.396663
library(knitr)
library(kableExtra)
inversa_R <- solve(R)
kable(inversa_R,
      caption = "Inversa de la matriz de correlación (R⁻¹)",
      digits = 4,
      row.names = TRUE) %>%
  kable_styling(full_width = TRUE, position = "right")
Inversa de la matriz de correlación (R⁻¹)
lotsize sqrft bdrms
lotsize 1.0372 -0.1610 -0.0558
sqrft -0.1610 1.4187 -0.7320
bdrms -0.0558 -0.7320 1.3967

VIF a partir de la inversa de R

VIFs <- diag(inversa_R)

VIFs
##  lotsize    sqrft    bdrms 
## 1.037211 1.418654 1.396663

Cálculo de los VIF’s usando “performance”.

library(performance)
## Warning: package 'performance' was built under R version 4.5.3
VIFs <- multicollinearity(modelo_estimado, verbose = FALSE)

VIFs
## # Check for Multicollinearity
## 
## Low Correlation
## 
##     Term  VIF    VIF 95% CI adj. VIF 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]
library(see)
## Warning: package 'see' was built under R version 4.5.3
plot(VIFs)

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

library(car)
## Cargando paquete requerido: carData
## 
## Adjuntando el paquete: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:psych':
## 
##     logit
VIFs_car <- vif(modelo_estimado)

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)