7.4 Supuestos de MCO en regresión múltiple

En el modelo de regresión múltiple se amplian los tres supuestos de mínimos cuadrados del modelo de regresión simple (ver Capítulo 5) y se agrega un cuarto supuesto. Estos supuestos se presentan en el Concepto clave 6.4. No se entrará en detalles sobre los supuestos 1-3, ya que sus ideas se generalizan fácilmente al caso de regresores múltiples. La atención se centrará en el cuarto supuesto. Este supuesto descarta la correlación perfecta entre regresores.

Concepto clave 6.4

Los supuestos de mínimos cuadrados en el modelo de regresión múltiple

El modelo de regresión múltiple viene dado por

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_1 X_{2i} + \dots + \beta_k X_{ki} + u_i \ , \ i=1,\dots,n. \]

Los supuestos de MCO en el modelo de regresión múltiple son una extensión de los realizados para el modelo de regresión simple:

  1. Los regresores \((X_{1i}, X_{2i}, \dots, X_{ki}, Y_i) \ , \ i=1,\dots,n\), se constuye de manera que la suposición i.i.d. se mantiene.
  2. \(u_i\) es un término de error con media condicional cero dados los regresores; es decir,

\[ E(u_i\vert X_{1i}, X_{2i}, \dots, X_{ki}) = 0. \]

  1. Los valores atípicos grandes son poco probables, formalmente \(X_{1i},\dots,X_{ki}\) y \(Y_i\) tienen cuartos momentos finitos.
  2. No existe multicolinealidad perfecta.

Multicolinealidad

Multicolinealidad significa que dos o más regresores en un modelo de regresión múltiple están fuertemente correlacionados. Si la correlación entre dos o más regresores es perfecta; es decir, un regresor puede escribirse como una combinación lineal del otro(s), se tiene multicolinealidad perfecta. Si bien la multicolinealidad fuerte en general es desagradable, ya que hace que la varianza del estimador MCO sea grande (se discutirá esto a profundidad más adelante), la presencia de multicolinealidad perfecta hace que sea imposible resolver el estimador MCO; es decir, el modelo no puede estimarse en primer lugar.

La siguiente sección presenta algunos ejemplos de multicolinealidad perfecta y demuestra cómo lm() los trata.

Ejemplos de multicolinealidad perfecta

¿Cómo reacciona R si se intenta estimar un modelo con regresores perfectamente correlacionados?

La función lm produce una advertencia en la primera línea de la sección del coeficiente de la salida (1 no definido debido a singularidades) e ignorará el regresor que se supone que es una combinación lineal del otro(s). Considere el siguiente ejemplo donde se agrega otra variable FracEL, la fracción de estudiantes de inglés, a CASchools donde las observaciones son valores escalados de las observaciones para inglés y lo se usa como regresor junto con STR e inglés en un modelo de regresión múltiple. En este ejemplo, inglés y FracEL son perfectamente colineales. El código R es el siguiente.

# definir la fracción de estudiantes de inglés        
CASchools$FracEL <- CASchools$english / 100

# estimar el modelo
mult.mod <- lm(score ~ STR + english + FracEL, data = CASchools) 

# obtener un resumen del modelo
summary(mult.mod)                                                 
#> 
#> Call:
#> lm(formula = score ~ STR + english + FracEL, data = CASchools)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -48.845 -10.240  -0.308   9.815  43.461 
#> 
#> Coefficients: (1 not defined because of singularities)
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
#> STR          -1.10130    0.38028  -2.896  0.00398 ** 
#> english      -0.64978    0.03934 -16.516  < 2e-16 ***
#> FracEL             NA         NA      NA       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 14.46 on 417 degrees of freedom
#> Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
#> F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

La fila FracEL en la sección de coeficientes del resultado consta de entradas NA, ya que FracEL se excluyó del modelo.

Si se tuviera que calcular el MCO a mano, se encontraría con el mismo problema, ¡pero nadie estaría ayudando! El cálculo simplemente falla. ¿Por qué es esto? Tome el siguiente ejemplo:

Suponga que desea estimar un modelo de regresión lineal simple con una constante y un regresor único \(X\). Como se mencionó anteriormente, para que exista una multicolinealidad perfecta, \(X\) tiene que ser una combinación lineal de los otros regresores. Dado que el único otro regresor es una constante (piense en el lado derecho de la ecuación del modelo como \(\beta_0 \times 1 + \beta_1 X_i + u_i\) de modo que \(\beta_1\) siempre se multiplica por \(1\) para cada observación), \(X\) también tiene que ser constante. Por \(\hat\beta_1\) se tiene

