Regresión Dummy

Author

Rubén Cool

Introducción

En la regresión lineal múltiple clásica, tus predictores \(x_1, x_2, \dots, x_k\) son continuos o cuantitativos (como edad, horas de estudio, temperatura o ingresos). Pero en la vida real, muchos factores clave son cualitativos o categóricos:

  • Turno escolar: Matutino vs. Vespertino.

  • Nivel educativo: Licenciatura, Maestría, Doctorado.

  • Tipo de institución: Pública vs. Privada.

  • Tratamiento: Grupo de Control vs. Grupo de Tratamiento.

No podemos simplemente poner en la ecuación \(x = \text{Pública}\), porque la matemática exige números. Tampoco es correcto asignar valores arbitrarios como \(\text{Pública} = 1\), \(\text{Privada} = 2\), \(\text{Mixta} = 3\) en una sola variable, porque eso asumiría implícitamente que la distancia entre Privada y Pública es la misma que entre Mixta y Privada, e impondría un orden numérico artificial.

Ocasionalmente se debe introducir variables categóricas (o cualitativas) con dos o más categorias. Por ejemplo, género (hombre, mujer), estado civil, partido politico. Estas se pueden representar como variables dummy o indicadoras. Estas variables toman dos (DOS) valores usualmente, cero (0) y uno (1). Los dos valores significan que la observación pertenece a una de dos categorias. Las variables dummy o indicadoras sirven para identificar categorías o clase a las que pertenecen las observaciones.

En la práctica, la variable respuesta \(y\) suele depender no solo de covariables continuas \(x\), sino también de atributos cualitativos \(D\) con \(r\) niveles (ej. tipo de riesgo, región geográfica o especie).

Ejemplo 1

Consideremos una muestra reducida de \(6\) pólizas de seguro donde analizamos el Costo del Siniestro (\(y\)), la Edad del Vehículo (\(x\)) y la Zona de Cobertura (Zona A y Zona B):

library(knitr)
library(car)
Warning: package 'car' was built under R version 4.6.1
Cargando paquete requerido: carData
Warning: package 'carData' was built under R version 4.6.1
library(stats)
# Creación del conjunto de datos
datos <- data.frame(
  Poliza = 1:6,
  y = c(12, 15, 18, 22, 26, 30),
  x = c(2, 3, 5, 8, 9, 11),
  Zona = c("A", "A", "A", "B", "B", "B"),
  D_A = c(1, 1, 1, 0, 0, 0),
  D_B = c(0, 0, 0, 1, 1, 1))

kable(datos, align = "c", caption = "Tabla de Datos Sintéticos")
Tabla de Datos Sintéticos
Poliza y x Zona D_A D_B
1 12 2 A 1 0
2 15 3 A 1 0
3 18 5 A 1 0
4 22 8 B 0 1
5 26 9 B 0 1
6 30 11 B 0 1

Escenario 1: Multicolinealidad Perfecta (Trampa de la Variable Dummy)

Intentamos estimar el modelo especificando la constante (\(\mathbf{1}_n\)), la covariable continua \(x\) y ambas variables indicadoras \(D_A\) y \(D_B\):

\(y_i = \beta_0 + \beta_1 x_i + \gamma_A D_{i,A} + \gamma_B D_{i,B} + \varepsilon_i\)

Matriz de Diseño \(\mathbf{X}\) y Matriz de Momentos \(\mathbf{X}^t\mathbf{X}\)

# Construcción de la matriz de diseño con ambas dummies
X_perf <- cbind(
  Intercepto = 1,
  X = datos$x,
  D_A = datos$D_A,
  D_B = datos$D_B)
X_perf
     Intercepto  X D_A D_B
[1,]          1  2   1   0
[2,]          1  3   1   0
[3,]          1  5   1   0
[4,]          1  8   0   1
[5,]          1  9   0   1
[6,]          1 11   0   1
t(X_perf)
           [,1] [,2] [,3] [,4] [,5] [,6]
Intercepto    1    1    1    1    1    1
X             2    3    5    8    9   11
D_A           1    1    1    0    0    0
D_B           0    0    0    1    1    1
XtX_perf <- t(X_perf) %*% X_perf
print(XtX_perf)
           Intercepto   X D_A D_B
