Regresión Discontínua
Francisco Villarroel
2022-12-08
Conceptualización sobre los diseños de Regresión Discontínua (RDD)
Las regresiones discontínuas se formularon en los años 60’s a partir de los trabajos de Donald Campbell, psicologo educacional. Si bien realizó numerosos trabajos con esta metodología, no tuvo mayor alcance hasta hace 20 años atrás.
La popularidad del método aumentó con los trabajos de Angrist y Lavy (1999) sobre los efectos de los tamaños de los cursos usando una condición particular de los colegios en Israel. La razón por la cual las regresiones discontínuas no tuvieron mucho impacto en el comienzo es porque en econometría primero se tuvo que aceptar el potencial outcome framework, es decir se tuvo que montar un aparataje teórico y conceptual para abrir espacios a la inferencia causal antes de que se pudiese aceptar la represión discontínua como un método válido y valioso.
Representación gráfica del RDD
En el DAG puesto a la izquierda tenemos una variable X que es contínua y corresponde a unidades de tratamiento asignadas a un tratamiento D, que impactan en una variable dependiente Y. Además observamos que hay una covariable U sobre Y y que además puede afectar a X. Para resolver este problema de cofounder de X se asigna un “cutoff” o nivel de umbral donde desde cierto punto en adelante podrá ser asignado a tratamiento y bajo etse cutoff será control. Eso puede ser desde el nivel de grados alcohólicos para ser detenido, puntaje en una prueba para obtener beneficios, puntaje de clasificación social para postulara beneficios del Estado, etc.
Esta variale X con cutoff se le denomina como Running variable. Lo importante de la regresión discontínua para el análisis causal, según cunningham, es que “podemos identificar efectos causales en cuyos sujetos que puntúan en un vecindario cercano certa de un cutoff \(c_o\). Espcíficamente, como mostraremos, el efecto causal promedio de esta subpoblación que es identificada como \(X \to c_0\)”.
Los supuestos principales de las regresiones discontínuas son los siguientes:
- Supuesto de Contínuidad: formalmente significa que “los resultados potenciales esperados son continuos en el punto de corte. Si los resultados potenciales esperados son continuos en el punto de corte, se excluyen necesariamente las intervenciones que compiten entre sí al mismo tiempo.”
De esta forma en el DAG de la derecha la relación \(X \to Y\) está eliminada pues el cutoff la controla. En opinión de Cunningham la continuidad siempre debería ser la hipótesis nula y la alternativa mostrar algún tipo de discontinuidad. “Si ves una tortuga en el poste, sabes que no ha llegado ahí por sí sola”
Al saber la probabilidad de ser asignado al tratamiento, se calculan las LATE o Local Average Treatment effect. Al no haber “soporte común” se trabaja en base a extrapolación, lo que significa que estamos comparando valores diferentes en la variable de ejecución. En este sentdido, “todos los métodos utilizados por la RDD son formas de manejar el sesgo de la extrapolación de la forma más limpia posible”
Ejemplo: “The Close Election”
En este ejemplo de replicación se utilizará una característica específica de las elecciones de Estados Unidos: los ganadores son declarados como tal cuando el candidato tiene un mínimo de votos. Entonces, “en la medida que carreras muy cercanas representen asignaciones exógenas en una victoria partidista, que vamos a discutir después, que uno puede elegir esas elecciones cercanas para”identificar el efecto causal en el ganador en una variedad de resultados”
A partire de estas premisas y para testear elementos clásicos de la economìa política, se analizará si los votantes pueden afectar la política pública. En esto existen dos visiones: la teoría de la convergencia y de la divergencia; la primera señala que la heterogeneidad ideológica de los votante fuerza a cada candidato a moderar su posición para tener éxito, mientras que la teoría de la divergencia, el cual en la iomposibilidad de generar commitment, los candidados simplemente persiguen sus propias políticas públicas de su preferencia.
library(tidyverse)
library(haven)
library(estimatr)
lmb_data <-read_dta("data.dta")
lmb_subset <- lmb_data %>%
filter(lagdemvoteshare>.48 & lagdemvoteshare<.52)
lm_1 <- lm_robust(score ~ lagdemocrat, data = lmb_subset, clusters = id)
lm_2 <- lm_robust(score ~ democrat, data = lmb_subset, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat, data = lmb_subset, clusters = id)
modelsummary::msummary(list(lm_1,lm_2,lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 31.196*** | 18.747*** | 0.242*** |
| (1.334) | (0.843) | (0.020) | |
| lagdemocrat | 21.284*** | 0.484*** | |
| (1.951) | (0.029) | ||
| democrat | 47.706*** | ||
| (1.356) | |||
| Num.Obs. | 915 | 915 | 915 |
| R2 | 0.115 | 0.578 | 0.235 |
| AIC | 8795.5 | 8117.3 | 1088.6 |
| BIC | 8809.9 | 8131.8 | 1103.0 |
| RMSE | 29.49 | 20.36 | 0.44 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
En esta tabla vemos que los resultados sin similares a los de Lee, Moretti y Butler, sin embargo son distintos en otros sentidos. Lo importnate acà es que se usaron regresiones cercanas al cutoff para estimar el efecto. Estas son regresiones locales en el sentido que usan datos cercanos al cutoff,es decir estamos estimando el coeficiente de \(D_t\) a la derecha del cutoff.
En cambio su usamos todos los datos obtenemos lo siquiente:
#using all data (note data used is lmb_data, not lmb_subset)
lmb_data <-read_dta("data.dta")
lm_1 <- lm_robust(score ~ lagdemocrat, data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat, data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat, data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 23.539*** | 17.576*** | 0.120*** |
| (0.337) | (0.263) | (0.004) | |
| lagdemocrat | 31.506*** | 0.818*** | |
| (0.484) | (0.005) | ||
| democrat | 40.763*** | ||
| (0.418) | |||
| Num.Obs. | 13588 | 13588 | 13588 |
| R2 | 0.227 | 0.376 | 0.676 |
| AIC | 129786.7 | 126881.1 | 3902.1 |
| BIC | 129809.3 | 126903.7 | 3924.7 |
| RMSE | 28.69 | 25.78 | 0.28 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Cuando utilizamos todos los datos vemos que variable lagdemocrat sube 10 puntos respecto de la anterior, pero los efectos en el momento mismo (variable “democrat”) se vuelve más pequeño, finalmente la incumbencia crece fuertemente.
A continuación se hará la misma regresión incluyendo los controles para la variable independiente y además se agregará la centralización de la running variable aplicando una corrección de -0.5, para que el centro se convierta en 0, así los valores negativos son demócratas y los positivos son republicanos.
lmb_data <-read_dta("data.dta")
lmb_data <- lmb_data %>%
mutate(demvoteshare_c = demvoteshare - 0.5)
lm_1 <- lm_robust(score ~ lagdemocrat + demvoteshare_c, data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat + demvoteshare_c, data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat + demvoteshare_c, data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 22.883*** | 11.034*** | 0.212*** |
| (0.443) | (0.336) | (0.005) | |
| lagdemocrat | 33.451*** | 0.552*** | |
| (0.848) | (0.010) | ||
| demvoteshare_c | −5.626** | −48.938*** | 0.773*** |
| (1.898) | (1.642) | (0.019) | |
| democrat | 58.502*** | ||
| (0.656) | |||
| Num.Obs. | 13577 | 13577 | 13577 |
| R2 | 0.227 | 0.424 | 0.735 |
| AIC | 129674.5 | 125683.0 | 1153.2 |
| BIC | 129704.5 | 125713.1 | 1183.3 |
| RMSE | 28.68 | 24.76 | 0.25 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
lmb_data <- lmb_data %>%
mutate(demvoteshare_c = demvoteshare - 0.5)
lm_1 <- lm_robust(score ~ lagdemocrat + demvoteshare_c, data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat + demvoteshare_c, data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat + demvoteshare_c, data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 22.883*** | 11.034*** | 0.212*** |
| (0.443) | (0.336) | (0.005) | |
| lagdemocrat | 33.451*** | 0.552*** | |
| (0.848) | (0.010) | ||
| demvoteshare_c | −5.626** | −48.938*** | 0.773*** |
| (1.898) | (1.642) | (0.019) | |
| democrat | 58.502*** | ||
| (0.656) | |||
| Num.Obs. | 13577 | 13577 | 13577 |
| R2 | 0.227 | 0.424 | 0.735 |
| AIC | 129674.5 | 125683.0 | 1153.2 |
| BIC | 129704.5 | 125713.1 | 1183.3 |
| RMSE | 28.68 | 24.76 | 0.25 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Los resultados siguen siendo bastante similares a los de Lee, Moretti y Butler, pero los efectos siguen siendo muy diferentes.
Ahora necesitamos una línea que marque las discontinuidades, lo cual se hace como una interacción entre la running variable con el tratamiento asignado.
lm_1 <- lm_robust(score ~ lagdemocrat*demvoteshare_c,
data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat*demvoteshare_c,
data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat*demvoteshare_c,
data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 31.435*** | 16.816*** | 0.287*** |
| (0.541) | (0.419) | (0.008) | |
| lagdemocrat | 30.508*** | 0.526*** | |
| (0.817) | (0.010) | ||
| demvoteshare_c | 66.042*** | −5.683* | 1.403*** |
| (3.161) | (2.611) | (0.044) | |
| lagdemocrat × demvoteshare_c | −96.475*** | −0.849*** | |
| (3.853) | (0.049) | ||
| democrat | 55.431*** | ||
| (0.637) | |||
| democrat × demvoteshare_c | −55.152*** | ||
| (3.219) | |||
| Num.Obs. | 13577 | 13577 | 13577 |
| R2 | 0.267 | 0.434 | 0.749 |
| AIC | 128963.7 | 125441.1 | 443.1 |
| BIC | 129001.3 | 125478.7 | 480.7 |
| RMSE | 27.94 | 24.54 | 0.25 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
En esta tabal vemos los resultados con las interacciones de variables que, si bien redujo los coeficientes, sigue siendo mayor a los que se utilizan sòlo en el intervalo de 0.05 a 0.5. Por lo tanto ahora se realiza la regresión con una variable cuadrática
lmb_data <- lmb_data %>%
mutate(demvoteshare_sq = demvoteshare_c^2)
lm_1 <- lm_robust(score ~ lagdemocrat*demvoteshare_c + lagdemocrat*demvoteshare_sq,
data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat*demvoteshare_c + democrat*demvoteshare_sq,
data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat*demvoteshare_c + lagdemocrat*demvoteshare_sq,
data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 33.547*** | 15.606*** | 0.330*** |
| (0.714) | (0.575) | (0.013) | |
| lagdemocrat | 13.030*** | 0.322*** | |
| (1.286) | (0.018) | ||
| demvoteshare_c | 134.977*** | −23.850*** | 2.798*** |
| (9.786) | (6.713) | (0.196) | |
| demvoteshare_sq | 212.127*** | −41.729** | 4.294*** |
| (22.763) | (14.686) | (0.456) | |
| lagdemocrat × demvoteshare_c | 57.055*** | 0.091 | |
| (15.412) | (0.241) | ||
| lagdemocrat × demvoteshare_sq | −641.851*** | −8.804*** | |
| (31.331) | (0.517) | ||
| democrat | 44.402*** | ||
| (0.909) | |||
| democrat × demvoteshare_c | 111.896*** | ||
| (9.781) | |||
| democrat × demvoteshare_sq | −229.954*** | ||
| (19.546) | |||
| Num.Obs. | 13577 | 13577 | 13577 |
| R2 | 0.371 | 0.456 | 0.822 |
| AIC | 126894.8 | 124919.4 | −4226.4 |
| BIC | 126947.4 | 124972.0 | −4173.8 |
| RMSE | 25.89 | 24.07 | 0.21 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La inclusión de la variable cuadrática disminuye considerablemente el coeficiente. Así, el efecto sobre el voto contemporáneo es menor que el de encuentran Lee, Moretti y Butler.
El punto es que seguimos estimando regresiones globales y el coefieciente sigue siendo mayor, lo que implicaría que hay outliers muy fuertes que hacen que los coeficiente suban. Por esto se va a liminar denuevo el anàlisis a una ventana más pequeña.
lmb_data <- lmb_data %>%
filter(demvoteshare > .45 & demvoteshare < .55) %>%
mutate(demvoteshare_sq = demvoteshare_c^2)
lm_1 <- lm_robust(score ~ lagdemocrat*demvoteshare_c + lagdemocrat*demvoteshare_sq,
data = lmb_data, clusters = id)
lm_2 <- lm_robust(score ~ democrat*demvoteshare_c + democrat*demvoteshare_sq,
data = lmb_data, clusters = id)
lm_3 <- lm_robust(democrat ~ lagdemocrat*demvoteshare_c + lagdemocrat*demvoteshare_sq,
data = lmb_data, clusters = id)
modelsummary::msummary(list(lm_1, lm_2, lm_3), stars = TRUE)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 37.121*** | 21.437*** | 0.418*** |
| (0.969) | (1.819) | (0.013) | |
| lagdemocrat | 7.347*** | 0.167*** | |
| (1.587) | (0.020) | ||
| demvoteshare_c | 830.925*** | 450.846** | 15.699*** |
| (20.956) | (161.352) | (0.228) | |
| demvoteshare_sq | 5333.335*** | 7878.904** | 91.607*** |
| (838.325) | (2995.192) | (10.893) | |
| lagdemocrat × demvoteshare_c | −156.876*** | 0.125 | |
| (35.740) | (0.357) | ||
| lagdemocrat × demvoteshare_sq | −10116.678*** | −188.329*** | |
| (1435.130) | (16.351) | ||
| democrat | 45.191*** | ||
| (2.679) | |||
| democrat × demvoteshare_c | −688.343** | ||
| (247.711) | |||
| democrat × demvoteshare_sq | −3887.820 | ||
| (4802.371) | |||
| Num.Obs. | 2387 | 2387 | 2387 |
| R2 | 0.445 | 0.563 | 0.774 |
| AIC | 21860.2 | 21290.4 | −74.4 |
| BIC | 21900.6 | 21330.9 | −34.0 |
| RMSE | 23.50 | 20.86 | 0.24 |
| Std.Errors | by: id | by: id | by: id |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
vemos entonces que ahora tenemos una muestra mayor a la del análisis inicial, y que la interacción de la variable cuadrática han disminuído el tamado estiamdo del voto futuro, aún cuando la muestra es mayor.
A modo de resumen: Se ajustó un modelo sin running variable, luego se agregó la running variable ontrudicida de varias formas (interactuándo con otras, como cuadrática, etc), y luego de ello se estiman los LATE de cada uno.
Al agregar controles como la running variable a la probabilida de voto demócrata, vemos que tenemos información diferente a lo de Lee, Moretti y Butler. Ese coeficiente de 21 inicial se atenúa, aunque sigue siendo significativo.
Lo importante ahora es la realización de métodos no-paramétricos, que en este caso significa la estimación un modelo que no asume una forma funcional para la relación entre la variable de resultado (\(Y\)) y la running variable (\(X\)). Esto se formaliza de la siguiente forma:
\[Y = f(X) + \varepsilon\]
Para revisar esto hacemos el gráfico correspondiente:
#aggregating the data
lmb_data <-read_dta("data.dta")
categories <- lmb_data$lagdemvoteshare
demmeans <- split(lmb_data$score, cut(lmb_data$lagdemvoteshare, 100)) %>%
lapply(mean) %>%
unlist()
agg_lmb_data <- data.frame(score = demmeans, lagdemvoteshare = seq(0.01,1, by = 0.01))
#plotting
lmb_data <- lmb_data %>%
mutate(gg_group = case_when(lagdemvoteshare > 0.5 ~ 1, TRUE ~ 0))
ggplot(lmb_data, aes(lagdemvoteshare, score)) +
geom_point(aes(x = lagdemvoteshare, y = score), data = agg_lmb_data) +
stat_smooth(aes(lagdemvoteshare, score, group = gg_group), method = "lm",
formula = y ~ x + I(x^2)) +
xlim(0,1) + ylim(0,100) +
geom_vline(xintercept = 0.5)
Vemos que este gráfico es bastante similar al del paper original.
Podemos probar otros ajustes como el modelo lineal, aunque sólo sea por fines de muestra:
ggplot(lmb_data, aes(lagdemvoteshare, score)) +
geom_point(aes(x = lagdemvoteshare, y = score), data = agg_lmb_data) +
stat_smooth(aes(lagdemvoteshare, score, group = gg_group), method = "loess") +
xlim(0,1) + ylim(0,100) +
geom_vline(xintercept = 0.5)
ggplot(lmb_data, aes(lagdemvoteshare, score)) +
geom_point(aes(x = lagdemvoteshare, y = score), data = agg_lmb_data) +
stat_smooth(aes(lagdemvoteshare, score, group = gg_group), method = "lm") +
xlim(0,1) + ylim(0,100) +
geom_vline(xintercept = 0.5)
Al ver este ajuste lineal vemos la gran influencia que tienen los outliers en las dos regresiones lineales, lo que hace que la influencia del punto de corte esté sobreestimado.
Es importante decir que si no hay una tendencia clara en lso polinomios, investigaciones igual se serias simplemente reportan un ajuste lineal. Por esto en ocasiones se ha recomentado hacer regresiones lineares locales o también regresiones polinómicas ponderadas por kernel.
library(tidyverse)
library(stats)
smooth_dem0 <- lmb_data %>%
filter(democrat == 0) %>%
select(score, demvoteshare)
smooth_dem0 <- as_tibble(ksmooth(smooth_dem0$demvoteshare, smooth_dem0$score,
kernel = "box", bandwidth = 0.1))
smooth_dem1 <- lmb_data %>%
filter(democrat == 1) %>%
select(score, demvoteshare) %>%
na.omit()
smooth_dem1 <- as_tibble(ksmooth(smooth_dem1$demvoteshare, smooth_dem1$score,
kernel = "box", bandwidth = 0.1))
ggplot() +
geom_smooth(aes(x, y), data = smooth_dem0) +
geom_smooth(aes(x, y), data = smooth_dem1) +
geom_vline(xintercept = 0.5)
Otro de los elementos importantes es la elección del ancho de banda para estimar el efecto causal que puede tener un tratamiento. En palabras de Cunningham “A la hora de elegir el ancho de banda, siempre hay un equilibrio entre el sesgo y la varianza: cuanto más corta es la ventana, menor es el sesgo, pero al tener menos datos, aumenta la varianza de la estimación”.
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 13577
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 5480 8097
## Eff. Number of Obs. 2112 1893
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 0.086 0.086
## BW bias (b) 0.141 0.141
## rho (h/b) 0.609 0.609
## Unique Obs. 2770 3351
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 46.491 1.241 37.477 0.000 [44.060 , 48.923]
## Robust - - 31.425 0.000 [43.293 , 49.052]
## =============================================================================Un último elemento relevante es que las regresiones discontínuas funcionan bien cuando se tiene una gran cantidad de datos. Se aplica una prueba de densidad de McCrary para estimación de densidad polinómica local para observar si ha habido manipulación de la running variable.
library(tidyverse)
library(rddensity)
library(rdd)
DCdensity(lmb_data$demvoteshare, cutpoint = 0.5)
## [1] 0.5925909density <- rddensity(lmb_data$demvoteshare, c = 0.5)
rdplotdensity(density, lmb_data$demvoteshare)
## $Estl
## Call: lpdensity
##
## Sample size 5480
## Polynomial order for point estimation (p=) 2
## Order of derivative estimated (v=) 1
## Polynomial order for confidence interval (q=) 3
## Kernel function triangular
## Scaling factor 0.40357984678845
## Bandwidth method user provided
##
## Use summary(...) to show estimates.
##
## $Estr
## Call: lpdensity
##
## Sample size 8097
## Polynomial order for point estimation (p=) 2
## Order of derivative estimated (v=) 1
## Polynomial order for confidence interval (q=) 3
## Kernel function triangular
## Scaling factor 0.596346493812611
## Bandwidth method user provided
##
## Use summary(...) to show estimates.
##
## $Estplot