\(\hat\beta_1 = \frac{\sum_{i = 1}^n (X_i - \bar{X})(Y_i - \bar{Y})} { \sum_{i=1}^n (X_i - \bar{X})^2} = \frac{\widehat{Cov}(X,Y)}{\widehat{Var}(X)}. \tag{6.7}\)

La varianza del regresor \(X\) está en el denominador. Dado que la varianza de una constante es cero, no se puede calcular esta fracción y \(\hat{\beta}_1\) no está definida.

Nota: En este caso especial, el denominador en (6.7) también es igual a cero. ¿Se puede mostrar eso?

Se pueden considerar dos ejemplos más en los que la selección de regresores induce una multicolinealidad perfecta. Primero, se debe suponer que se tiene la intención de analizar el efecto del tamaño de la clase en el puntaje de la prueba usando una variable ficticia que identifica las clases que no son pequeñas (\(NS\)). Se define que una escuela tiene el atributo \(NS\) cuando la proporción promedio de estudiantes por maestro de la escuela es de al menos \(12\),

\[ NS = \begin{cases} 0, \ \ \ \text{si STR < 12} \\ 1 \ \ \ \text{de lo contrario.} \end{cases} \]

Se agrega la columna correspondiente a CASchools y se estima un modelo de regresión múltiple con covariables computadora e inglés.

# si STR menor a 12, NS = 0, en caso contrario NS = 1
CASchools$NS <- ifelse(CASchools$STR < 12, 0, 1)

# estimar el modelo
mult.mod <- lm(score ~ computer + english + NS, data = CASchools)

# obtener un resumen del modelo
summary(mult.mod)                                                  
#> 
#> Call:
#> lm(formula = score ~ computer + english + NS, data = CASchools)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -49.492  -9.976  -0.778   8.761  43.798 
#> 
#> Coefficients: (1 not defined because of singularities)
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 663.704837   0.984259 674.319  < 2e-16 ***
#> computer      0.005374   0.001670   3.218  0.00139 ** 
#> english      -0.708947   0.040303 -17.591  < 2e-16 ***
#> NS                  NA         NA      NA       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 14.43 on 417 degrees of freedom
#> Multiple R-squared:  0.4291, Adjusted R-squared:  0.4263 
#> F-statistic: 156.7 on 2 and 417 DF,  p-value: < 2.2e-16

Nuevamente, el resultado de summary(mult.mod) indica que la inclusión de NS en la regresión haría inviable la estimación. ¿Que pasó aquí? Este es un ejemplo en el que se cometió un error lógico al definir el regresor NS: Al observar más de cerca \(NS\), la medida redefinida para el tamaño de la clase, revela que no existe una sola escuela con \(STR < 12\) por tanto, \(NS\) es igual a uno para todas las observaciones. Se puede verificar esto imprimiendo el contenido de CASchools\$NS o usando la función table(), vea ?Table.

table(CASchools$NS)
#> 
#>   1 
#> 420

Cascools$NS es un vector de \(420\), mientras que el conjunto de datos analizado incluye \(420\) observaciones. Esto obviamente viola la suposición 4 del Concepto clave 6.4: Las observaciones para el intercepto son siempre \(1\),

\[\begin{align*} intercept = \, & \lambda \cdot NS \end{align*}\]

\[\begin{align*} \begin{pmatrix} 1 \\ \vdots \\ 1 \end{pmatrix} = \, & \lambda \cdot \begin{pmatrix} 1 \\ \vdots \\ 1 \end{pmatrix} \\ \Leftrightarrow \, & \lambda = 1. \end{align*}\]

Dado que los regresores se pueden escribir como una combinación lineal entre sí, se enfrenta a multicolinenalidad perfecta y R excluye NS del modelo. Por lo tanto, el mensaje de advertencia es: ¡Piense cuidadosamente sobre cómo se relacionan los regresores en sus modelos!

Otro ejemplo de multicolinealidad perfecta se conoce como la trampa variable ficticia. Esto puede ocurrir cuando se utilizan múltiples variables ficticias como regresores. Un caso común para esto es cuando se utilizan variables ficticias para ordenar los datos en categorías mutuamente exclusivas. Por ejemplo, suponiendo que se tiene información espacial que indica si una escuela está ubicada en el norte, al oeste, al sur o al este de los EE. UU. Esto permite crear las variables ficticias