Intercepto          6  38   3   3
X                  38 304  10  28
D_A                 3  10   3   0
D_B                 3  28   0   3

Dado que \(\mathbf{1}_n = \mathbf{D}_A + \mathbf{D}_B\), las columnas son linealmente dependientes y el determinante es exacto a cero:

det_perf <- det(XtX_perf)
cat("Determinante de X'X:", det_perf, "\n")
Determinante de X'X: -1.865175e-14

Al ser \(\det(\mathbf{X}^t\mathbf{X}) = 0\), la matriz no es invertible (solve(XtX_perf) falla) y los parámetros no son identificables por MCO.

Escenario 2: Multicolinealidad Imperfecta (Matriz mal condicionada)

Para solucionar la trampa de la dummy, omitimos la Zona A (\(D_A\)) como categoría de referencia:

\(y_i = \beta_0 + \beta_1 x_i + \gamma_B D_{i,B} + \varepsilon_i\)

\(\gamma_B\) representa la diferencia constante (el ajuste o desplazamiento vertical) en el valor esperado de la variable respuesta \(Y\) entre la Zona B y la Zona A (categoría de referencia), manteniendo constante la edad del vehículo \(X\).

Important

Si una variable categórica tiene \(r\) categorías (niveles), debes incluir exactamente \(r - 1\) variables dummy en la regresión. La categoría omitida se convierte automáticamente en la Categoría de Referencia.

Interpretación del coeficiente \(\gamma_B\)

  • Para la Zona A (\(D_{i,B} = 0\), categoría de referencia):

    \(E(Y \mid X_i, \text{Zona A}) = \beta_0 + \beta_1 X_i + \gamma_B(0) = \mathbf{\beta_0 + \beta_1 X_i}\)

  • Para la Zona B (\(D_{i,B} = 1\)):

    \(E(Y \mid X_i, \text{Zona B}) = \beta_0 + \beta_1 X_i + \gamma_B(1) = \mathbf{(\beta_0 + \gamma_B) + \beta_1 X_i}\)

Restando ambas ecuaciones:

\(E(Y \mid X_i, \text{Zona B}) - E(Y \mid X_i, \text{Zona A}) = \gamma_B\)

  • Efecto de Grupo Purificado: \(\gamma_B\) no es simplemente la diferencia entre el costo promedio total de la Zona B y la Zona A. Es el efecto neto o adicional de estar en la Zona B sobre el costo del siniestro \(Y\), una vez que ya se ha controlado y corregido el efecto producido por la edad del vehículo (\(X\)).

  • Signo del Coeficiente:

    • Si \(\gamma_B > 0\), la Zona B presenta un costo de siniestro sistemáticamente mayor que la Zona A para vehículos de una misma edad.

    • Si \(\gamma_B < 0\), la Zona B representa un menor costo esperado que la Zona A a igualdad de edad del vehículo

Ahora vamos a ilustrar los resultados de la matriz \((X^tX){-1}\) considerando el modelo \(y_i = \beta_0 + \beta_1 x_i + \gamma_B D_{i,B} + \varepsilon_i\)

# Matriz de diseño reducida
X_imp <- cbind(
  Intercepto = 1,
  X = datos$x,
  D_B = datos$D_B
)

X_imp
     Intercepto  X D_B
[1,]          1  2   0
[2,]          1  3   0
[3,]          1  5   0
[4,]          1  8   1
[5,]          1  9   1
[6,]          1 11   1
t(X_imp)
           [,1] [,2] [,3] [,4] [,5] [,6]
Intercepto    1    1    1    1    1    1
X             2    3    5    8    9   11
D_B           0    0    0    1    1    1
XtX_imp <- t(X_imp) %*% X_imp

# Inversa de X'X
inv_XtX <- solve(XtX_imp)
print(round(inv_XtX, 4))
           Intercepto       X     D_B
Intercepto     1.5238 -0.3571  1.8095
X             -0.3571  0.1071 -0.6429
D_B            1.8095 -0.6429  4.5238
det(XtX_imp)
[1] 84

Observa la diagonal principal de \((\mathbf{X}^t\mathbf{X})^{-1}\):

  • \(\text{Var}(\hat{\beta}_0) = 1.5238 \, \sigma^2\)
  • \(\text{Var}(\hat{\beta}_1) = 0.1071 \, \sigma^2\)
  • \(\text{Var}(\hat{\gamma}_B) = \mathbf{4.5238 \, \sigma^2}\)

