Multicolinealidad

En el análisis de regresión múltiple, la naturaleza de las relaciones entre las variables predictoras y la variable de respuesta suelen ser de especial interés. Algunas preguntas frecuentes son:

¿Cuál es la importancia relativa de los efectos de las distintas variables predictoras?
¿Cuál es la magnitud del efecto de una determinada variable predictora sobre la variable de respuesta?
¿Puede eliminarse alguna variable predictora del modelo porque tiene poco o ningún efecto sobre la variable de respuesta?
¿Debería considerarse la posible inclusión de variables predictoras aún no incluidas en el modelo?

Si las variables predictoras incluidas en el modelo:

No están correlacionadas entre sí y,
No están correlacionadas con otras variables predictoras relacionadas con la variable de respuesta pero omitidas en el modelo.

Se pueden dar respuestas relativamente sencillas a estas preguntas. Desgraciadamente, en muchas situaciones no experimentales de los negocios, la economía y las ciencias sociales y biológicas, las variables predictoras o explicativas tienden a estar correlacionadas entre sí y con otras variables que están relacionadas con la variable de respuesta pero que no se incluyen en el modelo. Por ejemplo, en una regresión que pretende analizar el gasto alimentario familiar teniendo en cuenta la renta familiar, el ahorro familiar y la edad del cabeza de familia, las variables explicativas estarán correlacionadas entre sí. Además, también estarán correlacionadas con otras variables socioeconómicas no incluidas en el modelo que sí afectan al gasto alimentario familiar, como el tamaño de la familia.

Cuando las variables predictoras están correlacionadas entre sí, se dice que existe intercorrelación o multicolinealidad entre ellas. Se pueden observar algunos problemas con la multicolinealidad, como sigue.

Variables predictoras no correlacionadas

La Tabla 1 contiene datos que pretenden analizar el efecto del tamaño del equipo de trabajo y las bonificaciones en la productividad del equipo.

Tabla 1: Datos modelo regresoras no correlacionadas

Observación \(i\)	\(x_{i1}\): Tamaño del equipo	\(x_{i2}\): Bonificación	\(y_i\): productividad del equipo
1	4	2	42
2	4	2	39
3	4	3	48
4	4	3	51
5	6	2	49
6	6	2	53
7	6	3	61
8	6	3	60

Las variables predictoras \(x_1\) y \(x_2\) no están correlacionadas, \(\rho(x_1,x_2)\) = 0, como se muestra a continuación:

#Datos
tamano <- c(4,4,4,4,6,6,6,6)
bonificacion <- c(2,2,3,3,2,2,3,3)
productividad <- c(42,39,48,51,49,53,61,60)

# Correlación
correlacion<-cor(tamano, bonificacion)

Si hacemos el ajuste de \(\beta_k\) para los siguientes tres modelos:

Variable respuesta: \(y\): productividad, variables regresoras: \(x_1\): tamaño del equipo, \(x_2\): bonificación.

#Datos
tamano <- c(4,4,4,4,6,6,6,6)
bonificacion <- c(2,2,3,3,2,2,3,3)
productividad <- c(42,39,48,51,49,53,61,60)

# Modelo
modelo <- lm(productividad~tamano+bonificacion)
resumen<-summary(modelo)
betas <- resumen$coefficients[,1]
betas

(Intercept) tamano bonificacion 0.375 5.375 9.250

Obteniendo la siguiente función de regresión:

\[E[y_i] = 0.375 + 5.375x_{i1}+9.25x_{i2}\]

Variable respuesta: \(y\): productividad, variable regresora: \(x_1\): tamaño del equipo.

#Datos
tamano <- c(4,4,4,4,6,6,6,6)
bonificacion <- c(2,2,3,3,2,2,3,3)
productividad <- c(42,39,48,51,49,53,61,60)

# Modelo
modelo <- lm(productividad~tamano)
resumen<-summary(modelo)
betas <- resumen$coefficients[,1]
betas

(Intercept) tamano 23.500 5.375

Obteniendo la siguiente función de regresión:

\[E[y_i] =23.5 + 5.375x_{i1}\]

Variable respuesta: \(y\): productividad, variable regresora:\(x_2\): bonificación.

#Datos
tamano <- c(4,4,4,4,6,6,6,6)
bonificacion <- c(2,2,3,3,2,2,3,3)
productividad <- c(42,39,48,51,49,53,61,60)