\[\begin{align*} North_i =& \begin{cases} 1 \ \ \text{Si se encuentra en el norte} \\ 0 \ \ \text{de lo contrario} \end{cases} \\ West_i =& \begin{cases} 1 \ \ \text{Si se encuentra en el oeste} \\ 0 \ \ \text{de lo contrario} \end{cases} \\ South_i =& \begin{cases} 1 \ \ \text{Si se encuentra en el sur} \\ 0 \ \ \text{de lo contrario} \end{cases} \\ East_i =& \begin{cases} 1 \ \ \text{Si se encuentra en el este} \\ 0 \ \ \text{de lo contrario}. \end{cases} \end{align*}\]

Dado que las regiones son mutuamente excluyentes, por cada escuela \(i=1,\dots,n\) se tiene

\[ North_I + West_i + South_I + East_I = 1. \]

Se enfrenta con problemas al tratar de estimar un modelo que incluye una constante y las cuatro variables ficticias en el modelo; por ejemplo,

\[ TestScore = \beta_0 + \beta_1 \times STR + \beta_2 \times english + \beta_3 \times North_i + \beta_4 \times West_i + \beta_5 \times South_i + \beta_6 \times East_i + u_i \tag{6.8}\]

desde entonces, para todas las observaciones \(i=1,\dots,n\) el término constante es una combinación lineal de las variables ficticias:

\[\begin{align} intercept = \, & \lambda_1 \cdot (North + West + South + East) \\ \begin{pmatrix} 1 \\ \vdots \\ 1\end{pmatrix} = \, & \lambda_1 \cdot \begin{pmatrix} 1 \\ \vdots \\ 1\end{pmatrix} \\ \Leftrightarrow \, & \lambda_1 = 1 \end{align}\]

y se tiene una perfecta multicolinealidad. Por lo tanto, la “trampa de la variable ficticia” significa no prestar atención e incluir falsamente dummies exhaustivas y una constante en un modelo de regresión.

¿Cómo la función lm() utiliza una regresión como (6.8)? Primero se deben generar algunos datos categóricos artificiales y agregar una nueva columna llamada directions a CASchools y ver cómo lm() se comporta cuando se le pide que estime el modelo.

# sembrar la semilla para la reproducibilidad
set.seed(1)

# generar datos artificiales en la ubicación
CASchools$direction <- sample(c("West", "North", "South", "East"), 
                              420, 
                              replace = T)

# estimar el modelo
mult.mod <- lm(score ~ STR + english + direction, data = CASchools)

# obtener un resumen del modelo
summary(mult.mod)                                                 
#> 
#> Call:
#> lm(formula = score ~ STR + english + direction, data = CASchools)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -49.603 -10.175  -0.484   9.524  42.830 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    684.80477    7.54130  90.807  < 2e-16 ***
#> STR             -1.08873    0.38153  -2.854  0.00454 ** 
#> english         -0.65597    0.04018 -16.325  < 2e-16 ***
#> directionNorth   1.66314    2.05870   0.808  0.41964    
#> directionSouth   0.71619    2.06321   0.347  0.72867    
#> directionWest    1.79351    1.98174   0.905  0.36598    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 14.5 on 414 degrees of freedom
#> Multiple R-squared:  0.4279, Adjusted R-squared:  0.421 
#> F-statistic: 61.92 on 5 and 414 DF,  p-value: < 2.2e-16

Observe que R resuelve el problema por sí solo al generar e incluir las variables ficticias directionNorth, directionSouth y directionWest pero omitiendo directionEast. Por supuesto, la omisión de todas las demás ficticias lograría lo mismo. Otra solución sería excluir la constante e incluir todas las variables ficticias.

¿Significa esto que se pierde la información sobre las escuelas ubicadas en el Este? Afortunadamente, este no es el caso: la exclusión de directEast simplemente altera la interpretación de las estimaciones de coeficientes en las variables ficticias restantes de absoluto a relativo. Por ejemplo, la estimación del coeficiente en directionNorth establece que, en promedio, los puntajes de las pruebas en el norte son aproximadamente \(1.61\) puntos más altos que en el este.

