Capítulo 8 Regresión lineal

8.1 Resumen

En este capítulo:

  • Aprendemos a estimar modelos de regresión lineal en R usando la función básica lm().
  • Interpretamos y presentamos los resultados usando tablas y figuras.
  • Realizamos diagnósticos a los modelos que estimamos para afinar nuestros análisis. El objetivo final es utilizar técnicas estadísticas para evaluar nuestros modelos teóricos.

Principales conceptos: regresión lineal, coeficiente, significancia estadística, tabla de regresión, predicciones (valores esperados), diagnósticos. Funciones clave: lm(), modelsummary(), ggpredict().

8.1.1 Librerías

Vamos a utilizar las siguientes librerías:

library(tidyverse)
library(lmtest) # diagnosticos de modelos lineales
library(sandwich) # errores estandar robustos
library(ggeffects) # efectos y predicciones en modelos de regresion
library(modelsummary) # tablas de regresion
theme_set(theme_light())

8.1.2 Datos

Debemos descargar los siguientes archivos de datos y guardarlos en la carpeta /data de nuestro proyecto:

  • Datos de ingreso: link. Para descargar, hacer click derecho, “Guardar como…”.

8.2 Modelos de regresión lineal en R

Supongamos que tenemos una pregunta sobre la relación (causal) entre el ingreso de un individuo y su nivel educativo. Nuestra hipótesis es que a mayor número de años de educación, mayor ingreso anual en promedio. Buscamos y encontramos datos al respecto (y los guardamos en el archivo riverside_final.csv, disponible aquí y los cargamos en R:

# cargar datos
riverside <- read_csv("data/riverside_final.csv")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   edu = col_double(),
##   income = col_double(),
##   senior = col_double(),
##   gender = col_double(),
##   party = col_double()
## )

Reorganizamos un poco los datos -recodificamos unas variables- y miramos las primeras filas:

# recodificar variables categóricas
riverside <- riverside %>%
  mutate(
    gender = factor(gender, levels = c(0, 1), labels = c("Mujer", "Hombre")),
    party = factor(party, levels = c(0, 1, 2), labels = c("Rep.", "Dem.", "Ind."))
  )

# imprimir los datos a la consola
head(riverside)
## # A tibble: 6 x 5
##     edu income senior gender party
##   <dbl>  <dbl>  <dbl> <fct>  <fct>
## 1     8  26430      9 Mujer  Dem. 
## 2     8  37449      7 Hombre Rep. 
## 3    10  34182     16 Mujer  Dem. 
## 4    10  25479      1 Mujer  Ind. 
## 5    10  47034     14 Hombre Rep. 
## 6    12  37656     14 Hombre Rep.

Vemos que tenemos información sobre ingreso, nivel educativo, experiencia laboral, género y afiliación partidista para 32 individuos. Podemos empezar a trabajar.

8.3 Análisis exploratorio y visualización

Empezamos haciendo un poco de análisis exploratorio de datos. Primero, veamos un resumen de las variables usando la función summary():

summary(riverside)
##       edu         income          senior         gender    party   
##  Min.   : 8   Min.   :25479   Min.   : 1.00   Mujer :18   Rep.:10  
##  1st Qu.:12   1st Qu.:44513   1st Qu.: 9.75   Hombre:14   Dem.:15  
##  Median :16   Median :55830   Median :15.00               Ind.: 7  
##  Mean   :16   Mean   :53742   Mean   :14.81                        
##  3rd Qu.:20   3rd Qu.:62717   3rd Qu.:20.25                        
##  Max.   :24   Max.   :82726   Max.   :27.00

Usemos la función datasummary_skim() de la librería modelsummary para ver un resumen similar, pero con información adicional:

datasummary_skim(riverside)
Unique (#) Missing (%) Mean SD Min Median Max
edu 13 0 16.0 4.4 8.0 16.0 24.0
income 32 0 53742.1 14553.0 25479.0 55830.0 82726.0
senior 19 0 14.8 6.9 1.0 15.0 27.0

Esto solo nos muestra un resumen de las variables numéricas. Si tenemos variables categóricas (definidas como tipo factor o character en R), lo especificamos en la función:

datasummary_skim(riverside, type = "categorical")
N %
gender Mujer 18 56.2
Hombre 14 43.8
party Rep.  10 31.2
Dem. 15 46.9
Ind. 7 21.9

modelsummary también sirve para ver correlaciones entre variables numéricas:

datasummary_correlation(riverside)
edu income senior
edu 1 . .
income .79 1 .
senior .34 .58 1

Como nuestras principales variables de interés (ingreso y educación) son numéricas, un gráfico de dispersión de ambas es apropiado para evaluar la relación entre ellas:

riverside %>%
  ggplot(aes(edu, income)) +
  geom_point() +
  scale_y_continuous(labels = scales::dollar) +
  labs(x = "Educación (años)", y = "Ingreso anual (USD)")

Parece que hay una relación lineal positiva entre los años de educación y el ingreso anual de un individuo. Sin embargo, después de revisar la literatura, creemos que la experiencia laboral y el género de un individuo también impactan su ingreso. Así que exploramos la relación entre la experiencia profesional (en años) de cada individuo y su ingreso:

riverside %>%
  ggplot(aes(senior, income)) +
  geom_point() +
  scale_y_continuous(labels = scales::dollar) +
  labs(x = "Experiencia laboral (años)", y = "Ingreso anual (USD)")

Vemos un patrón similar. Además, podemos mirar la relación entre género -una variable categórica binaria en nuestros datos- e ingreso:

riverside %>%
  ggplot(aes(gender, income)) +
  geom_boxplot() + # caja y bigote
  scale_y_continuous(labels = scales::dollar) +
  labs(x = "Género", y = "Ingreso anual (USD)")

Para que un análisis juicioso de las causas del ingreso a nivel individual debe considerar más de un aspecto.

8.4 Análisis de regresión

Basados en nuestra teoría y la literatura, creemos que el ingreso de un individuo \(i\) es una función lineal del nivel de educación, la experiencia laboral y el género de dicho individuo. Entonces, la ecuación del modelo de regresión que vamos a estimar es el siguiente:

\[ \hat{\text{Ingreso}}_i = \hat{\alpha} + \hat{\beta}_1 \text{Educación}_i + \hat{\beta}_2 \text{Experiencia}_i + \hat{\beta}_3 \text{Hombre}_i + \hat{u}_i \]

8.4.1 Estimación

En R, estimamos modelos de regresión lineal usando la función lm(). Esta función toma como principal argumento una fórmula de la forma y ~ x + z, donde y es la variable dependiente y x y z son variables independientes.

Estimemos un modelo simple (\(\text{Ingreso}_i = \alpha + \hat{\beta} \text{Educación}_i\)) y guardémoslo como un objeto:

modelo_simple <- lm(income ~ edu, data = riverside)

Así de fácil estimamos nuestra primera regresión lineal. Si queremos explorar los resultados del modelo (los coeficientes estimados, por ejemplo), podemos usar la función summary():

summary(modelo_simple)
## 
## Call:
## lm(formula = income ~ edu, data = riverside)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15808  -5783   2088   5127  18379 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11321.4     6123.2   1.849   0.0743 .  
## edu           2651.3      369.6   7.173 5.56e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8978 on 30 degrees of freedom
## Multiple R-squared:  0.6317, Adjusted R-squared:  0.6194 
## F-statistic: 51.45 on 1 and 30 DF,  p-value: 5.562e-08

Vemos varios elementos, incluyendo coeficientes (Estimate), errores estándar (Std. Error) y p-values (Pr(>|t|)) de cada variable, más unas estadísticas del ajuste del modelo en general (como el \(R^2\)). Estos elementos se pueden extraer individualmente. Por ejemplo, si quiero ver solo los coeficientes:

coef(modelo_simple)
## (Intercept)         edu 
##   11321.379    2651.297

Por ahora, nos vamos a enfocar en estos coeficientes. Recordemos brevemente la interpretación básica de los coeficientes de un modelo de regresión lineal:

  • El intercepto (\(\alpha\)) es el valor esperado de \(Y\) (ingreso, en este caso) cuando todas las variables independientes \(= 0\).
  • El coeficiente de regresión de una variable independiente (\(\beta\)) es el incremento en el valor esperado de \(Y\) (ingreso) asociado a un incremento de una unidad en \(X\) (educación).

Así, basados en este modelo, concluimos que un aumento de 1 año en el nivel educativo de un individuo se asocia con un incremento de aproximadamente 2651 USD anuales en su ingreso.

Podemos estimar modelos de regresión con variables independientes categóricas, por ejemplo binarias (también llamadas variables dummy o variables indicador). Usemos la variable gender que toma dos valores en los datos: “Hombre” y “Mujer”. Estimemos el modelo con la misma función, pero cambiando la variable independiente, y veamos los resultados:

modelo_genero <- lm(income ~ gender, data = riverside)
summary(modelo_genero)
## 
## Call:
## lm(formula = income ~ gender, data = riverside)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23459 -10673   4280   8085  22807 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     48938       3225  15.175 1.29e-15 ***
## genderHombre    10980       4876   2.252   0.0318 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13680 on 30 degrees of freedom
## Multiple R-squared:  0.1446, Adjusted R-squared:  0.1161 
## F-statistic: 5.072 on 1 and 30 DF,  p-value: 0.03179

Recordemos que cuando se trata de variables independientes binarias, la interpretación de los coeficientes cambia levemente:

  • \(\alpha\) es el valor esperado de \(Y\) cuando el valor de \(X\) es la categoría de referencia - en este caso, \(\hat{Y}\) cuando \(X = \text{Mujer}\).
  • \(\beta\) es la diferencia en el valor esperado de \(Y\) para la(s) otra(s) categoría(s) de la variable categórica. En este caso, cuando \(X = 1\) (hombre), \(\hat{Y} = \hat{\alpha} + \hat{\beta}\). Cuando \(X = 0\) (mujer), entonces \(\hat{Y} = \hat{\alpha}\).

Así, basados en este modelo, concluimos que un aumento de 1 año en el nivel educativo de un individuo se asocia con un incremento de aproximadamente 2651 USD anuales en su ingreso.

El uso de variables independientes categóricas binarias es fácilmente extendible a variables categóricas multinominales u ordinales, pues podemos conceptualizarlas como una serie de variables dummy. La e

Veamos qué pasa cuando utilizamos la variable “identificación partidista” (party) como variable independiente. Esta variable toma tres valores: “Rep.”, “Dem.”, o “Ind.”. La ecuación sería:

\[ \hat{\text{Ingreso}}_i = \hat{\alpha} + \hat{\beta}_1 \text{Dem.}_i + \hat{\beta}_2 \text{Ind.}_i + \hat{u}_i \]

Donde “Dem.” e “Ind.” son variables dummy que pueden tomar los valores 0 o 1, para quienes pertecen a ese grupo (1) o no (0). Como antes, la categoría omitida es la de referencia y el intercepto la “absorbe”: \(\hat\alpha\) es el ingreso esperado para los “Rep.”.

En la práctiva, no tenemos que crear estas variable “dummy” manualmente; nuestro software las crea automáticamente. Entonces, estimamos el modelo en R de la siguiente manera y obtenemos estos resultados:

modelo_party <- lm(income ~ party, data = riverside)
summary(modelo_party)
## 
## Call:
## lm(formula = income ~ party, data = riverside)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -35353  -7729   2414   8438  26639 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    52588       4588  11.463 2.73e-12 ***
## partyDem.      -1384       5923  -0.234    0.817    
## partyInd.       8244       7149   1.153    0.258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14510 on 29 degrees of freedom
## Multiple R-squared:  0.07035,    Adjusted R-squared:  0.006236 
## F-statistic: 1.097 on 2 and 29 DF,  p-value: 0.3472

Vemos dos coeficientes para party: uno para individuos que se identifican como “Dem.” y otro para “Ind.”. No hay coeficiente para “Rep.”, lo cual indicica que es la categoría de referencia o base. La interpretación es análoga:

  • \(\alpha\) es el valor esperado de \(Y\) cuando el valor de \(X\) es la categoría de referencia - en este caso, \(\hat{Y}\) cuando \(X = \text{Rep.}\).
  • \(\beta_1\) y \(\beta_2\) son la diferencia en el valor esperado de \(Y\) cuando comparamos las categorías correspondientes con la de referencia. En este caso, cuando \(X = \text{Dem.}\), \(\hat{Y} = \hat{\alpha} + \hat{\beta}_1\). Cuando \(X = \text{Ind.}\) (Rep.), entonces \(\hat{Y} = \hat{\alpha} + \hat{\beta}_2\).

¿Pero podemos decir que estas asociaciones es una buena aproximación al efecto de la educación o del género sobre el ingreso? Sabemos que hay variables omitidas, lo cual introduce sesgo. Podemos mejorar la situación usando modelos de regresión múltiples.

8.4.1.1 Regresión múltiple

Sin embargo, dado nuestro modelo teórico, nos interesa estimar un modelo de regresión múltiple con más de una variable independiente. Para esto, agregamos más variables independientes a la fórmula con el operador +. Estimamos el modelo, lo guardamos como un objeto e imprimimos a la consola un resumen de los resultados:

modelo_ingreso <- lm(income ~ edu + senior + gender + party, data = riverside)
summary(modelo_ingreso)
## 
## Call:
## lm(formula = income ~ edu + senior + gender + party, data = riverside)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11262.9  -4865.4   -854.5   3935.1  12100.1 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5092.8     4750.9   1.072 0.293579    
## edu            2153.4      293.8   7.329 8.78e-08 ***
## senior          707.2      179.0   3.951 0.000531 ***
## genderHombre   8084.2     2586.6   3.125 0.004331 ** 
## partyDem.      -805.0     2939.3  -0.274 0.786348    
## partyInd.      2558.3     3326.3   0.769 0.448760    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6356 on 26 degrees of freedom
## Multiple R-squared:   0.84,  Adjusted R-squared:  0.8093 
## F-statistic: 27.31 on 5 and 26 DF,  p-value: 1.427e-09

8.4.2 Interpretación

Esta información es suficiente para interpretar tres elementos básicos: la magnitud, dirección y significancia estadística de los coeficientes. La única diferencia es que ahora interpretamos los coeficientes “manteniendo todo lo demás constante” (ceteris paribus) o “controlando por” las demás variables. Así, basados en la tabla anterior:

  • El coeficiente edu es el cambio en income asociado a un incremento de 1 unidad en edu, ceteris paribus.
  • El coeficiente senior es el cambio en income asociado a un incremento de 1 unidad en senior, ceteris paribus.
  • El coeficiente genderHombre es la diferencia en income entre hombres y mujeres (la categoría de referencia), ceteris paribus.
  • El coeficiente partyDem. es la diferencia en income entre demócratas y republicanos (la categoría de referencia), ceteris paribus.
  • El coeficiente partyInd. es la diferencia en income entre independientes y republicanos (la categoría de referencia), ceteris paribus.

El signo (positivo o negativo) en cada coeficiente indica si la relación es positiva o negativa. Así, mientras que un aumento en el nivel educativo tiene una asociación positiva con ingreso (si la educación aumenta, el ingreso aumenta), los demócratas tienen menos ingreso, comparados con los republicanos.

El valor p (o p-value) para cada coeficiente nos indica la probabilidad de observar ese coeficiente (\(\hat\beta\), que estimamos) si el coeficiente “real” (\(\beta\), en la población) es igual a 0. Convencionalmente, si es probabilidad es baja -típicamente \(p < 0.05\)- decimos que el coeficiente asociado es estadísticamente significativo.

Volviendo a nuestros resultados, podemos hacer la siguiente interpretación en términos sustanciales que:

  • Manteniendo las demás variables constantes, un año adicional de educación se asocia en promedio con un incremento de 2153 USD en el ingreso anual.
  • Manteniendo las demás variables constantes, un año adicional de experiencia laboral se asocia en promedio con un incremento de 7.07 USD en el ingreso anual.
  • Manteniendo las demás variables constantes, el ingreso anual de un hombre es en promedio 8084 USD mayor que el de una mujer.
  • Manteniendo las demás variables constantes, el ingreso anual de un demócrata es en promedio -805 USD menor que el de un republicano.
  • Manteniendo las demás variables constantes, el ingreso anual de un independiente es en promedio 8084 USD mayor que el de un republicano.

Los tres primeros coeficientes son estadísticamente significativos (\(p < 0.05\)), lo cual indica que la probabilidad de observar estos valores si la hipótesis nula (\(\beta = 0\)) fuese cierta es menos del 5%. O sea, es tan poco probable observar estos coeficientes, que ya que los observamos, debe ser que la hipótesis nula no representa el mundo real. En el mundo en que la hipótesis nula es cierta, \(\beta = 0\) – el p-value nos dice si el \(\beta\) que observamos/estimamos encaja en un mundo en que \(\beta = 0\). No podemos decir lo mismo de los coeficientes asociados a la indentidad partidista; por tanto, no hay una relación clara entre esta variable y el ingreso.

–>

–> –> –> –> –> –> –> –> –> –> –> –> –>

8.4.3 Diagnósticos

El modelo de regresión lineal por mínimos cuadrados ordinarios (MCO u OLS) es el estimador lineal insesgado de mínima varianza (“OLS is BLUE”) cuando se cumplen una serie de supuestos, entre ellos, supuestos sobre los errores \(u\) y su distribución. R nos permite realizar unos diagnósticos para evaluar si estos supuestos se cumplen.

8.4.3.1 Linealidadad

En OLS, asumimos que \(Y\) es una función lineal de las variables independientes. Evaluamos este supuesto visualmente con una gráfica de los valores esperados del modelo (las predicciones) vs. los residuos del mismo, buscando si hay una relación entre ambos:

plot(modelo_ingreso, which = 1)

Si vemos una relación clara (¡miren la línea roja!), hay un problema de linealidad. Violar este supuesto es un problema de específicación del modelo, o sea qué variables incluimos y de qué forma entran. En este caso, podríamos transformar nuestras variables independientes con logaritmos o funciones exponenciales.

8.4.3.2 Normalidad

OLS asume que los residuos (errores) del modelo están distribuidos normalmente. Nuevamente, realizamos una evaluación visual del supuesto (una gráfica Q-Q), buscando que los puntos no se desvíen mucho de la línea diagonal, especialmente en los extremos:

plot(modelo_ingreso, which = 2)

Otra opción es usar un prueba estadística, como el test de Shapiro o el Kolmogorov-Smirnoff, vía las funciones shapiro.test(). Si obtenemos un p-value menor a 0.05, rechazamos la hipótesis nula de normalidad en la distribución de los residuos:

modelo_ingreso %>%
  resid() %>%
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.97118, p-value = 0.5326

Cuando vemos violaciones muy severas de este supuesto, puede que haya un problema con nuestra variable dependiente. Quizás en realidad sea una medida binaria, categórica u ordinal con pocos niveles. En estos casos, podemos buscar medidas alternativas u técnicas distintas a la regresión lineal por OLS.

8.4.3.3 Independencia

Un supuesto muy importante es que las observaciones son sacadas de una muestra de manera independiente; en otras palabras, que no hay patrones de dependencia espacial, temporal o multinivel entre observaciones. Este es un supuesto que podemos evaluar simplemente conociendo nuestros datos y teoría. ¿Mis casos son países que observo año tras año? ¿O quizás espero que lo que sucede en un municipio afecte a los vecinos? Esas son situaciones de autocorrelación temporal y/o espacial. Otra forma de evaluar este supuesto es viendo mapas o series de tiempo de nuestra variable dependiente en búsqueda de tendencias.

Sin embargo, también existen pruebas estadísticas para evaluar el supuesto de independencia. El test de Durbin-Watson, implementado con la función dwtest() de la librería lmtest, es una de estas pruebas:

dwtest(modelo_ingreso)
## 
##  Durbin-Watson test
## 
## data:  modelo_ingreso
## DW = 1.7049, p-value = 0.1433
## alternative hypothesis: true autocorrelation is greater than 0

Un p-value por encima de 0.05 es consistente con la hipótesis nula de no autocorrelación.

Violaciones al supuesto de independencia se solucionan estimando un modelo distinto; en algunas ocasiones, muy distinto. Estos son temas avanzados que no cubrimos aquí, pero más adelante discutimos los estimadores de efectos fijos para analizar datos en panel.

8.4.3.4 Homoesquedasticidad

El supuesto de homoesquedasticidad significa que la varianza de los residuos es constante. Violar este supuesto -los errores están correlacionados entre sí, por ejemplo, y la varianza de los mismos no es constante- implica que los errores estándar estarían mal estimados (más pequeños de lo que deberían ser), dándonos una falsa sensación de confianza.

Para evaluar este supuesto, podemos construir un gráfico de valores esperados de \(Y\) versus los residuos estandarizados (en realidad, la raíz cuadrada de los valores absolutos de estos errores). Cuando se cumple el supuesto, no debemos observar ninguna relación clara entre ambos valores:

plot(modelo_ingreso, which = 3)

Buscamos que a) la línea roja sea aproximadamente horizontal y b) la dispersión de los puntos no cambie mucho en función de los valores esperados. Aquí parece que tenemos un problema de heteroesquedasticidad.