El valor elevado en la posición \((3,3)\) se debe a que la covariable continua \(x\) y la dummy \(D_B\) están fuertemente asociadas en la muestra (la Zona B agrupa únicamente a los vehículos con \(X \ge 8\)).

Ajustamos el modelo mediante lm() y calculamos el Factor de Inflación de la Varianza (VIF):

modelo_sub <- lm(y ~ x + D_B, data = datos)
summary(modelo_sub)

Call:
lm(formula = y ~ x + D_B, data = datos)

Residuals:
         1          2          3          4          5          6 
 1.166e-15  7.500e-01 -7.500e-01 -1.000e+00  7.500e-01  2.500e-01 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   7.5000     1.1819   6.346  0.00792 **
x             2.2500     0.3134   7.180  0.00557 **
D_B          -2.5000     2.0364  -1.228  0.30711   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9574 on 3 degrees of freedom
Multiple R-squared:  0.9881,    Adjusted R-squared:  0.9802 
F-statistic: 124.8 on 2 and 3 DF,  p-value: 0.001295
# Cálculo de VIF
vif(modelo_sub)
       x      D_B 
6.785714 6.785714

Dado que el \(VIF\) de la variable dummy del grupo B es mayor a \(5\), el modelo de regresión presenta un problema grave de multicolinealidad, por lo que se elimina la variable del modelo.

Una variable categórica con más de 2 niveles (\(k > 2\))

Supongamos que queremos predecir el Sueldo mensual (\(y\)) de un egresado universitario según sus Años de experiencia (\(x\)) y su Nivel máximo de estudios:

  1. Licenciatura

  2. Maestría

  3. Doctorado

Como hay \(r = 3\) niveles, creamos \(r - 1 = 2\) variables dummy y dejamos Licenciatura como categoría base (\(D_1=0, D_2=0\)):

  • \(D_1 = 1\) si tiene Maestría, \(0\) en otro caso.

  • \(D_2 = 1\) si tiene Doctorado, \(0\) en otro caso.

El Modelo

\(y = \beta_0 + \beta_1 x + \beta_2 D_1 + \beta_3 D_2 + \varepsilon\)

Veamos las ecuaciones esperadas para cada grupo:

  • Licenciatura (\(D_1=0, D_2=0\)): \(E[Y] = \beta_0 + \beta_1 X\)

  • Maestría (\(D_1=1, D_2=0\)): \(E[Y] = (\beta_0 + \beta_2) + \beta_1 X\)

  • Doctorado (\(D_1=0, D_2=1\)): \(E[Y] = (\beta_0 + \beta_3) + \beta_1 X\)

Interpretación de Coeficientes:

  • \(\beta_0\): Sueldo base estimado para alguien con Licenciatura y 0 años de experiencia.

  • \(\beta_2\): Diferencia esperada en el sueldo entre quien tiene Maestría vs. quien tiene Licenciatura (misma experiencia).

  • \(\beta_3\): Diferencia esperada en el sueldo entre quien tiene Doctorado vs. quien tiene Licenciatura (misma experiencia).

Nota importante: Si quieres comparar Doctorado vs. Maestría directamante, la diferencia esperada sería \((\beta_3 - \beta_2)\).

Planteamiento Formal de la Prueba Global (Prueba \(F\))

Imaginemos el modelo completo con una variable continua \(x\) (Horas de estudio) y una variable categórica de 3 niveles representada por 2 dummies: \(D_1\) (Maestría) y \(D_2\) (Doctorado):

\(Y = \beta_0 + \beta_1 X + \beta_2 D_1 + \beta_3 D_2 + \varepsilon\)

La prueba de hipótesis global evalúa si al menos uno de los predictores del modelo tiene un efecto distinto de cero:

  • Hipótesis Nula (\(H_0\)): \(\beta_1 = \beta_2 = \beta_3 = 0\)

    (Ninguna variable cuantitativa ni categoría explicativa ayuda a predecir \(Y\)).

  • Hipótesis Alternativa (\(H_1\)): Al menos un \(\beta_j \neq 0\) para \(j \in \{1, 2, 3\}\).

    (Al menos una variable o categoría aporta información significativa para predecir \(Y\)).