Un último ejemplo considera el caso en el que una relación lineal perfecta surge de regresores redundantes. Suponiendo que se tiene un regresor \(PctES\), el porcentaje de hablantes de inglés en la escuela donde

\[ PctES = 100 - PctEL\]

y tanto \(PctES\) como \(PctEL\) se incluyen en un modelo de regresión. Un regresor es redundante ya que el otro transmite la misma información. Dado que obviamente este es un caso en el que los regresores se pueden escribir como una combinación lineal, se termina con una multicolinealidad perfecta, nuevamente.

Haciendo esto en R.

# porcentaje de hablantes de inglés
CASchools$PctES <- 100 - CASchools$english

# estimar el modelo
mult.mod <- lm(score ~ STR + english + PctES, data = CASchools)

# obtener un resumen del modelo
summary(mult.mod)                                                 
#> 
#> Call:
#> lm(formula = score ~ STR + english + PctES, data = CASchools)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -48.845 -10.240  -0.308   9.815  43.461 
#> 
#> Coefficients: (1 not defined because of singularities)
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
#> STR          -1.10130    0.38028  -2.896  0.00398 ** 
#> english      -0.64978    0.03934 -16.516  < 2e-16 ***
#> PctES              NA         NA      NA       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 14.46 on 417 degrees of freedom
#> Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
#> F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

Una vez más, lm() se niega a estimar el modelo completo usando MCO y excluye PctES.

Lo anterior se debe a la multicolinealidad perfecta, sus consecuencias para el estimador MCO en modelos generales de regresión múltiple se pueden demostrar utilizando notación matricial.

Multicolinealidad imperfecta

A diferencia de la multicolinealidad perfecta, la multicolinealidad imperfecta es, hasta cierto punto, un problema menor. De hecho, la multicolinealidad imperfecta es la razón por la que se está interesado en estimar modelos de regresión múltiple en primer lugar: El estimador MCO permite aislar las influencias de los regresores correlacionados en la variable dependiente. Si no fuera por estas dependencias, no habría razón para recurrir a un enfoque de regresión múltiple y simplemente se podría trabajar con un modelo de regresor único. Sin embargo, este rara vez es el caso en las aplicaciones. Ya se sabe que ignorar las dependencias entre regresores que influyen en la variable de resultado tiene un efecto adverso sobre los resultados de la estimación.

Entonces, ¿cuándo y por qué es un problema la multicolinealidad imperfecta? Suponga que tiene el modelo de regresión

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + u_i \tag{6.9} \]

y está interesado en estimar \(\beta_1\), el efecto en \(Y_i\) de un cambio de una unidad en \(X_{1i}\), mientras se mantiene constante \(X_{2i}\). No sabe que el verdadero modelo incluye \(X_2\). Sigue un razonamiento y agrega \(X_2\) como una covariable al modelo para abordar un posible sesgo de variable omitida. Está seguro de que \(E(u_i\vert X_{1i}, X_{2i})=0\) y que no existe razón alguna para sospechar una violación de los supuestos 2 y 3 del Concepto clave 6.4. Si \(X_1\) y \(X_2\) están altamente correlacionadas, el MCO tendrá dificultades para estimar con precisión \(\beta_1\). Eso significa que aunque \(\hat\beta_1\) es un estimador consistente e imparcial para \(\beta_1\), tiene una gran variación debido a que \(X_2\) está incluido en el modelo. Si los errores son homocedásticos, este problema se puede comprender mejor con la fórmula de la varianza de \(\hat\beta_1\) en el modelo (6.9):

\[ \sigma^2_{\hat\beta_1} = \frac{1}{n} \left( \frac{1}{1-\rho^2_{X_1,X_2}} \right) \frac{\sigma^2_u}{\sigma^2_{X_1}}. \tag{6.10} \]

Primero, si \(\rho_{X_1,X_2}=0\), es decir, si no existe correlación entre ambos regresores, incluir \(X_2\) en el modelo no influye en la varianza de \(\hat\beta_1\).

En segundo lugar, si \(X_1\) y \(X_2\) están correlacionados, \(\sigma^2_{\hat\beta_1}\) es inversamente proporcional a \(1-\rho^2_{X_1,X_2}\) así que cuanto más fuerte sea la correlación entre \(X_1\) y \(X_2\), el menor es \(1-\rho^2_{X_1,X_2}\) y, por lo tanto, mayor es la varianza de \(\hat\beta_1\).