Un test de Breusch-Pagan nos permitiría ponerle más precisión a este diagnóstico. Usamos la función bptest(), nuevamente de la librería lmtest. Veamos:

bptest(modelo_ingreso)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ingreso
## BP = 1.7593, df = 5, p-value = 0.8813

La hipótesis nula en esta prueba es homoesquedasticidad. Si la rechazamos con un p-value bajo (menor a 0.05), tenemos heteroesquedasticidad.

Una de las razones por las que podemos encontrar heteroesquedasticidad en un modelo es que tenemos observaciones con mucha influencia sobre los resultados; podríamos decir que son atípicas. La \(D\) de Cook es una estadística que nos permite medir la influencia de una observación en un modelo de regresión lineal. Vamos a calcular esta estadística para cada una de las observaciones que utilizamos en nuestro modelo y vamos a colocarlas junto a dichas observaciones en un tibble:

datos_modelo <- bind_cols(
  modelo_ingreso$model,
  d_cook = cooks.distance(modelo_ingreso)
)
head(datos_modelo)
##   income edu senior gender party      d_cook
## 1  26430   8      9  Mujer  Dem. 0.002473517
## 2  37449   8      7 Hombre  Rep. 0.006452342
## 3  34182  10     16  Mujer  Dem. 0.007670930
## 4  25479  10      1  Mujer  Ind. 0.062101165
## 5  47034  10     14 Hombre  Rep. 0.005507870
## 6  37656  12     14 Hombre  Rep. 0.089874223

