Este texto se basa en los siguientes materiales:

> library(tidyverse)
> library(broom)
> library(car)

Colinealidad entre regresores

Pensemos brevemente en dos de las variables de nuestro modelo: t_hogar y nivel_ed_agg. ¿Hay correlación entre ambas? Es una pregunta releveante.

Esto expresa un problema potencial bastante común en la regresión múltiple: la correlación entre las variables predictoras. Decimos que las dos variables predictoras son colineales (pronunciadas como colineales) cuando están correlacionadas, y esta multicolinealidad complica la estimación del modelo. Si bien es imposible evitar que surja la multicolinealidad en los datos de observación, los experimentos generalmente se diseñan para evitar que los predictores sean multicolineales.

Volvamos brevemente a nuestro ejemplo antropométrico para ilustrar esta situación de forma más clara.

> df <- read_delim("https://raw.githubusercontent.com/rmcelreath/rethinking/master/data/Howell1.csv",
+     delim = ";")
> 
> df <- df %>%
+     mutate(male = as.factor(case_when(male == 0 ~ "No", TRUE ~ "Yes")))

Filtremos los menores de 18 años:

> df_menores <- df %>%
+     filter(age < 18)

Si recordamos bien, la altura y la edad (y el peso) estaban fuertemente correlacionadas

> df_menores %>%
+     ggplot(aes(x = age, height)) + geom_point() + theme_minimal()

En estos casos, cuando dos regresores están altamente correlacionados, utilizamos el término colinealidad. La presencia de colinealidad puede plantear problemas en el contexto de regresión, ya que puede ser difícil separar los efectos individuales de las variables colineales en la respuesta. En otras palabras, dado que la edad y la altura tienden a aumentar o disminuir juntas, puede ser difícil determinar cómo cada uno por separado se asocia con la variable dependiente (peso).

Corramos una regresión con todos los predictores. Es decir, vamos a tratar de modelar el peso en función de la edad, la altura y el sexo.

Un problema de la existencia de colinealidad es que hace que las inferencias de los parámetros (\(\beta\)) sean menos precisas. Ante la presencia de la colinealidad, los errores estándar de las estimaciones se incrementan:

Entrenemos una regresión con las tres variables independientes:

> model1 <- df_menores %>%
+     lm(weight ~ ., data = .)
> tidy(model1)
## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -11.4      1.82       -6.29 2.22e- 9
## 2 height         0.243    0.0247      9.81 1.34e-18
## 3 age            0.431    0.118       3.64 3.53e- 4
## 4 maleYes        0.538    0.422       1.28 2.04e- 1

Ahora, eliminenos la edad:

> model2 <- df_menores %>%
+     lm(weight ~ . - age, data = .)
> tidy(model2)
## # A tibble: 3 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -17.2     0.935     -18.4   5.86e-44
## 2 height         0.328   0.00831    39.4   3.75e-93
## 3 maleYes        0.239   0.427       0.559 5.77e- 1

Noten como los errores estándar (que vamos a definir con precisión más adelante pero que por ahora podemos asociarlos a la incerteza muestral sobre los parámetros) de height cambian fuertemente en ambos modelos: en el primer modelo es de 0.025; en el segundo, es 0.008.

Para evitar tal situación, es deseable identificar y abordar posibles problemas de colinealidad al ajustar el modelo. Una forma sencilla de detectar la colinealidad es observar la matriz de correlación de los predictores. Una celda de esta matriz con valores grandes indica un par de variables altamente correlacionadas.

Desafortunadamente, no todos los problemas de colinealidad pueden ser detectado por inspección de la matriz de correlación: es posible que exista colinealidad entre tres o más variables incluso si no hay un par de variables tiene una correlación particularmente alta. A esta situación la llamamos multicolinealidad. En lugar de inspeccionar la matriz de correlación, una mejor manera de evaluar la multicolinealidad es calcular el factor de inflación de la varianza (VIF, por sus siglas en inglés).