En tercer lugar, aumentar el tamaño de la muestra ayuda a reducir la varianza de \(\hat\beta_1\). Por supuesto, esto no se limita al caso de dos regresores: En regresiones múltiples, la multicolinealidad imperfecta infla la varianza de uno o más estimadores de coeficientes. Es una cuestión empírica qué estimaciones de coeficientes se ven gravemente afectadas por esto y cuáles no. Cuando el tamaño de la muestra es pequeño, a menudo uno se enfrenta a la decisión de aceptar la consecuencia de agregar un gran número de covariables (mayor varianza) o utilizar un modelo con pocos regresores (posible sesgo de variable omitida). Esto se llama compensación de sesgo-varianza.

En resumen, las consecuencias indeseables de la multicolinealidad imperfecta generalmente no son el resultado de un error lógico cometido por el investigador (como suele ser el caso de la multicolinealidad perfecta) sino más bien un problema que está vinculado a los datos utilizados, el modelo que se va a estimar y la pregunta de investigación en cuestión.

Estudio de simulación: Multicolinealidad imperfecta

Resulta importante realizar un estudio de simulación para ilustrar los problemas esbozados anteriormente.

  1. Usando (6.9) como proceso de generación de datos y al elegir \(\beta_0 = 5\), \(\beta_1 = 2.5\) y \(\beta_2 = 3\), así como \(u_i\) (que es un término de error distribuido como \(\mathcal{N}(0,5)\)). En un primer paso, se toma una muestra de los datos del regresor de una distribución normal bivariada:

\[ X_i = (X_{1i}, X_{2i}) \overset{i.i.d.}{\sim} \mathcal{N} \left[\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 10 & 2.5 \\ 2.5 & 10 \end{pmatrix} \right] \]

Es sencillo ver que la correlación entre \(X_1\) y \(X_2\) en la población es bastante baja:

\[ \rho_{X_1,X_2} = \frac{Cov(X_1,X_2)}{\sqrt{Var(X_1)}\sqrt{Var{(X_2)}}} = \frac{2.5}{10} = 0.25 \]

  1. A continuación, se estima el modelo (6.9) y se guardan las estimaciones para \(\beta_1\) y \(\beta_2\). Esto se repite \(10000\) veces con un bucle for, por lo que termina con una gran cantidad de estimaciones que permiten describir las distribuciones de \(\hat\beta_1\) y \(\hat\beta_2\).

  2. Se repiten los pasos 1 y 2, pero aumentando la covarianza entre \(X_1\) y \(X_2\) de \(2.5\) a \(8.5\) de manera que la correlación entre los regresores sea alta:

\[ \rho_{X_1,X_2} = \frac{Cov(X_1,X_2)}{\sqrt{Var(X_1)}\sqrt{Var{(X_2)}}} = \frac{8.5}{10} = 0.85 \]

  1. Para evaluar el efecto sobre la precisión de los estimadores de aumentar la colinealidad entre \(X_1\) y \(X_2\), se estiiman y comparan las varianzas de \(\hat\beta_1\) y \(\hat\beta_2\) .
# cargar paquetes
library(MASS)
library(mvtnorm)

# establecer el número de observaciones
n <- 50

# inicializar vectores de coeficientes
coefs1 <- cbind("hat_beta_1" = numeric(10000), "hat_beta_2" = numeric(10000))
coefs2 <- coefs1

# sembrar semilla
set.seed(1)

# muestreo y estimación de bucles
for (i in 1:10000) {
  
  # para cov(X_1, X_2) = 0.25
  X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
  u <- rnorm(n, sd = 5)
  Y <- 5 + 2.5 * X[, 1] + 3 * X[, 2] + u
  coefs1[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
  
  # para cov(X_1, X_2) = 0.85
  X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
  Y <- 5 + 2.5 * X[, 1] + 3 * X[, 2] + u
  coefs2[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
  
}

# obtener estimaciones de varianza
diag(var(coefs1))
#> hat_beta_1 hat_beta_2 
#> 0.05674375 0.05712459
diag(var(coefs2))
#> hat_beta_1 hat_beta_2 
#>  0.1904949  0.1909056

Se está interesado en las variaciones que son los elementos diagonales. Se puede ver que debido a la alta colinealidad, las varianzas de \(\hat\beta_1\) y \(\hat\beta_2\) se han más que triplicado, lo que significa que es más difícil estimar con precisión los coeficientes verdaderos.