Ahora, podemos buscar y eliminar las observaciones con una \(D\) de Cook alta. Dos reglas informales sugieren que observaciones con valores por encima de 1 o con valores por encima de 3 veces la media pueden tener alta influencia. Vamos a eliminar estas observaciones, reestimar el modelo y volver a evaluar el supuesto de homoesquedasticidad:

datos_modelo %>%
  filter(d_cook < 3*mean(d_cook, na.rm = TRUE)) %>%
  lm(income ~ edu + senior + gender, data = .) %>%
  plot(., which = 3)

Parece que no mejoró mucho la situación.

8.4.3.4.1 Errores estándar “robustos”

Otra opción para lidiar con los problemas asociados con la heteroesquedasticidad es corregir los errores estándar después de estimar el modelo y diagnosticarlo. Después de todo, sabemos que ese es el efecto de violar el supuesto. Vamos a utilizar la función coeftest() de la librería sandwich para recalcular los errores estándar. El resultado son lo que llamamos “erorres estándar robustos”, “errores consistentes cono heteroesquedasticiad” o “errores de Huber-White”:

coeftest(modelo_ingreso, vcov = vcovHC)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   5092.82    4702.82  1.0829 0.2887809    
## edu           2153.40     244.77  8.7977 2.844e-09 ***
## senior         707.22     169.43  4.1741 0.0002965 ***
## genderHombre  8084.16    2493.22  3.2425 0.0032421 ** 
## partyDem.     -804.99    3087.15 -0.2608 0.7963348    
## partyInd.     2558.31    3471.46  0.7370 0.4677484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Noten que los errores estándar de las variables edu y senior ahora son más grandes. La práctica de estimar errores estándar robustos de distintos tipos es bastante común en Ciencia Política y Economía. Más adelante, los incorporamos a nuestra presentación de resultados.