El Estadístico de Prueba \(F\)

El estadístico se calcula exactamente igual a través de la Tabla de Análisis de Varianza (ANOVA):

\(F = \frac{\text{MSR}}{\text{MSE}} = \frac{\text{SSR} / k}{\text{SSE} / (n - k - 1)}\)

Donde:

  • \(k\): Número total de parámetros ajustados sin contar el intercepto (en el ejemplo, \(k = 3\): una continua + dos dummies).

  • \(n\): Tamaño de la muestra.

La Prueba \(F\) Parcial (o por bloques)

Aunque la prueba global nos dice si el modelo entero funciona, en modelos con variables dummy suele ser mucho más interesante responder a una pregunta más específica:

¿Aporta algo la variable categórica en su conjunto (con todas sus \(dummies\)), después de haber tomado en cuenta las variables cuantitativas?

No es lo mismo evaluar \(\beta_2 = 0\) individualmente (que solo compara la categoría Maestría contra Licenciatura) que evaluar si la variable “Nivel de Estudios” como un todo es relevante.

Para hacer esto se utiliza una Prueba \(F\) Parcial (o de modelos anidados), comparando:

  1. Modelo Reducido (sin la categórica):

    \(Y = \beta_0 + \beta_1 X + \varepsilon\)

  2. Modelo Completo (con la categórica completa):

    \(Y = \beta_0 + \beta_1 X + \beta_2 D_1 + \beta_3 D_2 + \varepsilon\)

Hipótesis para la prueba del bloque dummy:

  • \(H_0\): \(\beta_2 = \beta_3 = 0\) (El nivel educativo no tiene ningún efecto sobre el sueldo, una vez controlada la experiencia \(X\)).

  • \(H_1\): Al menos un \(\beta_j \neq 0\) para \(j \in \{2, 3\}\). (El nivel educativo sí influye en el sueldo).

El estadístico \(F\) parcial compara la suma de cuadrados de los residuos de ambos modelos:

\(F_{\text{parcial}} = \frac{(\text{SSE}_{\text{reducido}} - \text{SSE}_{\text{completo}}) / q}{\text{MSE}_{\text{completo}}}\)

(Donde \(q\) es el número de variables dummy en el bloque que estamos probando, en este caso \(q = 2\)).

Ejemplo 2

Se realiza un estudio para evaluar el Puntaje final (\(Y\)) obtenido por 100 estudiantes en un examen de Estadística, en función de dos variables:

  1. \(x_1\) (Horas de estudio): Variable cuantitativa continua (horas semanales dedicadas).

  2. \(x_2\) (Tipo de escuela): Variable cualitativa nominal con dos niveles (Pública y Privada).

Nos interesa no solo evaluar si el tipo de escuela afecta el puntaje promedio inicial, sino también determinar si el impacto de cada hora adicional de estudio varía según el tipo de institución (término de interacción)

library(ggplot2)
Warning: package 'ggplot2' was built under R version 4.6.1
# 1. Simulación de datos
set.seed(42)
n <- 100
horas <- runif(n, 1, 20)
escuela <- sample(c("Pública", "Privada"), n, replace = TRUE)

# Creamos el modelo verdadero con interacción
# Puntaje = 40 + 1.5*Horas + 5*DummyPrivada + 0.8*(Horas * DummyPrivada) + error
dummy_privada <- ifelse(escuela == "Privada", 1, 0)
puntaje <- 40 + 1.5 * horas + 5 * dummy_privada + 0.8 * (horas * dummy_privada) + rnorm(n, sd = 4)

datos <- data.frame(puntaje, horas, escuela)

# 2. Ajuste del modelo en R (R crea las dummies automáticamente si es factor)
datos$escuela <- factor(datos$escuela, levels = c("Pública", "Privada")) # Pública como base



modelo_reducido <- lm(puntaje ~ horas, data = datos)
modelo_completo <- lm(puntaje ~ horas + escuela, data = datos)

# Prueba F parcial para la variable categórica "nivel_estudios"
anova(modelo_reducido)
Analysis of Variance Table

Response: puntaje
          Df Sum Sq Mean Sq F value    Pr(>F)    
horas      1  13466 13466.4  211.49 < 2.2e-16 ***
Residuals 98   6240    63.7                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova( modelo_completo)
Analysis of Variance Table

