library(tidyverse) # ggplot(), %>%, mutate(), and friends
library(broom) # Convert models to data frames
library(scales) # Format numbers with functions like comma(), percent(), and dollar()
library(modelsummary) # Create side-by-side regression tablesEjemplop de Diff-in-Diff
Antecedentes del programa
En 1980, Kentucky aumentó el límite de ingresos semanales cubiertos por la compensación laboral. Queremos saber si esta nueva política provocó que los trabajadores pasaran más tiempo desempleados. Si las prestaciones no son lo suficientemente generosas, los trabajadores podrían demandar a las empresas por lesiones laborales, mientras que unas prestaciones excesivamente generosas podrían causar problemas de riesgo moral e inducir a los trabajadores a ser más imprudentes en el trabajo o a alegar lesiones fuera del trabajo.
La principal variable de resultado que nos interesa es log_duration(en los datos originales como ldurat, pero la renombramos para que sea más legible), o la duración registrada (en semanas) de las prestaciones de compensación laboral. La registramos porque la variable está bastante sesgada: la mayoría de las personas están desempleadas durante algunas semanas, y algunas lo están durante un tiempo prolongado. La política se diseñó para que el aumento del límite no afectara a los trabajadores con bajos ingresos, pero sí a los de altos ingresos, por lo que utilizamos a los trabajadores con bajos ingresos como grupo de control y a los de altos ingresos como grupo de tratamiento.
Los datos se incluyen en el paquete R {wooldridge} como injuryconjunto de datos. Si instala el paquete, lo carga con [nombre del archivo library(wooldridge)] y lo ejecuta ?injuryen la consola, podrá ver todos sus detalles. Para que practique la carga de datos desde archivos externos, exporté los datos de lesiones como un archivo CSV (usando [nombre del archivo write_csv(injury, "injury.csv")]) y los incluí aquí.
Éstas son las columnas principales de las que nos ocuparemos por ahora:
durat(que renombraremos comoduration): Duración de los beneficios de desempleo, medida en semanasldurat(que renombraremos alog_duration): Versión registrada dedurat(log(durat))after_1980(que renombraremos comoafter_1980): Variable indicadora que marca si la observación ocurrió antes (0) o después (1) del cambio de política en 1980. Este es nuestro tiempo (o variable antes/después )highearn: Variable indicadora que marca si la observación tiene un ingreso bajo (0) o alto (1). Esta es nuestra variable de grupo (o de tratamiento/control ).
Cargar y limpiar datos
Primero, descarguemos la base de datos:
# Load the data.
# It'd be a good idea to click on the "injury_raw" object in the Environment
# panel in RStudio to see what the data looks like after you load it
injury_raw <- read_csv("data/injury.csv")A continuación, limpiaremos injury_rawun poco los datos para limitarlos solo a las observaciones de Kentucky. También cambiaremos el nombre de algunas columnas con nombres extraños:
injury <- injury_raw %>%
filter(ky == 1) %>%
# The syntax for rename is `new_name = original_name`
rename(duration = durat, log_duration = ldurat,
after_1980 = afchnge)Análisis exploratorio de datos
En primer lugar, podemos observar la distribución de los beneficios de desempleo entre las personas con ingresos altos y bajos (nuestros grupos de control y tratamiento):
ggplot(data = injury, aes(x = duration)) +
# binwidth = 8 makes each column represent 2 months (8 weeks)
# boundary = 0 make it so the 0-8 bar starts at 0 and isn't -4 to 4
geom_histogram(binwidth = 8, color = "white", boundary = 0) +
facet_wrap(vars(highearn))
La distribución está muy sesgada: la mayoría de las personas en ambos grupos reciben entre 0 y 8 semanas de beneficios (¡y un puñado con más de 180 semanas! ¡eso es 3,5 años!).
Si utilizamos el logaritmo de la duración, podemos obtener una distribución menos sesgada que funciona mejor con los modelos de regresión:
ggplot(data = injury, mapping = aes(x = log_duration)) +
geom_histogram(binwidth = 0.5, color = "white", boundary = 0) +
# Uncomment this line if you want to exponentiate the logged values on the
# x-axis. Instead of showing 1, 2, 3, etc., it'll show e^1, e^2, e^3, etc. and
# make the labels more human readable
# scale_x_continuous(labels = trans_format("exp", format = round)) +
facet_wrap(vars(highearn))
También deberíamos comprobar la distribución del desempleo antes y después del cambio de política. Copiar y pegar uno de los fragmentos del histograma y modificar las facetas:
ggplot(data = injury, mapping = aes(x = log_duration)) +
geom_histogram(binwidth = 0.5, color = "white", boundary = 0) +
facet_wrap(vars(after_1980))
Las distribuciones parecen más normales, pero no podemos apreciar fácilmente ninguna diferencia entre los grupos antes/después y los de tratamiento/control. Sin embargo, podemos graficar los promedios. Hay varias maneras de hacerlo.
Puedes usar una stat_summary()capa para que ggplot calcule estadísticas de resumen, como promedios. Aquí simplemente calculamos la media:
ggplot(injury, aes(x = factor(highearn), y = log_duration)) +
geom_point(size = 0.5, alpha = 0.2) +
stat_summary(geom = "point", fun = "mean", size = 5, color = "red") +
facet_wrap(vars(after_1980))
Pero también podemos calcular la media y el intervalo de confianza del 95%:
ggplot(injury, aes(x = factor(highearn), y = log_duration)) +
stat_summary(geom = "pointrange", size = 1, color = "red",
fun.data = "mean_se", fun.args = list(mult = 1.96)) +
facet_wrap(vars(after_1980))
¡Ya podemos empezar a ver el clásico gráfico de diferencias! Parece que, después de 1980, las personas con altos ingresos tuvieron, en promedio, un desempleo más prolongado.
También podemos usar group_by()y summarize()para calcular las medias de los grupos antes de enviar los datos a ggplot. Prefiero esto porque me da más control sobre los datos que estoy graficando:
plot_data <- injury %>%
# Make these categories instead of 0/1 numbers so they look nicer in the plot
mutate(highearn = factor(highearn, labels = c("Low earner", "High earner")),
after_1980 = factor(after_1980, labels = c("Before 1980", "After 1980"))) %>%
group_by(highearn, after_1980) %>%
summarize(mean_duration = mean(log_duration),
se_duration = sd(log_duration) / sqrt(n()),
upper = mean_duration + (1.96 * se_duration),
lower = mean_duration + (-1.96 * se_duration))
ggplot(plot_data, aes(x = highearn, y = mean_duration)) +
geom_pointrange(aes(ymin = lower, ymax = upper),
color = "darkgreen", size = 1) +
facet_wrap(vars(after_1980))
O, graficado en el formato más estándar diff-in-diff:
ggplot(plot_data, aes(x = after_1980, y = mean_duration, color = highearn)) +
geom_pointrange(aes(ymin = lower, ymax = upper), size = 1) +
# The group = highearn here makes it so the lines go across categories
geom_line(aes(group = highearn))
Diferencia en diferencia a mano
Podemos encontrar esa diferencia exacta completando la tabla 2x2 antes/después del tratamiento/control:
| Before 1980 | After 1980 | ∆ | |
|---|---|---|---|
| Low earners | A | B | B − A |
| High earners | C | D | D − C |
| ∆ | C − A | D − B | (D − C) − (B − A) |
Una combinación de group_by()y summarize()hace que esto sea realmente fácil:
diffs <- injury %>%
group_by(after_1980, highearn) %>%
summarize(mean_duration = mean(log_duration),
# Calculate average with regular duration too, just for fun
mean_duration_for_humans = mean(duration))
diffs
## # A tibble: 4 × 4
## # Groups: after_1980 [2]
## after_1980 highearn mean_duration mean_duration_for_humans
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0 1.13 6.27
## 2 0 1 1.38 11.2
## 3 1 0 1.13 7.04
## 4 1 1 1.58 12.9Podemos sacar cada uno de estos números de la tabla con algunos filter()s y pull():
before_treatment <- diffs %>%
filter(after_1980 == 0, highearn == 1) %>%
pull(mean_duration)
before_control <- diffs %>%
filter(after_1980 == 0, highearn == 0) %>%
pull(mean_duration)
after_treatment <- diffs %>%
filter(after_1980 == 1, highearn == 1) %>%
pull(mean_duration)
after_control <- diffs %>%
filter(after_1980 == 1, highearn == 0) %>%
pull(mean_duration)
diff_treatment_before_after <- after_treatment - before_treatment
diff_treatment_before_after
## [1] 0.2
diff_control_before_after <- after_control - before_control
diff_control_before_after
## [1] 0.0077
diff_diff <- diff_treatment_before_after - diff_control_before_after
diff_diff
## [1] 0.19La estimación de diferencias en diferencias es de 0,19, lo que significa que el programa provoca un aumento de la duración del desempleo de 0,19 semanas registradas. Sin embargo, el uso de semanas registradas no tiene sentido, por lo que debemos interpretarlo con porcentajes (¡ aquí tienes una guía práctica!; este es el Ejemplo B, donde se registra la variable dependiente/de resultado). Recibir el tratamiento (es decir, tener ingresos altos después del cambio de política) provoca un aumento del 19 % en la duración del desempleo.
ggplot(diffs, aes(x = as.factor(after_1980),
y = mean_duration,
color = as.factor(highearn))) +
geom_point() +
geom_line(aes(group = as.factor(highearn))) +
# If you use these lines you'll get some extra annotation lines and
# labels. The annotate() function lets you put stuff on a ggplot that's not
# part of a dataset. Normally with geom_line, geom_point, etc., you have to
# plot data that is in columns. With annotate() you can specify your own x and
# y values.
annotate(geom = "segment", x = "0", xend = "1",
y = before_treatment, yend = after_treatment - diff_diff,
linetype = "dashed", color = "grey50") +
annotate(geom = "segment", x = "1", xend = "1",
y = after_treatment, yend = after_treatment - diff_diff,
linetype = "dotted", color = "blue") +
annotate(geom = "label", x = "1", y = after_treatment - (diff_diff / 2),
label = "Program effect", size = 3)
# Here, all the as.factor changes are directly in the ggplot code. I generally
# don't like doing this and prefer to do that separately so there's less typing
# in the ggplot code, like this:
#
# diffs <- diffs %>%
# mutate(after_1980 = as.factor(after_1980), highearn = as.factor(highearn))
#
# ggplot(diffs, aes(x = after_1980, y = avg_durat, color = highearn)) +
# geom_line(aes(group = highearn))Diferencias en diferencias con regresión
Calcular todos los elementos a mano de esa manera es tedioso, así que podemos usar la regresión. Recuerde que debemos incluir variables indicadoras para el tratamiento/control y para el antes/después, así como la interacción entre ambos. Así es como se ve la ecuación matemática:
\[ \begin{aligned} \log(\text{duration}) = &\alpha + \beta \ \text{highearn} + \gamma \ \text{after_1980} + \\ & \delta \ (\text{highearn} \times \text{after_1980}) + \epsilon \end{aligned} \]
El coeficiente \(\delta\) es el efecto que nos importa al final: el estimador de diferencias en diferencias.
model_small <- lm(log_duration ~ highearn + after_1980 + highearn * after_1980,
data = injury)
tidy(model_small)
## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.13 0.0307 36.6 1.62e-263
## 2 highearn 0.256 0.0474 5.41 6.72e- 8
## 3 after_1980 0.00766 0.0447 0.171 8.64e- 1
## 4 highearn:after_1980 0.191 0.0685 2.78 5.42e- 3El coeficiente para highearn:after_1980es el mismo que el que obtuvimos manualmente, ¡como debe ser! Es estadísticamente significativo, por lo que podemos estar bastante seguros de que no es 0.
Diferencias en diferencias con regresión + controles
Una ventaja de usar la regresión para las diferencias en diferencias es que podemos incluir variables de control para aislar el efecto. Por ejemplo, las reclamaciones de los trabajadores de la construcción o la manufactura tienden a tener una duración mayor que las de los trabajadores de otros sectores. O quizás quienes reclaman lesiones de espalda tienden a tener reclamaciones más largas que quienes reclaman lesiones en la cabeza. También podríamos querer controlar los datos demográficos de los trabajadores, como el género, el estado civil y la edad.
Estimemos una versión ampliada del modelo de regresión básico con las siguientes variables adicionales:
malemarriedagehosp(1 = hospitalizado)indust(1 = fabricación, 2 = construcción, 3 = otros)injtype(1-8; categorías para diferentes tipos de lesiones)lprewage(registro de salario antes de presentar una reclamación)
Importante : industy injtypeaparecen en el conjunto de datos como números (1-3 y 1-8), pero en realidad son categorías. Debemos indicarle a R que los trate como categorías (o factores); de lo contrario, asumirá que puede tener un tipo de lesión de 3.46 o algo imposible.
# Convert industry and injury type to categories/factors
injury_fixed <- injury %>%
mutate(indust = as.factor(indust),
injtype = as.factor(injtype))model_big <- lm(log_duration ~ highearn + after_1980 + highearn * after_1980 +
male + married + age + hosp + indust + injtype + lprewage,
data = injury_fixed)
tidy(model_big)
## # A tibble: 18 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.53 0.422 -3.62 2.98e- 4
## 2 highearn -0.152 0.0891 -1.70 8.86e- 2
## 3 after_1980 0.0495 0.0413 1.20 2.31e- 1
## 4 male -0.0843 0.0423 -1.99 4.64e- 2
## 5 married 0.0567 0.0373 1.52 1.29e- 1
## 6 age 0.00651 0.00134 4.86 1.19e- 6
## 7 hosp 1.13 0.0370 30.5 5.20e-189
## 8 indust2 0.184 0.0541 3.40 6.87e- 4
## 9 indust3 0.163 0.0379 4.32 1.60e- 5
## 10 injtype2 0.935 0.144 6.51 8.29e- 11
## 11 injtype3 0.635 0.0854 7.44 1.19e- 13
## 12 injtype4 0.555 0.0928 5.97 2.49e- 9
## 13 injtype5 0.641 0.0854 7.51 7.15e- 14
## 14 injtype6 0.615 0.0863 7.13 1.17e- 12
## 15 injtype7 0.991 0.191 5.20 2.03e- 7
## 16 injtype8 0.434 0.119 3.65 2.64e- 4
## 17 lprewage 0.284 0.0801 3.55 3.83e- 4
## 18 highearn:after_1980 0.169 0.0640 2.64 8.38e- 3
# Extract just the diff-in-diff estimate
diff_diff_controls <- tidy(model_big) %>%
filter(term == "highearn:after_1980") %>%
pull(estimate)Tras controlar diversos controles demográficos, la estimación de diferencias en diferencias es menor (0,17), lo que indica que la política provocó un aumento del 17 % en la duración de las semanas de desempleo tras una lesión laboral. Es menor porque las demás variables independientes ahora explican parte de la variación en log_duration.
Comparación de resultados
Podemos colocar los coeficientes del modelo uno al lado del otro para comparar el valor de highearn:after_1980 cuando cambiamos el modelo.
modelsummary(list("Simple" = model_small, "Full" = model_big))| Simple | Full | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 1.126*** | -1.528*** |
| (0.031) | (0.422) | |
| highearn | 0.256*** | -0.152+ |
| (0.047) | (0.089) | |
| after_1980 | 0.008 | 0.050 |
| (0.045) | (0.041) | |
| highearn × after_1980 | 0.191** | 0.169** |
| (0.069) | (0.064) | |
| male | -0.084* | |
| (0.042) | ||
| married | 0.057 | |
| (0.037) | ||
| age | 0.007*** | |
| (0.001) | ||
| hosp | 1.130*** | |
| (0.037) | ||
| indust2 | 0.184*** | |
| (0.054) | ||
| indust3 | 0.163*** | |
| (0.038) | ||
| injtype2 | 0.935*** | |
| (0.144) | ||
| injtype3 | 0.635*** | |
| (0.085) | ||
| injtype4 | 0.555*** | |
| (0.093) | ||
| injtype5 | 0.641*** | |
| (0.085) | ||
| injtype6 | 0.615*** | |
| (0.086) | ||
| injtype7 | 0.991*** | |
| (0.191) | ||
| injtype8 | 0.434*** | |
| (0.119) | ||
| lprewage | 0.284*** | |
| (0.080) | ||
| Num.Obs. | 5626 | 5347 |
| R2 | 0.021 | 0.190 |
| R2 Adj. | 0.020 | 0.187 |