Hay otras soluciones al problema de heteroesquedasticidad, como aplicar weighted least squares (WLS) o usar técnicas como el bootstrap. Por el momento, procedemos con nuestro análisis, estimamos errores robustos si es el caso y reportamos nuestros resultados.

8.5 Predicciones y efectos

Hasta ahora, hemos discutido cómo interpretar los coeficientes de un modelo de regresión simple y múltiple con distintos tipos de variables independientes (numéricas, binarias y multinominales), utilizando la función summary() para ver estos resultados. Pero podemos refinar un poco más la interpretación usando el modelo de regresión como herramienta predictiva, lo cual nos permite ver mejor cambios en la variable dependiente cuando las variables independientes toman distintos valores.

Una de las herramientas más potentes que nos da un modelo de regresión lineal es la capacidad de hacer predicciones. Una predicción no es más que el valor esperado de \(Y\) para determinados valores de las variables independientes. Hallamos estos valores esperados o predicciones utilizando los coeficientes estimados y reemplazando los valores de las variables independientes en la ecuación del modelo de regresión lineal.

Calculemos un valor esperado “a mano” usando nuestro modelo. Primero, vamos a extraer los coeficientes del modelo usando coef() y pluck() - están en el orden en que aparecen en el resumen del modelo:

alfa <- coef(modelo_ingreso) %>% pluck(1)
beta_edu <- coef(modelo_ingreso) %>% pluck(2)
beta_exp <- coef(modelo_ingreso) %>% pluck(3)
beta_hombre <- coef(modelo_ingreso) %>% pluck(4)
beta_dem <- coef(modelo_ingreso) %>% pluck(5)
beta_ind <- coef(modelo_ingreso) %>% pluck(6)

Ahora, vamos a seleccionar valores de las variables independientes que nos interesan:

intercepto <- 1 # siempre es 1
edu <- 10 # 10 años de educacion
exp <- max(riverside$senior, na.rm = TRUE) # ejecutiva
gen <- 0 # mujer
dem <- 1 # demócrata
ind <- 0

Finalmente, reemplazamos para hallar el valor esperado de \(Y\):

intercepto + beta_edu*edu + beta_exp+exp + beta_hombre*gen + beta_dem*dem + beta_ind*ind
## [1] 21464.26

Comparemos con nuestra predicción para un hombre, manteniendo las demás características igual:

gen <- 1
intercepto + beta_edu*edu + beta_exp+exp + beta_hombre*gen + beta_dem*dem + beta_ind*ind
## [1] 29548.42

¡La diferencia entre ambos valores esperados es exactamente el coeficiente de la variable género, 8084.1579866!

Con esta información, podemos hacer inferencias sobre la relación entre las variables de interés. En otras palabras, podemos ver “efectos”, asumiendo que hemos superado los distintos obstáculos a la inferencia causal (especialmente, que no hay sesgo por variables omitidas).

Podemos automatizar, o si se quiere “tercerizar”, este proceso. Digamos que estamos interesados en conocer el ingreso esperado de una mujer con educación y experiencia promedio en comparación con un hombre con las mismas cualificaciones. Vamos a usar la librería ggeffects, específicamente la función ggpredict(). Tomamos el modelo, lo pasamos a la función y especificamos las variables y si queremos, los valores que nos interesan. La función mantiene las variables numéricas no especificadas constantes en la media y mantiene las categóricas en su valor de referencia (es importante entonces recodificar bien):

modelo_ingreso %>%
  ggpredict(terms = "gender")
## # Predicted values of income
## 
## gender | Predicted |               95% CI
## -----------------------------------------
## Mujer  |  50023.02 | [44635.48, 55410.56]
## Hombre |  58107.17 | [53825.13, 62389.22]
## 
## Adjusted for:
## *    edu = 16.00
## * senior = 14.81
## *  party =  Rep.

Vemos entonces el ingreso esperado para hombres y mujeres con 16 años de educación, 14.81 años de experiencia y con afiliación republicana.

O de pronto nos interesa comparar el valor esperado del ingreso de individuos cuando se tiene poca y mucha experiencia. Usamos la opción [minmax] para ver el mínimo y el máximo de la variable independiente de interés:

modelo_ingreso %>%
  ggpredict(terms = "senior[minmax]")
## # Predicted values of income
## 
## senior | Predicted |               95% CI
## -----------------------------------------
##      1 |  40254.50 | [33166.81, 47342.19]
##     27 |  58642.30 | [51619.79, 65664.80]
## 
## Adjusted for:
## *    edu = 16.00
## * gender = Mujer
## *  party =  Rep.

En este caso, vemos el ingreso esperado para individuos con poca y mucha experiencia (1 y 27 años) con las siguientes características: 16 años de mujer, mujeres, republicanas.

Podemos incluir más variables independientes (hasta cuatro en ggeffects):

modelo_ingreso %>%
  ggpredict(terms = c("gender", "senior[minmax]"))
## # Predicted values of income
## 
## # senior = 1
## 
## gender | Predicted |               95% CI
## -----------------------------------------
## Mujer  |  40254.50 | [33166.81, 47342.19]
## Hombre |  48338.66 | [41537.76, 55139.56]
## 
## # senior = 27
## 
## gender | Predicted |               95% CI
## -----------------------------------------
## Mujer  |  58642.30 | [51619.79, 65664.80]
## Hombre |  66726.45 | [61007.34, 72445.57]
## 
## Adjusted for:
## *   edu = 16.00
## * party =  Rep.

Esta es una herramienta potente para interpretar los resultados de un modelo de regresión lineal.

8.6 Presentar resultados