El VIF es el cociente de la varianza de \(\beta_{j}\) al ajustar el modelo completo dividido por el varianza de \(\beta_{j}\) si se ajusta por sí solo. El valor más pequeño posible para VIF es 1, lo que indica la ausencia total de colinealidad. Típicamente en la práctica hay una pequeña cantidad de colinealidad entre los predictores. Como regla aproximada, un valor VIF que excede 5 o 10 indica una cantidad problemática de colinealidad.

El VIF para cada variable se puede calcular usando la fórmula

\[VIF(\beta_{j}) = \frac{1}{1-R^2_{X_{j}|X_{-j}}}\]

donde \(R^2_{X_{j}|X_{-j}}\) es el \(R^2\) de una regresión de \(X_{j}\) contra todos los otros predictores. Si \(R^2_{X_{j}|X_{-j}}\) es cercano a 1, entonces, hay colinealidad, por lo cual \(VIF\) será elevado.

> vif(model1)
##   height      age     male 
## 9.473102 9.437404 1.043606

Vemos cómo height y age tienen valores bastante cercanos a nuestro límite.

Ante el problema de la colinealidad, existen dos soluciones sencillas. La primera es eliminar una de las variables problemáticas de la regresión. Esto generalmente se puede hacer sin mucho compromiso con la regresión. ajuste, ya que la presencia de colinealidad implica que la información que este variable proporciona sobre la respuesta es redundante en presencia de la otras variables.

La segunda solución puede ser combinar las varaibles colineales en un solo predictor. Una opción es generar alguna tipología o algún índice. Veremos en el módulo 3 algunas formas de reducción de dimensionalidad (Análisis de Componentes Principales -PCA o Análisis de Correspondencias múltiples - MCA) son formas habituales de lidiar con la multicolinealidad.

¿Qué pasa en nuestro modelo de determinación de ingresos? ¿Qué variables podrían tener algún grado de colinealidad?

> enes <- read_rds('./data/ENES_Personas_M1_EOW.rds')
> 
> enes <- enes %>% 
+    mutate(v109 = case_when(
+                   v109=='Varón' ~ 'Masculino',
+                   TRUE ~ 'No masculino'),
+           nivel_ed_agg = case_when(
+                   nivel_ed == 'Menores de 5 años' | 
+                   nivel_ed == 'Sin instrucción (incluye nunca asistió o sólo asistió a sala de 5)' |
+                   nivel_ed == 'Primaria/EGB incompleto' | 
+                   nivel_ed == 'Primaria/EGB completo' |
+                   nivel_ed == 'Educación especial' | nivel_ed == 'NS/NR'~ '0_Bajo',
+      
+                 nivel_ed == 'Secundario/Polimodal incompleto' | 
+                 nivel_ed == 'Secundario/Polimodal completo' ~ '1_Medio',
+     
+                 nivel_ed == 'Terciario incompleto' | 
+                 nivel_ed == 'Terciario completo' |
+                 nivel_ed == 'Universitario incompleto' | 
+                 nivel_ed == 'Universitario completo' ~ '2_Alto'
+    ),
+    class_eow_agg = case_when(
+          class_eow == 'Managers' | class_eow == 'Supervisores' ~  'Managers/superv.',
+          class_eow == 'Trabajadores' ~ 'Trabajadores',
+          class_eow == 'Pequeña burguesía' ~ 'Pequeña burguesía',
+          class_eow == 'Inactivo, desocupado o menor' ~ 'Inactivo, desocupado o menor',
+          class_eow == 'Empleadores' ~ 'Empleadores'
+   )
+    )

Estimemos la regresión:

> lm_enes <- enes %>%
+     filter(estado == "Ocupado") %>%
+     lm(v213b ~ t_hogar + v108 + v109 + nivel_ed_agg + class_eow_agg, data = .)

Y calculemos el VIF:

> vif(lm_enes)
##                   GVIF Df GVIF^(1/(2*Df))
## t_hogar       1.089919  1        1.043992
## v108          1.161298  1        1.077635
## v109          1.063069  1        1.031052
## nivel_ed_agg  1.190976  2        1.044662
## class_eow_agg 1.088522  3        1.014237

Pese a haber cierto grado de colinealidad (esperable, por cierto) ninguno supera el valor límite de 5.