Response: puntaje
          Df  Sum Sq Mean Sq F value    Pr(>F)    
horas      1 13466.4 13466.4  720.51 < 2.2e-16 ***
escuela    1  4427.0  4427.0  236.86 < 2.2e-16 ***
Residuals 97  1812.9    18.7                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modelo_reducido, modelo_completo)
Analysis of Variance Table

Model 1: puntaje ~ horas
Model 2: puntaje ~ horas + escuela
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     98 6240.0                                  
2     97 1812.9  1      4427 236.86 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis del Rendimiento Vehicular (mtcars)

Analizaremos la relación entre el rendimiento del combustible y las características del vehículo, incorporando una variable dummy para el tipo de transmisión.

Queremos modelar el rendimiento de los automóviles en millas por galón (\(y\): mpg) en función de una covariable cuantitativa: el peso del vehículo en miles de libras (\(x\): wt), y una variable cualitativa: el tipo de transmisión (am), donde:

  • am = 0: Transmisión Automática (Grupo de Referencia)

  • am = 1: Transmisión Manual (Variable Dummy \(D_{am}\))

library(knitr)
library(ggplot2)
library(dplyr)
# Copia de mtcars
datos <- mtcars

class(mtcars$am)
[1] "numeric"
# Convertir am en factor con etiquetas legibles
datos$am_factor <- factor(datos$am, levels = c(0, 1), labels = c("Automatico", "Manual"))

# Visualizamos las primeras filas
head(select(datos, mpg, wt, am_factor))
                   mpg    wt  am_factor
Mazda RX4         21.0 2.620     Manual
Mazda RX4 Wag     21.0 2.875     Manual
Datsun 710        22.8 2.320     Manual
Hornet 4 Drive    21.4 3.215 Automatico
Hornet Sportabout 18.7 3.440 Automatico
Valiant           18.1 3.460 Automatico
modelo_mtcars <- lm(mpg ~ wt + am_factor, data = datos)
summary(modelo_mtcars)

Call:
lm(formula = mpg ~ wt + am_factor, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5295 -2.3619 -0.1317  1.4025  6.8782 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     37.32155    3.05464  12.218 5.84e-13 ***
wt              -5.35281    0.78824  -6.791 1.87e-07 ***
am_factorManual -0.02362    1.54565  -0.015    0.988    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.098 on 29 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7358 
F-statistic: 44.17 on 2 and 29 DF,  p-value: 1.579e-09
vif(modelo_mtcars)
       wt am_factor 
 1.921413  1.921413

Análisis del VIF

Dado que ningún valor del \(VIF\) sobrepasa el 5, entonces no se presenta un problema severo de multicolinealidad. El modelo lineal es:

\(\text{mpg}_i = \beta_0 + \beta_1 \text{wt}_i + \gamma \, D_{am,i} + \varepsilon_i\)

  • \(\hat{\beta}_0\) (Intercepto = 37.32): Es la estimación del rendimiento esperado (en millas por galón) para un vehículo automático (am_factor = Automático) cuyo peso (\(wt\)) fuera hipotéticamente de cero.

  • \(\hat{\beta}_1\) (Peso = -5.35): Por cada 1,000 libras adicionales de peso, la estimación del rendimiento de combustible disminuye en promedio 5.35 millas por galón, manteniendo constante el tipo de transmisión.

  • \(\hat{\gamma}\) (Manual = -0.02): A igual peso del vehículo, los autos con transmisión manual tienen una estimación promedio de rendimiento de combustible -0.02 millas por galón mayor en comparación con los vehículos automáticos (el grupo de referencia).

# Generar predicciones en los datos
datos$predichos <- predict(modelo_mtcars)

ggplot(datos, aes(x = wt, y = mpg, color = am_factor)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_line(aes(y = predichos), linewidth = 1.2) +
  labs(
    title = "Rendimiento de Combustible vs Peso del Vehículo",
    subtitle = "Modelo de Regresión Lineal con Variable Dummy de Transmisión",
    x = "Peso del Vehículo (miles de lbs)",
    y = "Rendimiento (mpg)",
    color = "Transmisión"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

modelo_mtcars1 <- lm(mpg ~ wt , data = datos)
summary(modelo_mtcars1)

Call:
lm(formula = mpg ~ wt, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10