Miremos dos formas principales de presentar los resultados de una regresión: tablas y gráficas. Vamos a utilizar las librerías ggeffects y modelsummary para esto. Hay muchas más opciones, como margins, sjPlot, gtsummary, stargazer, etc.

8.6.1 Tablas de regresión

Casi todos los artículos que utilizan modelos de regresión presentan los resultados en tablas. En ellas, típicamente vemos los coeficientes, intervalos de confianza o errores estándar, la significancia de los coeficientes (con p-values), el número de observaciones y algunas medidas del ajuste del modelo, como el \(R^2\).

Para hacer tablas de regresión, usamos la función modelsummary() de modelsummary. Le pasamos nuestro modelo y la función arroja una tabla formateada:

modelsummary(modelo_ingreso)
Model 1
(Intercept) 5092.825
(4750.857)
edu 2153.403
(293.821)
senior 707.223
(179.005)
genderHombre 8084.158
(2586.555)
partyDem. -804.993
(2939.325)
partyInd. 2558.308
(3326.328)
Num.Obs. 32
R2 0.840
R2 Adj. 0.809
AIC 658.6
BIC 668.9
Log.Lik. -322.311
F 27.308

Los argumentos de la función nos permiten cambiar muchos de los elementos de la tabla, como ponerle etiquetas a los coeficientes, presentar intervalos de confianza en vez de errores estándar y agregar un título:

modelsummary(
  modelo_ingreso,
  coef_map = c(
    "(Intercept)" = "(Intercepto)",
    "edu" = "Educación",
    "senior" = "Experiencia",
    "genderHombre" = "Género: hombre"
  ),
  statistic = "conf.int",
  title = "Resultados del análisis."
)
Tab. 8.1: Resultados del análisis.
Model 1
(Intercepto) 5092.825
[-4672.702, 14858.351]
Educación 2153.403
[1549.445, 2757.361]
Experiencia 707.223
[339.273, 1075.173]
Género: hombre 8084.158
[2767.418, 13400.898]
Num.Obs. 32
R2 0.840
R2 Adj. 0.809
AIC 658.6
BIC 668.9
Log.Lik. -322.311
F 27.308

Una tabla de regresión puede presentar varios modelos lado a lado para efectos de comparación. Si ponemos los modelos que estimamos antes en una lista (y, como bono, les damos nombres), podemos incluirlos en una misma tabla:

lista_modelos <- list(
  "Educación" = modelo_simple,
  "Género" = modelo_genero,
  "Completo" = modelo_ingreso
)
modelsummary(
  lista_modelos, # varios modelos
  coef_map = c(
    "(Intercept)" = "(Intercepto)",
    "edu" = "Educación",
    "senior" = "Experiencia",
    "genderHombre" = "Género: hombre",
    "partyDem." = "Afil.: dem.",
    "partyInd." = "Afil.: ind."
  )
)
Educación Género Completo
(Intercepto) 11321.379 48938.167 5092.825
(6123.235) (3224.912) (4750.857)
Educación 2651.297 2153.403
(369.623) (293.821)
Experiencia 707.223
(179.005)
Género: hombre 10980.476 8084.158
(4875.609) (2586.555)
Afil.: dem. -804.993
(2939.325)
Afil.: ind. 2558.308
(3326.328)
Num.Obs. 32 32 32
R2 0.632 0.145 0.840
R2 Adj. 0.619 0.116 0.809
AIC 677.3 704.3 658.6
BIC 681.7 708.7 668.9
Log.Lik. -335.655 -349.137 -322.311
F 51.452 5.072 27.308

Hay varias opciones adicionales. Por ejemplo, podemos indicar niveles de significancia estadística con la opción stars = TRUE. Adicionalmente, si queremos reportar errores estándar robustos, está la opción vcov = ... que trabaja en conjunto con las funciones de la librería sandwich. Así, para errores estándar robustos “básicos”:

modelsummary(
  lista_modelos, # varios modelos
  coef_map = c(
    "(Intercept)" = "(Intercepto)",
    "edu" = "Educación",
    "senior" = "Experiencia",
    "genderHombre" = "Género: hombre",
    "partyDem." = "Afil.: dem.",
    "partyInd." = "Afil.: ind."
  ),
  stars = TRUE,
  vcov = "robust"
)
## Warning: In version 0.8.0 of the `modelsummary` package, the default significance markers produced by the `stars=TRUE` argument were changed to be consistent with R's defaults.
## This warning is displayed once per session.
Educación Género Completo
(Intercepto) 11321.379+ 48938.167*** 5092.825
(5880.924) (3217.157) (4702.821)
Educación 2651.297*** 2153.403***
(350.350) (244.769)
Experiencia 707.223***
(169.432)
Género: hombre 10980.476* 8084.158**
(5087.413) (2493.216)
Afil.: dem. -804.993
(3087.151)
Afil.: ind. 2558.308
(3471.459)
Num.Obs. 32 32 32
R2 0.632 0.145 0.840
R2 Adj. 0.619 0.116 0.809
AIC 677.3 704.3 658.6
BIC 681.7 708.7 668.9
Log.Lik. -335.655 -349.137 -322.311
F 51.452 5.072 27.308
Std. Errors Robust Robust Robust
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Si, en cambio, tenemos datos en panel (muchas unidades and varios periodos de tiempo), podemos usar errores estándar robustos por cluster o grupo. Miremos rápidamente con unos datos simulados en los cuales observamos 500 empresas (firm) a lo largo de diez años year, midiendo unas variables x y y. Como creemos que los residuos/errores pueden estar correlacionados entre empresas y años, incluimos una corrección con base en las variables relevantes:

data("PetersenCL")
mod_firmas <- lm(x ~ y, data = PetersenCL)
modelsummary(
  mod_firmas,
  # esta formula indica que variables vamos a usar para los clusters
  vcov = ~ firm + year
)
Model 1
(Intercept) -0.002
(0.030)
y 0.201
(0.009)
Num.Obs. 5000
R2 0.208
R2 Adj. 0.208
AIC 12952.1
BIC 12971.6
Log.Lik. -6473.041
F 1310.740
Std. Errors C: firm + year

Por último, podemos exportar una tabla a un documento Word:

modelsummary(
  modelo_ingreso,
  output = "output/tabla_modelo_ingreso.docx"
)

8.6.2 Gráficas

Visualizar los resultados de modelos, especialmente las predicciones y la incertidumbre de estos, puede ayudar mucho a entenderlos mejor y comunicar nuestros resultados. Vamos a usar la función ggpredict de ggeffects. También podríamos usar sjPlot o hacerlo completamente “a mano”.

Como mencionamos arriba, podemos usar ggpredict() para ver los “efectos” de distintas variables, entendidos como el cambio en la variable dependiente para distintos valores de las independientes. Para esto, especificamos las variables independientes que nos interesan con terms =. Esta función mantiene constantes los factores (variables categóricas) en su categoría de referencia y las variables numéricas según su media.

