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:

> df_menores %>%
+         select(age, height) %>%
+         cor(method = "pearson")
##             age   height
## age    1.000000 0.943313
## height 0.943313 1.000000

Vemos, entonces, que age y height presentan una alta correlación lineal (medida por el R de Perason).

Desafortunadamente, no todos los problemas de colinealidad pueden ser detectados por inspección de la matriz de correlación. Esta solo nos brinda informació sobre correlaciones bivariadas (pares de variables). Sin embargo, 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 Variance Inflation Factor).

El VIF es el cociente de la varianza de \(\beta_{j}\) al ajustar el modelo completo dividido por la 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 siguiente 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.

Conviene detenerse en qué es esa regresión auxiliar: se trata de hacer una regresión de \(X_{j}\) contra los demás predictores, sin usar la variable dependiente. Su \(R^2\) responde a la pregunta: ¿cuánto de la variación de \(X_{j}\) ya está contenida en el resto del modelo? Si la respuesta es “casi toda”, esa variable aporta poca información propia y su coeficiente queda mal identificado.

La lógica es la siguiente: se toma la variable \(X_{i}\) y se hace una regresión utilizándola como variable dependiente. Esta regresión intenta “explicar” su comportamiento en función de todas las demás variables independientes del modelo original. Así, si \(X_{i}\) no puede ser explicada por el resto de los predictores, entonces \(R^2_{X_{j}|X_{-j}}\) será cercano a 0 y el \(VIF\) será cercano a 1. Por el contrario, di \(R^2_{X_{j}|X_{-j}}\) es cercano a 1, entonces, hay colinealidad, por lo cual \(VIF\) será elevado.

Podemos calcularlo a mano para ver que no hay nada mágico en la función vif(). Tomemos height y regresémosla contra los otros dos predictores del modelo:

> aux <- lm(height ~ age + male, data=df_menores)
> 
> r2_aux <- summary(aux)$r.squared
> r2_aux
## [1] 0.894438
> 
> 1 / (1 - r2_aux)
## [1] 9.473102

Como el VIF multiplica la varianza del coeficiente, sobre el error estándar el efecto es \(\sqrt{VIF}\). Eso permite traducir cada valor a algo directamente interpretable: cuánto más ancho es el intervalo de confianza respecto del caso sin colinealidad.

\(R^2_{X_{j}\vert X_{-j}}\) VIF Multiplicador del error estándar
0.00 1.0 \(\times\) 1.00 (ausencia de colinealidad)
0.50 2.0 \(\times\) 1.41
0.80 5.0 \(\times\) 2.24
0.89 9.5 \(\times\) 3.08
0.90 10.0 \(\times\) 3.16

Es decir: un VIF de 10 no significa que el modelo esté “mal”, sino que el intervalo de confianza de ese coeficiente es unas 3 veces más ancho de lo que sería si la variable fuera independiente del resto.

Veamos ahora la función:

> 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. Noten que el valor de height coincide con el que calculamos a mano más arriba.

Estos números cierran el ejercicio que hicimos al principio. Pasemos los VIF a escala de error estándar:

> sqrt(vif(model1))
##   height      age     male 
## 3.077841 3.072036 1.021570

El error estándar de height debería ser aproximadamente 3 veces mayor en el modelo completo que en un modelo sin colinealidad. Eso es exactamente lo que observamos al comparar model1 con model2: 0.0247 contra 0.0083.

El VIF, entonces, no informa nada que no pudiéramos ver comparando modelos a mano. Su ventaja es que lo hace de una sola vez, para todos los predictores y sin necesidad de reestimar nada. En un modelo con quince variables, la comparación manual es inviable.

Antes de pasar a las soluciones, conviene fijar un criterio: la colinealidad no deteriora el ajuste ni la capacidad predictiva del modelo. Las predicciones dentro del rango de los datos siguen siendo tan buenas como antes. Lo único que se degrada es la posibilidad de atribuir un efecto a cada variable por separado. De ahí que la decisión dependa del objetivo:

  • Si el objetivo es predecir un VIF alto puede ignorarse.
  • Si el objetivo es interpretar coeficientes hay que resolverlo.

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 el ajuste de la regresión, ya que la presencia de colinealidad implica que la información que esa variable proporciona sobre la respuesta es redundante en presencia de las otras variables.

La segunda solución puede ser combinar las variables 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.

Una aclaración sobre el GVIF

Habrán notado que esta salida tiene tres columnas y la anterior una sola. La diferencia es que acá hay predictores categóricos de más de dos categorías. Cuando eso ocurre, la variable entra al modelo como varias dummies y no tiene un único coeficiente, así que car::vif() devuelve el GVIF (VIF generalizado), que evalúa el bloque completo de dummies:

  • Df: los grados de libertad, es decir la cantidad de dummies que genera la variable. nivel_ed_agg tiene 3 categorías, entonces Df = 2; class_eow_agg tiene 4, entonces Df = 3.
  • GVIF: crece con la cantidad de categorías, así que no es comparable con los umbrales de 5 y 10.
  • GVIF^(1/(2*Df)): es la versión corregida por Df, y es la que hay que mirar. Como está en escala de error estándar (no de varianza), se compara con \(\sqrt{5} = 2.24\) y \(\sqrt{10} = 3.16\); equivalentemente, se la puede elevar al cuadrado y comparar contra 5 y 10 directamente.

Hagamos esa última transformación para poder leer la salida con los umbrales de siempre:

> vif(lm_enes)[, 3]^2
##       t_hogar          v108          v109  nivel_ed_agg class_eow_agg 
##      1.089919      1.161298      1.063069      1.091319      1.028677

Todos los valores quedan muy por debajo de 5. Para las variables con Df = 1, GVIF coincide con el VIF común y las tres columnas dicen lo mismo.