# Modelo
modelo <- lm(productividad~bonificacion)
resumen<-summary(modelo)
betas <- resumen$coefficients[,1]
betas

(Intercept) bonificacion 27.25 9.25

Obteniendo la siguiente función de regresión:

\[E[y_i] =27.25 + 9.25x_{i2}\]

Un resultado importante a tener en cuenta de los resultados anteriores es que el coeficiente de regresión \(\hat{\beta_1}\) es el mismo, \(\hat{\beta_1}=5.375\) es el mismo, aunque se incluya o no, la variable regresora \(x_2\). Los mismo pasa con \(\hat{\beta_2}=9.25\). Lo anterior es el resultado de tener variables regresoras no correlacionadas. Así, cuando las variables predictoras no están correlacionadas, los efectos que les atribuye un modelo de regresión de primer orden son los mismos independientemente de cuál de ellas se incluya en el modelo.

Variables predictoras correlacionadas

Para empezar a observar la naturaleza del problema de multicolinealidad, se empleará un ejemplo de dos variables predictoras perfectamente correlacionadas. Los datos se muestran en la Tabla 2:

Tabla 2: Datos variables regresoras correlacionadas

\(i\)	\(x_{i1}\)	\(x_{i2}\)	\(y_i\)
1	2	6	23
2	8	9	83
3	6	8	63
4	10	10	103

Las variables predictoras \(x_1\) y \(x_2\) están correlacionadas, \(\rho(x_1,x_2)\) = 1, como se muestra a continuación:

#Datos
y <- c(23,83,63,103)
x1 <- c(2,8,6,10)
x2 <- c(6,9,8,10)

# Correlación
correlacion<-cor(x1, x2)
correlacion

[1] 1

Si se desea ajustar un modelo mediante mínimos cuadrados en R obtenemos lo siguiente:

#Datos
y <- c(23,83,63,103)
x1 <- c(2,8,6,10)
x2 <- c(6,9,8,10)

#Modelo
modelo <- lm(y~x1+x2)
summary(modelo)

Call: lm(formula = y ~ x1 + x2)

Residuals: 1 2 3 4 1.470e-15 -5.515e-15 -1.834e-16 4.228e-15

Coefficients: (1 not defined because of singularities) Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.000e+00 6.065e-15 4.946e+14 <2e-16 x1 1.000e+01 8.493e-16 1.177e+16 <2e-16 x2 NA NA NA NA

Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 5.024e-15 on 2 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 1.386e+32 on 1 and 2 DF, p-value: < 2.2e-16

La frase Coefficients: (1 not defined because of singularities) aparece debido a que uno de los coeficientes del modelo de regresión no pudo ser estimado debido a una colinealidad.

Sin embargo si tomamos los dos siguientes modelos de ejemplo:

\[\begin{align} \hat{y} = -87 + x_1 + 18x_2\\ \\ \hat{y} = -7 + 9x_1 + 2x_2 \end{align}\]

Observamos su comportamiento en la Figura 1:

Figura 1: Modelos ajustados

Si se calculan los valores ajustados para ambos modelos se obtienen los valores mostrados en Tabla 3

#Datos
y <- c(23,83,63,103)
x1 <- c(2,8,6,10)
x2 <- c(6,9,8,10)

#Modelo
modelo1 <- -87+x1+18*x2
modelo2 <- -7+9*x1+2*x2

modelo1

[1] 23 83 63 103

modelo2

[1] 23 83 63 103

Tabla 3: Valores ajustados para dos modelos

\(x_{i1}\)	\(x_{i2}\)	\(y_i\)	\(\hat{y_i}\):modelo 1	\(\hat{y_i}\): modelo 2
2	6	23	23	23
8	9	83	83	83
6	8	63	63	63
10	10	103	103	103

Cómo se observa en la Tabla 3 los modelos ajustan perfectamente a los datos (\(y_i = \hat{y_i}:modelo~1=\hat{y_i}:modelo~2\))

De hecho, se puede demostrar que infinitas funciones de respuesta se ajustarán perfectamente a los datos de este problema. La razón es que las variables predictoras \(x_1\) y \(x_2\) están perfectamente relacionadas. según la relación: \(x_2 =5 +0.5x_1\). Lo ocurrido en este ejemplo tiene dos implicaciones interesantes:

La perfecta relación entre \(x_1\), y \(x_2\) no inhibió nuestra capacidad para obtener un buen ajuste a los datos.
Dado que muchas funciones de respuesta diferentes proporcionan el mismo ajuste, no podemos interpretar ningún conjunto de coeficientes de regresión como reflejo de los efectos de las diferentes variables predictoras.

Efectos de la multicolinealidad

El hecho de que algunas o todas las variables predictoras estén correlacionadas entre sí no inhibe, en general, nuestra capacidad para obtener un buen ajuste ni tiende a afectar a las inferencias sobre las respuestas medias o las predicciones de nuevas observaciones, siempre que estas inferencias se realicen dentro de la región de observaciones.
La contrapartida en la vida real de que muchas funciones de regresión diferentes se ajusten igual de bien a los datos de nuestro ejemplo idealizado es que los coeficientes de regresión estimados tienden a tener una gran variabilidad muestral cuando las variables predictoras están muy correlacionadas. Así, los coeficientes de regresión estimados tienden a variar mucho de una muestra a otra cuando las variables predictoras están muy correlacionadas. Como resultado, es posible que sólo se disponga de información imprecisa sobre los verdaderos coeficientes de regresión individuales. De hecho, muchos de los coeficientes de regresión estimados pueden no ser estadísticamente significativos aunque exista una relación estadística definida entre la variable de respuesta y el conjunto de variables predictoras.
La interpretación común de un coeficiente de regresión como medida del cambio en el valor esperado de la variable de respuesta cuando la variable de predicción dada se incrementa en una unidad mientras que todas las demás variables de predicción se mantienen constantes no es totalmente aplicable cuando existe multicolinealidad. Puede ser conceptualmente factible pensar en variar una variable predictora y mantener constantes las demás, pero en la práctica puede no ser posible hacerlo en el caso de variables predictoras que estén muy correlacionadas. Por ejemplo, en un modelo de regresión para predecir el rendimiento de un cultivo a partir de la cantidad de precipitaciones y las horas de sol, la relación entre las dos variables predictoras hace que no sea realista considerar la posibilidad de variar una manteniendo constante la otra. Por lo tanto, la simple interpretación de los coeficientes de regresión como medida de los efectos marginales no suele estar justificada cuando las variables predictoras están muy correlacionadas.

Diagnóstico de la multicolinealidad - Factor de Inflación de la Varianza (VIF)

Un método formal para detectar la presencia de multicolinealidad es el uso de Factores de Inflación de la Varianza (\(VIF\)). Estos factores miden cuánto aumentan las varianzas de los coeficientes de regresión estimados en comparación con cuando las variables predictoras no están correlacionadas linealmente.

Para empezar a entender cómo se calculan los \(VIF\). Se puede recordar que la matriz de varianza-covarianza de los coeficientes de regresión estimados se calcula como se muestra en la Ecuación 1:

\[\begin{align} \sigma^2(\hat{\beta}) = \sigma^2(X^{'}X)^{-1} \end{align} \tag{1}\]

Para propósitos de medir el impacto de la multicolinealidad, es útil trabajar con el modelo de regresión estandarizado, que se obtiene transformando los datos como se muestra en la Ecuación 2:

\[\begin{align} Y_i^* = \frac{1}{\sqrt{n-1}} \left( \frac{Y_i - \bar{Y}}{S_Y} \right)\\ \\ X_{ik}^* = \frac{1}{\sqrt{n-1}} \left( \frac{X_{ik} - \bar{X_k}}{S_k} \right)\\ \\ k=1,2,...,p-1\\ \\ S_y = \sqrt{\frac{\sum_i(Y_i-\bar{Y})^2}{n-1}}\\ \\ S_k = \sqrt{\frac{\sum_i(X_{ik}-\bar{X_{k}})^2}{n-1}} \end{align} \tag{2}\]

Obteniendo el modelo estandarizado que se muestra en la Ecuación 3:

\[\begin{align} Y_i^* = \beta_1^*X_{i1}^* + \beta_2^*X_{i2}^*+...+\beta_{p-1}^*X_{i~~p-1}^*+\epsilon_i^* \end{align} \tag{3}\]

Los parámetros del modelo estandarizado se relacionan con los parámetros del modelo sin estandarizar como se muestra en la Ecuación 4

\[\begin{align} \beta_k = \left( \frac{Sy}{Sk} \right)\beta_k^*\\ \end{align} \tag{4}\]

Ejemplo de modelo estandarizado

Dwaine Studios, Inc. gestiona estudios de fotografía en \(21\) ciudades intermedias. Estos estudios están especializados en retratos de niños. La empresa está considerando expandirse a otras ciudades intermedias y desea investigar si las ventas en una ciudad pueden predecirse a partir del número de personas de \(16\) años o menos en la comunidad y el ingreso per cápita en la comunidad. En la Tabla 4 se muestran los datos de estas variables correspondientes al año más reciente de las \(21\) ciudades en las que opera actualmente Dwaine Studios.

Tabla 4: Datos ejemplo ventas de Dwaine Studios.

Observación \(i\)	\(x_1\): Población menor 16 años (miles de personas)	\(x_2\): Ingreso per cápita (miles de dólares)	\(y\): ventas (miles de dólares)
1	68.5	16.7	174.4
2	45.2	16.8	164.4
3	91.3	18.2	244.2
4	47.8	16.3	154.6
5	46.9	17.3	181.6
6	66.1	18.2	207.5
7	49.5	15.9	152.8
8	52.0	17.2	163.2
9	48.9	16.6	145.4
10	38.4	16.0	137.2
11	87.9	18.3	241.9
12	72.8	17.1	191.1
13	88.4	17.4	232.0
14	42.9	15.8	145.3
15	52.5	17.8	161.1
16	85.7	18.4	209.7
17	41.3	16.5	146.4
18	51.7	16.3	144.0
19	89.6	18.1	232.6
20	82.7	19.1	224.1
21	52.3	16.0	166.5

**Transformación usando Ecuación 2 de la variable respuesta, Tabla 5*

Tabla 5: Transformación de variable respuesta.

Observación \(i\)	\(Y_i\): ventas	\(Y_i^*\): ventas transformada
1	174.4	-0.0463679
2	164.4	-0.1081526
3	244.2	0.3848891
4	154.6	-0.1687016
5	181.6	-0.0018830
6	207.5	0.1581393
7	152.8	-0.1798228
8	163.2	-0.1155668
9	145.4	-0.2255435
10	137.2	-0.2762069
11	241.9	0.3706786
12	191.1	0.0568125
13	232.0	0.3095118
14	145.3	-0.2261613
15	161.1	-0.1285415
16	209.7	0.1717320
17	146.4	-0.2193650
18	144.0	-0.2341933
19	232.6	0.3132189
20	224.1	0.2607019
21	166.5	-0.0951778

Transformación usando Ecuación 2 de las variables regresoras, Tabla 6

Tabla 6: Transformación de variable respuesta.

Observación \(i\)	\(X_1\): población	\(X_1^*\)	\(X_2\): ingreso	\(X_2^*\)
1	68.5	0.0778281	16.7	-0.1020521
2	45.2	-0.2019757	16.8	-0.0790081
3	91.3	0.3516275	18.2	0.2436083
4	47.8	-0.1707529	16.3	-0.1942282
5	46.9	-0.1815608	17.3	0.0362120
6	66.1	0.0490071	18.2	0.2436083
7	49.5	-0.1503381	15.9	-0.2864043
8	52.0	-0.1203162	17.2	0.0131680
9	48.9	-0.1575433	16.6	-0.1250961
10	38.4	-0.2836352	16.0	-0.2633603
11	87.9	0.3107978	18.3	0.2666523
12	72.8	0.1294657	17.1	-0.0098760
13	88.4	0.3168022	17.4	0.0592561
14	42.9	-0.2295958	15.8	-0.3094484
15	52.5	-0.1143118	17.8	0.1514322
16	85.7	0.2843786	18.4	0.2896963
17	41.3	-0.2488098	16.5	-0.1481402
18	51.7	-0.1239188	16.3	-0.1942282
19	89.6	0.3312127	18.1	0.2205643
20	82.7	0.2483523	19.1	0.4510045
21	52.3	-0.1167136	16.0	-0.2633603

Ajuste del modelo con variables transformadas

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Transformación
sk_poblacion <- sqrt(sum((poblacion-mean(poblacion))^2)/(length(poblacion)-1))
sk_ingreso <- sqrt(sum((ingreso-mean(ingreso))^2)/(length(ingreso)-1))
poblacion_barra <- mean(poblacion)
ingreso_barra <- mean(ingreso)
poblacion_transf <- (1/sqrt(length(poblacion)-1))*((poblacion-poblacion_barra)/sk_poblacion)
ingreso_transf <- (1/sqrt(length(ingreso)-1))*((ingreso-ingreso_barra)/sk_ingreso)
sy <- sqrt(sum((ventas-mean(ventas))^2)/(length(ventas)-1))
y_barra <- mean(ventas)
y_transf <- (1/sqrt(length(ventas)-1))*((ventas-y_barra)/sy)

#Estimacion de parámetros
X <- matrix(c(rep(1, 21), poblacion_transf,ingreso_transf), nrow=21, ncol=3)
Y <- matrix(y_transf, ncol=1)
BETAS_TRANS <- solve(t(X)%*%X)%*%(t(X)%*%Y)
BETAS_TRANS

          [,1]

[1,] -1.225996e-17 [2,] 7.483670e-01 [3,] 2.511039e-01

Por lo que el modelo estimado sería:

\[Y_i^* = 0.7483670X_1^* + 0.2511039X_2^*\]

La matriz de varianza-covarianza de los coeficientes de regresión estimados estandarizados que establece que la matriz \(X^{'}X\) se denota \(r_{XX}^{-1}\), por lo que la matriz de varianza y covarianza quedaría como se muestra en la Ecuación 5:

\[\begin{align} \sigma^2(\hat{\beta^*}) = (\sigma^{*})^2r_{XX}^{-1} \end{align} \tag{5}\]

Donde \(r_{XX}\) es la matriz de coeficientes de correlaciones entre las variables \(X\). Definida como sigue, Ecuación 6:

\[\begin{align} r_{XX}= \begin{bmatrix} 1 & r_{12}&...& r_{1,p-1}\\ r_{21} & 1 &...& r_{2, p-1}\\ . & . &...&. \\ . & . &...&.\\ . & . &...&.\\ r_{p-1,1} & r_{p-1,2}&...&1 \end{bmatrix} \end{align} \tag{6}\]

\((\sigma^*)^2\) corresponde a la varianza del término del error para el modelo transformado.

De la Ecuación 5 se puede enotar que la varianza de \(\hat{\beta_k}~~;~~ k=1,2,...,p-1\) se puede expresar como se muestra en la Ecuación 7:

\[\begin{align} \sigma^2(\hat{\beta^*}) = (\sigma^{*})^2(VIF)_k \end{align} \tag{7}\]

\(VIF_k\) denota el \(k-ésimo\) elemento de la diagonal de la matriz \(r_{XX}^{-1}\)

El elemento de la diagonal \(VIF_k\) de la matriz \(r_{XX}^{-1}\) es llamado Factor de Inflación de la Varianza para \(\hat{\beta_k}\). El factor de inflación de la varianza \((VIF)_k\) se analiza de la siguiente manera:

\((VIF)_k=1\): la variable \(X_k\) no está linearmente relacionada con las otras variables \(X\).
\((VIF)_k > 1\): indica varianza inflada para \(\hat{\beta_k}\) como resultado de la multicolinealidad entre las variables \(X\).
\((VIF)_k \rightarrow \infty\): indica que \(X_k\) tiene correlación lineal perfecta con otro variable \(X\) dentro del modelo.

Una regla de decición para el \((VIF)_k\) podría ser la siguiente:

Si \((VIF)_k \leq 5\) no hay problemas de multicolinealidad.
Si \(5<(VIF)_k \leq 10\) hay problemas de multicolinealidad moderados.
\((VIF)_k > 10\) hay problemas de multicolinealidad graves.

Ejemplo de cálculo de \((VIF)_k\)

VIF mediante matriz de correlación \(r_{XX}\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Matriz de correlación
datos <- data.frame(poblacion, ingreso)
matriz_correlacion <- cor(datos)
matriz_correlacion

      poblacion   ingreso

poblacion 1.0000000 0.7812993 ingreso 0.7812993 1.0000000

#Matriz de correlación inversa
rxx_inv <- solve(matriz_correlacion)
rxx_inv

      poblacion   ingreso

poblacion 2.566924 -2.005536 ingreso -2.005536 2.566924

Se obtiene que \((VIF)_{población}=2.566924\), \((VIF)_{ingreso}=2.566924\).

VIF usando sentencias de R

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#VIF
library(car)
vif_k <- vif(modelo)
vif_k

poblacion ingreso 2.566924 2.566924