Le pasamos el modelo a la función `ggpredict()y luego pasamos los resultados de esta a plot(). A continuación, vemos el valor esperado (predicción) de la variable dependiente (ingreso anual) en todo el rango de la variable educación, controlando por género (categoría de referencia: mujer), experiencia (constante en la media) y afiliación partidista (republicano). Como ggeffects funciona con ggplot2, podemos modificar y agregar utilizando los trucos que ya conocemos:

modelo_ingreso %>%
  ggpredict(terms = "edu") %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de educación",
       subtitle = "Mujeres, controlando por experiencia laboral",
       x = "Educación (años)", y = "Ingreso (USD)")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Adicionalmente, podemos especificar dos variables independientes para visualizar aún más variación. Vemos cómo la “línea base” de ingreso (el intercepto) es distinto para hombres y mujeres. Así, vemos más claramente la diferencia en ingreso para hombres y mujeres, por ejemplo. Además, esto evidencia que incluir una variable categórica esencialmente crea un intercepto distinto para cada grupo, sin cambiar la pendiente (el coeficiente):

modelo_ingreso %>%
  ggpredict(terms = c("edu", "gender")) %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de educación y experiencia", 
       x = "Educación (años)", y = "Ingreso (USD)", color = "Género")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Si las dos variables son numéricas, ggpredict() automáticamente escoge tres valores del rango de la segunda variable que le pasamos (incluyendo la media) y grafica tres líneas:

modelo_ingreso %>%
  ggpredict(terms = c("edu", "senior")) %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de educación y experiencia",
       x = "Educación (años)", y = "Ingreso (USD)", color = "Experiencia")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Podemos espeficiar cuántas líneas queremos ver y los valores que debe tomar la segunda variable independiente:

modelo_ingreso %>%
  ggpredict(terms = c("edu", "senior [5, 25]")) %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de educación y experiencia",
       x = "Educación (años)", y = "Ingreso (USD)", color = "Experiencia")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

Finalmente, podemos cambiar el orden de las variables en terms = si nos interesa resaltar el efecto de una variable más que la otra:

modelo_ingreso %>%
  ggpredict(terms = c("gender", "edu [5, 15, 25]")) %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de género y educación",
       x = "Género", y = "Ingreso (USD)", color = "Educación")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

modelo_ingreso %>%
  ggpredict(terms = c("senior", "edu [5, 25]")) %>%
  plot() +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Valores esperados de ingreso para distintos valores de educación y experiencia",
       x = "Experiencia (años, log)", y = "Ingreso (USD)", color = "Educación")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

8.7 Conclusiones

Utilidad de los modelos de regresión lineal:

  • Predicción: si cumplimos con los supuestos básicos, podemos hacer predicciones puntuales precisas basados en la asociación entre variables.
  • Explicación: si cumplimos con más supuestos (evitamos problemas de endogeneidad y de sesgo por variables omitidas), podemos hacer inferencias sobre relaciones de tipo causal.

Otras posibilidades:

  • Variables dependientes categóricas: regresión logística y modelos por máxima verosimilitud (MLE o maximum likelihood estimation).
  • Experimentos: podemos usar OLS para encontrar la diferencia de medias.
  • Data Science: regresión como modelo predictivo que se puede “entrenar”.

Utilizar R para el análisis estadístico:

  • Curva de aprendizaje, pero grandes posibilidades.
  • Practicar, practicar y practicar.
  • Google es nuestro amigo.

8.8 Taller: regresión lineal

Encontrar una base de datos relevante para el proyecto de investigación. Limpiar y ordenar los datos. El código de esta entrega no debe incluir todo el procedimiento para cargar la base de datos original y el proceso de limpieza y organización de la misma.

8.8.1 Estimación

Estimar dos modelos de regresión lineal múltiple usando la función lm() de R. Los modelos estimados deben cumplir con los siguientes requisitos (1.0 punto):

  • La variable dependiente en ambos debe ser la variable dependiente del proyecto de investigación. Idealmente esta debe ser una variable de tipo continuo, pero excepcionalmente puede ser de tipo ordinal (si tiene \(\geq 7\) categorías ordenadas) o binaria/dummy (entendiendo que la regresión lineal no es el mejor modelo para este tipo de variables). En estos dos casos, deben asegurarse que las variables estén codificadas (temporalmente) como tipo numérico o dbl en R.
  • Las variables independientes deben ser a) factores potencialmente explicativos de la variable dependiente y b) variables de control. La selección de variables debe estar basada en su revisión de literatura, teoría e hipótesis.
  • Debe haber por lo menos una variable independiente numérica en uno de los dos modelos. Si no tienen una, deben buscar más datos.
  • Debe haber por lo menos una variable independiente categórica binaria en uno de los dos modelos. Si no tienen una, deben crearla a partir de otra.
  • Dado el punto 5 de este taller, es recomendable que por lo menos uno de los modelos tenga ambos tipos de variables.
  • Sugerencia: Si la literatura sugiere un debate entre dos enfoques teóricos, pueden estimar un modelo que represente cada enfoque (por ejemplo, un modelo solo con variables institucionales y otro solo con variables sociológicas), seguido de un modelo con todas las variables.
  • O pueden construir el modelo completo de manera progresiva: un modelo con solo variables independientes de interés teórico, seguido de uno que añade variables de control.

8.8.2 Comunicación de resultados

Presentar los resultados de los modelos estimados en una sola tabla de regresión, usando las funciones de la librería modelsummary (o alguna similar) para producirla, teniendo en cuenta las siguientes indicaciones (1.0 punto).

  • La tabla debe tener título y los coeficientes deben estar nombrados apropiadamente (o sea, con algo distinto al código de la variable en la base de datos).
  • Usar errores estándar robustos o robustos por cluster, dependiendo del tipo de datos (corte transversal o panel).
  • Nota: Pueden ejecutar ?modelsummary en la consola para explorar todas las opciones de esta funciones y cómo utilizarlas.

8.8.3 Diagnóstico

Evaluar si los modelos estimados cumplen el supuesto de homoesquedasticidad (\(\text{Var}(e_i) = \sigma^2\)) y si hay observaciones con un grado alto de influencia. Esta evaluación puede hacerse de manera visual o gráfica o usando tests estadísticos. Deben incluir un breve párrafo en el que consignen sus apreciaciones al respecto (1.0 punto).

8.8.4 Interpretación

Interpretar sustantivamente todos los coeficientes de los modelos incluídos en la tabla y la incertidumbre en torno a estos estimados (pueden usar p-values o intervalos de confianza). No es necesario interpretar el coeficiente del intercepto. Cuando decimos “sustantivamente” nos referimos a que las interpretaciones deben hacer especial énfasis en la forma en que los resultados de los modelos se conectan con la evaluación de las hipótesis propuestas en el proyecto (1.0 punto).

8.8.5 Interpretación gráfica

Usando uno de los dos modelos de regresión estimados en el primer punto, construir una gráfica con ggpredict() donde presenten valores esperados de la variable dependiente para el rango de valores de una variable independiente numérica y para todas las categorías de una variable categórica binaria (manteniendo constantes las demás variables). En otras palabras, la gráfica debe mostrar los valores esperados de la variable dependiente (con intervalos de confianza) para dos o más grupos (definidos por la variable categórica) a lo largo del rango de una variable independiente numérica (1.0 punto).

8.9 Ejercicios adicionales

8.9.1 Votos por Trump

Digamos que nos interesa estudiar los resultados de las elecciones presidenciales de 2016 en Estados Unidos. Para este proyecto, utilizamos una base de datos que incluye información sobre un conjunto de variables sociales y políticas a nivel de condado (municipio). Para esto, vamos a cargar la librería socviz que cuenta con datos sociopolíticos de Estados Unidos:

library(socviz) # datos y herramientas para visualizar datos sociales
## 
## Attaching package: 'socviz'
## The following object is masked _by_ '.GlobalEnv':
## 
##     edu

Específicamente, nos interesa la base de datos county_data, la cual limpiamos y reducimos:

# seleccionamos y transformamos datos
county_data_sub <- county_data %>% 
  as_tibble() %>%
  select(name, state, census_region, black, white, hh_income,
         per_dem_2016, per_gop_2016, partywinner12) %>%
  mutate(partywinner12 = fct_explicit_na(partywinner12),
         black = black/100, 
         white = white/100, 
         hh_income = log(hh_income))

sample_n(county_data_sub, 10)
## # A tibble: 10 x 9
##    name      state census_region black white hh_income per_dem_2016 per_gop_2016
##    <chr>     <fct> <fct>         <dbl> <dbl>     <dbl>        <dbl>        <dbl>
##  1 Kenai Pe~ AK    West          0.009 0.846      11.0        0.377        0.529
##  2 Allen Pa~ LA    South         0.232 0.718      10.6        0.228        0.743
##  3 St. Jose~ IN    Midwest       0.131 0.815      10.7        0.477        0.475
##  4 Wayne Co~ KY    South         0.018 0.963      10.3        0.179        0.797
##  5 Jefferso~ NE    Midwest       0.005 0.975      10.7        0.242        0.696
##  6 Terry Co~ TX    South         0.05  0.921      10.5        0.240        0.735
##  7 Warren C~ OH    Midwest       0.035 0.903      11.2        0.289        0.665
##  8 Mendocin~ CA    West          0.01  0.866      10.7        0.604        0.312
##  9 Atlantic~ NJ    Northeast     0.173 0.714      10.9        0.515        0.453
## 10 Bedford ~ TN    South         0.082 0.877      10.6        0.221        0.749
## # ... with 1 more variable: partywinner12 <fct>

Veamos el resumen de algunas de estas variables. La tabla a continuación (construida con datasummary() de modelsummary) muestra la media y desviación estándar (variables numéricas) y el número de observaciones por grupo (variables categóricas):

# estadisticas resumen
datasummary(
  # variables numericas
  (`Prop. afro` = black) +
    (`Prop. blanco` = white) +
    (`Ingreso promedio (log.)` = hh_income) +
    (`Prop. voto Dem.` = per_dem_2016) +
    (`Prop. voto Rep.` = per_gop_2016) ~ census_region *
    # funciones
    ((`Media` = Mean) + (`Desv. est.` = SD)),
  data = county_data_sub
)
Midwest
Northeast
South
West
Media Desv. est. Media Desv. est. Media Desv. est. Media Desv. est.
Prop. afro 0.03 0.05 0.06 0.07 0.17 0.18 0.02 0.02
Prop. blanco 0.93 0.10 0.89 0.10 0.79 0.18 0.86 0.16
Ingreso promedio (log.) 10.75 0.17 10.90 0.23 10.62 0.26 10.78 0.24
Prop. voto Dem. 0.29 0.12 0.44 0.15 0.32 0.16 0.34 0.17
Prop. voto Rep.  0.66 0.12 0.52 0.15 0.65 0.16 0.57 0.17

Podemos mirar las frecuencias de las variables categóricas:

datasummary(
  # variable categorica
  (`Partido 2012` = partywinner12) ~ 
    # funcion
    N,
  data = county_data_sub
)
Partido 2012 N
Democrat 686
Republican 2455
(Missing) 54

También una tabla cruzada con dos variables categóricas y porcentajes además de frecuencias usando el operador *:

datasummary(
  # variables categoricas
  (`Región` = census_region) * 
    (`Partido 2012` = partywinner12) + 1 ~ 
    # funciones
    N + Percent(),
  data = county_data_sub
)
Región Partido 2012 N Percent
Midwest Democrat 191 5.98
Republican 864 27.04
(Missing) 12 0.38
Northeast Democrat 124 3.88
Republican 93 2.91
(Missing) 9 0.28
South Democrat 251 7.86
Republican 1171 36.65
(Missing) 18 0.56
West Democrat 120 3.76
Republican 327 10.23
(Missing) 14 0.44
All 3195 100.00

O combinar variables numéricas y categóricas en una tabla de estadísticas descriptivas por grupos, usando nuevamente *:

datasummary(
  # variables numericas
  (`Prop. afro` = black) +
    (`Prop. blanco` = white) +
    (`Ingreso promedio (log.)` = hh_income) +
    (`Prop. voto Dem.` = per_dem_2016) +
    (`Prop. voto Rep.` = per_gop_2016) ~ 
    # variables categoricas
    census_region *
    # funciones
    ((`Media` = Mean) + (`Desv. est.` = SD)),
  data = county_data_sub
)
Midwest
Northeast
South
West
Media Desv. est. Media Desv. est. Media Desv. est. Media Desv. est.
Prop. afro 0.03 0.05 0.06 0.07 0.17 0.18 0.02 0.02
Prop. blanco 0.93 0.10 0.89 0.10 0.79 0.18 0.86 0.16
Ingreso promedio (log.) 10.75 0.17 10.90 0.23 10.62 0.26 10.78 0.24
Prop. voto Dem. 0.29 0.12 0.44 0.15 0.32 0.16 0.34 0.17
Prop. voto Rep.  0.66 0.12 0.52 0.15 0.65 0.16 0.57 0.17

Queremos entender la variación en la proporción de votos obtenidos por el candidato Demócrata (Donald Trump) en la elección de 2016 en cada condado. A continuación, vemos la distribución de esta variable:

county_data %>%
  ggplot(aes(per_gop_2016)) +
  geom_histogram() +
  labs(x = "Prop. voto Rep.", y = "Número de condados")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 54 rows containing non-finite values (stat_bin).

Nuestra teoría sugiere que esperamos ver una mayor proporción de votos por Trump en condados con una mayor proporción de individuos que se identifican como blancos. La siguiente gráfica de dispersión explora esta relación:

county_data %>%
  ggplot(aes(white, per_gop_2016)) +
  geom_point(alpha = 0.5) +
  labs(x = "Prop. blanco", y = "Prop. voto Rep.")
## Warning: Removed 54 rows containing missing values (geom_point).

Para evaluar esta hipótesis, estimamos un modelo de regresión lineal múltiple. La variable dependiente es la proporción de votos por el candidato Republicano en cada condado. La variable independiente de interés es la proporción de habitantes blancos en cada condado. Además, incluimos varias variables de control: el logaritmo natural del ingreso promedio (en USD) de los hogares en cada condado, una dummy que indica el partido ganador en 2012 por condado y la región donde se ubica el condado. La ecuación que estimamos es:

\[ \hat{\text{Prop. voto Rep.}} = \alpha + \hat{\beta}_1 \times \text{Prop. blanco}_i + \hat{\beta}_2 \times \text{log(Ingreso)}_i + \hat{\beta}_3 \times \text{Partido 2012}_i + \hat{\beta}_4 \times \text{Región}_i + \hat{u}_i \]

La siguiente tabla muestra los resultados del modelo. Incluye el coeficiente (\(\beta\)) estimado para cada variable y el intercepto, así como intervalos de confianza del 95% entre paréntesis y la significancia estadística de cada coeficiente, además de unas estadísticas resumen del modelo (como el \(\text{R}^2\)).

modelo_rep <- lm(
  per_gop_2016 ~ white + hh_income + partywinner12 + census_region,
  data = county_data
)

coefs <- c(
    "(Intercept)" = "Intercepto", 
    "white" = "Prop. Blanco", 
    "hh_income" = "Ingreso prom. (log.)",
    "partywinner12Republican" = "Partido 2012: Rep.", 
    "census_regionNortheast" = "Región: Nordeste", 
    "census_regionSouth" = "Región: Sur",
    "census_regionWest" = "Región: Oeste"
)

modelsummary(
  list(`Modelo 1` = modelo_rep), 
  stars = TRUE,
  coef_map = coefs
)
Modelo 1
Intercepto 0.256***
(0.011)
Prop. Blanco 0.004***
(0.000)
Ingreso prom. (log.) 0.000***
(0.000)
Partido 2012: Rep.  0.218***
(0.004)
Región: Nordeste -0.032***
(0.007)
Región: Sur 0.030***
(0.004)
Región: Oeste -0.043***
(0.005)
Num.Obs. 3141
R2 0.702
R2 Adj. 0.702
AIC -6544.1
BIC -6495.7
Log.Lik. 3280.056
F 1232.355
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

8.9.1.1 Interpretación

Para cada coeficiente, interpretar:

  • Dirección.
  • Magnitud.
  • Significancia estadística.

8.9.1.2 Predicciones

Con la ecuación y los resultados, podemos hacer predicciones sobre el valor esperado de \(Y\) (proporción de votos por Trump) para distintos valores de las variables independientes. ¿Cuál es el valor esperado de \(Y\) para los siguientes perfiles de condados?

  • 80% blanco, en la región Sur, en un condado en el que los Republicanos ganaron en 2012- las demás variables en sus medias (ver la tabla de estadísticas descriptivas).
  • 20% blanco, en el Nordeste, en un condado en el que los Demócratas ganaron en 2012- nuevamente, establecemos las demás variables en sus medias.

8.9.2 Duración de gobiernos democráticos

Algunos politólogos argumentan que la variación en la duración de gobiernos democráticos en el periodo posguerra se explica por características de los gobiernos de turnos: el nivel de apoyo parlamentario de los partidos en el gobierno, el número de partidos en el mismo y la disciplina que las organizaciones partidistas pueden imponer a sus miembros. En otras palabras, tenemos el siguiente modelo teórico de la duración de un gobierno:

\[ \text{Duración gob.} = f(\text{Apoyo parl. gob.} + \text{Núm. partidos gob.} + \text{Disciplina part.} + \text{Comp. estocástico}) \]

8.9.2.1 Datos

Como estudiantes de Ciencia Política, buscamos datos al respecto para evaluar este argumento y nos encontramos con una base de datos recopilada por Rob Franzese hace unos años a partir del Political Data Handbook of OECD Countries con datos para países de la OCDE. La base de datos se encuentra en el archivo OECD_country_data.xls aquí.

  • Cargue los datos en R usando las funciones y librerías apropiadas para un archivo .xls.

Afortunadamente, el archivo también incluye un diccionario de variables (en inglés) en la segunda pestaña del archivo. Pista: read_excel() tiene un argumento sheet =.

  • Identifique las variables relevantes para estimar el modelo especificado arriba y asegúrese que estén en el formato adecuado.

8.9.2.2 Regresión

La información contenida en esta base de datos nos permite decir algo sobre el argumento teórico sugerido al inicio.

  • Utilizando estos datos, estime un modelo de regresión lineal sencillo (sin utilizar variables transformadas) que evalúe el argumento, presente los resultados en una tabla de regresión e interpretelos de forma sustancial (en términos del argumento).

8.9.2.3 Críticas

Frente a estos resultados, una colega sugiere que lo que está sucediendo es un poco distinto. Ella dice: “la disciplina partidista no tiene un efecto independiente sobre la duración de los gobiernos. Más bien, condiciona el efecto del apoyo legislativo al gobierno y del número de partidos en la coalición sobre la duración”.

  • Especifique un modelo que tome en cuenta esta crítica, estímelo via OLS e interprete los efectos de las variables independientes sobre la dependiente:

8.9.2.4 Extra: comparar modelos

¿Cómo adjudicamos entre estos modelos? Otra colega sugiere que una prueba F, realizada por medio de un análisis de varianza (ANOVA) de los dos modelos, podría decirnos algo sobre cuál es mejor.

  • Realice una prueba F y concluya cuál modelo prefiere en términos estadísticos y teóricos:

8.9.2.5 Comunicar resultados

Es útil visualizar los resultados de modelos estadísticos como estos. Para el siguiente punto, seleccione el modelo que crea que se ajusta mejor a los datos y representa mejor la teoría.

  • Construya una gráfica que muestre el efecto del número de partidos en la coalición de gobierno y del apoyo legislativo al gobierno sobre la duración de los mismos cuando la disciplina partidista es alta y cuando es baja. Incluya intervalos de confianza en la gráfica.

8.9.3 Opinión pública sobre el presidente

Supongamos que queremos entender las percepciones y evaluaciones de individuos sobre los presidentes, particularmente los de USA. Por inspiración de las diosas politológicas, sabemos que el modelo “real” es:

\[ \text{Eval. presid.} = \alpha + \beta_1\text{Afroamericano} + \beta_2\text{ID partidista} + \beta_3\text{Ideología} + \beta_4\text{Eval. idiosincrática} + \beta_5\text{Eval. sociotrópica} + \epsilon \]

Esta es una versión de un modelo clásico en Ciencia Política, el cual considera que las evaluaciones individuales de los políticos son una función de posiciones ideológicos y percepciones de la economía. Queremos evaluar este modelo.

8.9.3.1 Datos

Nuevamente, buscamos datos y encontramos la encuesta del American National Election Study (ANES) para 2004, en la cual le preguntaron a los encuestados por su percepción del entonces presidente George W. Bush. Los datos se encuentran en el archivo nes_2004_data.csv aquí. Las principales variables son:

Nombre de variable Definición
thermom Sentimientos por Bush (1:100; 0=mal, 100=bien)
race Raza/etnicidad (10, 12, 13, 14 y 15=afroamericano)
partyid Identificación partidista (0=muy Demócrata, 3=independiente, 6=muy Republicano, 7=otro)
libcon ¿Liberal o conservador? (1=liberal, 3=moderado, 5=conservador, 7=no responde)
betteroff Situación económica familiar el último año (1-5; 1=peor, 3=igual, 5=mejor)
economy Situación económica del país el último año (1-5; 1=peor, 3=igual, 5=mejor)
edulevel Nivel educativo máximo (0-7; 0=sin bachillerato, 7=posgrado)
  • Cargue el archivo de datos y realice los procedimientos necesarios para asegurarse que las variables estén debidamente codificadas.

8.9.3.2 Modelos

Aprovechemos que sabemos que el modelo de arriba es el “real”.

  • Estime el modelo “real”, interprete los coeficientes y diga algo sobre la interpretación causal de los mismos.

  • Estime el siguiente modelo -llamémoslo “modelo 2”-, que tiene una variable menos y donde los coeficientes son estimados (\(\hat{\beta}\) en vez de \(\beta\)):

\[ \text{Eval. presid.} = \alpha + \hat{\beta_1}\text{Afroamericano} + \hat{\beta_2}\text{ID partidista} + \hat{\beta_3}\text{Eval. idiosincrática} + \hat{\beta_4}\text{Eval. sociotrópica} + \epsilon \]

  • Estime este otro modelo -llamémoslo algo creativo, como “modelo 3”-, que tiene una variable más y donde los coeficientes son estimados (\(\hat{\beta}\) en vez de \(\beta\)):

\[ \text{Eval. presid.} = \alpha + \hat{\beta_1}\text{Afroamericano} + \hat{\beta_2}\text{ID partidista} + \hat{\beta_3}\text{Ideología} + \hat{\beta_4}\text{Eval. idiosincrática} + \hat{\beta_5}\text{Eval. sociotrópica} + \hat{\beta_6}{Educación} + \epsilon \]

8.9.3.3 Extra: presentar resultados

  • Presente los resultados de los tres modelos en una tabla y saque conclusiones sobre el sesgo de los coeficientes en el modelo 2 y la eficiencia del modelo 3.

  • Presente gráficamente y discuta los efectos de todas las variables independientes en uno de los modelos.