# Se cargan los paquetes necesarios para todo el flujo de trabajo:
library(tidyverse)
library(readxl)
library(lubridate)
library(gtsummary)
library(gt)
library(fixest)
library (did)
library (broom)
library(scales)
ruta_base <- "Base_cancer_cuello_uterino.xlsx"
datos_crudos <- read_excel(ruta_base, skip = 5)
# Inspección inicial: dimensiones y nombres de variables, para confirmar
# que la lectura fue correcta y que se conservaron las 37 columnas esperadas.
glimpse(datos_crudos)
## Rows: 12,000
## Columns: 37
## $ id_registro <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
## $ id_mujer <dbl> 328, 2128, 2763, 4275, 2365, 1542, 2501, …
## $ id_municipio <dbl> 28, 22, 47, 59, 5, 7, 23, 28, 51, 28, 41,…
## $ id_ips <dbl> 107, 1, 72, 130, 135, 44, 149, 175, 92, 6…
## $ fecha_contacto <dttm> 2025-07-18, 2025-12-12, 2023-04-23, 2023…
## $ anio <dbl> 2025, 2025, 2023, 2023, 2025, 2022, 2025,…
## $ mes <dbl> 7, 12, 4, 1, 5, 12, 1, 9, 6, 10, 2, 10, 7…
## $ trimestre <dbl> 3, 4, 2, 1, 2, 4, 1, 3, 2, 4, 1, 4, 3, 3,…
## $ periodo_trimestre <chr> "2025T3", "2025T4", "2023T2", "2023T1", "…
## $ edad <dbl> 51, 35, 29, 23, 56, 42, 30, 51, 48, 24, 2…
## $ aseguramiento <chr> "Subsidiado", "No asegurada", "Contributi…
## $ quintil_SES_1bajo_5alto <dbl> 2, 1, 2, 2, 3, 1, 2, 2, 3, 1, 2, 1, 1, 1,…
## $ region <chr> "Caribe", "Andina", "Andina", "Caribe", "…
## $ ruralidad_muni <chr> "Intermedio", "Urbano", "Urbano", "Urbano…
## $ tamano_poblacional <chr> "<50k", "50-200k", "50-200k", "<50k", "50…
## $ tamizaje_prev_ult_3anios <dbl> 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0,…
## $ tipo_prueba <chr> "Citologia", "HPV", "Citologia", "Citolog…
## $ tamizaje_realizado <dbl> 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1,…
## $ resultado_positivo <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ seguimiento_completo <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ cin2_mas <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ cancer_cuello_invasor <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ estadio_cancer <chr> "NA", "NA", "NA", "NA", "NA", "NA", "NA",…
## $ umbral_edad_programa <dbl> 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 3…
## $ running_edad_centrada <dbl> 21, 5, -1, -7, 26, 12, 0, 21, 18, -6, -1,…
## $ elegible_por_edad <dbl> 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1,…
## $ recibio_componente_intensivo <dbl> 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1,…
## $ umbral_pobreza_priorizacion <dbl> 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 6…
## $ running_pobreza_centrada <dbl> -7.0, 33.0, 8.1, 16.9, 0.6, 8.0, 3.0, -7.…
## $ elegible_por_pobreza <dbl> 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1,…
## $ fecha_inicio_programa_muni <dttm> 2025-01-01, NA, NA, 2024-12-01, 2023-09-…
## $ indice_pobreza_0_100 <dbl> 53.0, 93.0, 68.1, 76.9, 60.6, 68.0, 63.0,…
## $ oferta_servicios_salud_0_1 <dbl> 0.58, 0.62, 0.29, 0.17, 0.47, 0.54, 0.49,…
## $ nivel_ips <chr> "III", "I", "III", "II", "III", "I", "II"…
## $ tipo_ips <chr> "Publica", "Privada", "Mixta", "Mixta", "…
## $ estrato_ips <chr> "Alta_capacidad", "Baja_capacidad", "Medi…
## $ puntaje_riesgo_0_100 <dbl> 75.4, 97.1, 56.5, 75.4, 71.7, 40.9, 66.5,…
# Revisión de valores faltantes por variable.
datos_crudos %>%
summarise(across(everything(), ~ sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "variable", values_to = "n_na") %>%
filter(n_na > 0) %>%
arrange(desc(n_na))
## # A tibble: 1 × 2
## variable n_na
## <chr> <int>
## 1 fecha_inicio_programa_muni 6091
# Se ajustan los tipos de variable
datos <- datos_crudos %>%
mutate(
fecha_contacto = as_date(fecha_contacto),
fecha_inicio_programa_muni = as_date(fecha_inicio_programa_muni),
aseguramiento = factor(aseguramiento),
region = factor(region),
ruralidad_muni = factor(ruralidad_muni,
levels = c("Rural", "Intermedio", "Urbano")),
tamano_poblacional = factor(tamano_poblacional,
levels = c("<50k", "50-200k", ">200k")),
tipo_prueba = factor(tipo_prueba),
nivel_ips = factor(nivel_ips, levels = c("I", "II", "III")),
tipo_ips = factor(tipo_ips),
estrato_ips = factor(estrato_ips,
levels = c("Baja_capacidad", "Media_capacidad",
"Alta_capacidad")),
estadio_cancer = factor(estadio_cancer, levels = c("I", "II", "III", "IV"))
)
# Verificación rápida de la estructura final
glimpse(datos)
## Rows: 12,000
## Columns: 37
## $ id_registro <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
## $ id_mujer <dbl> 328, 2128, 2763, 4275, 2365, 1542, 2501, …
## $ id_municipio <dbl> 28, 22, 47, 59, 5, 7, 23, 28, 51, 28, 41,…
## $ id_ips <dbl> 107, 1, 72, 130, 135, 44, 149, 175, 92, 6…
## $ fecha_contacto <date> 2025-07-18, 2025-12-12, 2023-04-23, 2023…
## $ anio <dbl> 2025, 2025, 2023, 2023, 2025, 2022, 2025,…
## $ mes <dbl> 7, 12, 4, 1, 5, 12, 1, 9, 6, 10, 2, 10, 7…
## $ trimestre <dbl> 3, 4, 2, 1, 2, 4, 1, 3, 2, 4, 1, 4, 3, 3,…
## $ periodo_trimestre <chr> "2025T3", "2025T4", "2023T2", "2023T1", "…
## $ edad <dbl> 51, 35, 29, 23, 56, 42, 30, 51, 48, 24, 2…
## $ aseguramiento <fct> Subsidiado, No asegurada, Contributivo, C…
## $ quintil_SES_1bajo_5alto <dbl> 2, 1, 2, 2, 3, 1, 2, 2, 3, 1, 2, 1, 1, 1,…
## $ region <fct> Caribe, Andina, Andina, Caribe, Caribe, P…
## $ ruralidad_muni <fct> Intermedio, Urbano, Urbano, Urbano, Inter…
## $ tamano_poblacional <fct> <50k, 50-200k, 50-200k, <50k, 50-200k, <5…
## $ tamizaje_prev_ult_3anios <dbl> 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0,…
## $ tipo_prueba <fct> Citologia, HPV, Citologia, Citologia, HPV…
## $ tamizaje_realizado <dbl> 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1,…
## $ resultado_positivo <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ seguimiento_completo <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ cin2_mas <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ cancer_cuello_invasor <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ estadio_cancer <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ umbral_edad_programa <dbl> 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 3…
## $ running_edad_centrada <dbl> 21, 5, -1, -7, 26, 12, 0, 21, 18, -6, -1,…
## $ elegible_por_edad <dbl> 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1,…
## $ recibio_componente_intensivo <dbl> 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1,…
## $ umbral_pobreza_priorizacion <dbl> 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 6…
## $ running_pobreza_centrada <dbl> -7.0, 33.0, 8.1, 16.9, 0.6, 8.0, 3.0, -7.…
## $ elegible_por_pobreza <dbl> 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1,…
## $ fecha_inicio_programa_muni <date> 2025-01-01, NA, NA, 2024-12-01, 2023-09-…
## $ indice_pobreza_0_100 <dbl> 53.0, 93.0, 68.1, 76.9, 60.6, 68.0, 63.0,…
## $ oferta_servicios_salud_0_1 <dbl> 0.58, 0.62, 0.29, 0.17, 0.47, 0.54, 0.49,…
## $ nivel_ips <fct> III, I, III, II, III, I, II, I, II, III, …
## $ tipo_ips <fct> Publica, Privada, Mixta, Mixta, Mixta, Pu…
## $ estrato_ips <fct> Alta_capacidad, Baja_capacidad, Media_cap…
## $ puntaje_riesgo_0_100 <dbl> 75.4, 97.1, 56.5, 75.4, 71.7, 40.9, 66.5,…
# PASO 1: Definir, a nivel municipal, si el municipio implementó el programa
# en algún momento del periodo observado (municipio_tratado = 1) o nunca lo
# implementó (municipio_tratado = 0, sirve como grupo de comparación puro).
datos <- datos %>%
mutate(
municipio_tratado = if_else(!is.na(fecha_inicio_programa_muni), 1, 0)
)
# Conteo de municipios tratados vs. no tratados (grupo de comparación)
datos %>%
distinct(id_municipio, municipio_tratado) %>%
count(municipio_tratado)
## # A tibble: 2 × 2
## municipio_tratado n
## <dbl> <int>
## 1 0 28
## 2 1 27
# PASO 2: Construir el periodo (año-trimestre) de la observación y el periodo
# de inicio del programa en formato comparable
datos <- datos %>%
mutate(
periodo_num = (anio - min(anio)) * 4 + trimestre,
anio_inicio_muni = year(fecha_inicio_programa_muni),
trimestre_inicio_muni = quarter(fecha_inicio_programa_muni),
periodo_inicio_num = if_else(
municipio_tratado == 1,
(anio_inicio_muni - min(anio)) * 4 + trimestre_inicio_muni,
0L # Convención del paquete 'did': 0 = nunca tratado
)
)
# Verificación: tabla de cohortes de entrada al programa (periodo_inicio_num)
datos %>%
filter(municipio_tratado == 1) %>%
distinct(id_municipio, periodo_inicio_num) %>%
count(periodo_inicio_num) %>%
arrange(periodo_inicio_num)
## # A tibble: 9 × 2
## periodo_inicio_num n
## <dbl> <int>
## 1 5 1
## 2 6 4
## 3 7 3
## 4 8 3
## 5 9 3
## 6 10 2
## 7 11 7
## 8 12 3
## 9 13 1
# PASO 3: Construir el tiempo relativo al tratamiento (event time)
datos <- datos %>%
mutate(
tiempo_relativo = if_else(
municipio_tratado == 1,
periodo_num - periodo_inicio_num,
NA_integer_
),
# Indicador simple post-intervención (1 = el contacto ocurrió en el
# municipio después de que el programa fuera activado en ese territorio)
post_programa = if_else(
municipio_tratado == 1 & periodo_num >= periodo_inicio_num, 1, 0
)
)
summary(datos$tiempo_relativo)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -12.000 -5.000 -1.000 -1.336 3.000 11.000 6091
# PASO 4: Definir la exposición individual efectiva, distinguiéndola de la
# exposición territorial.
datos <- datos %>%
mutate(
grupo_exposicion = case_when(
municipio_tratado == 1 & recibio_componente_intensivo == 1 ~ "Expuesta (componente intensivo)",
municipio_tratado == 1 & recibio_componente_intensivo == 0 ~ "Municipio con programa, sin componente",
municipio_tratado == 0 ~ "Municipio sin programa (control)"
)
)
datos %>% count(grupo_exposicion, sort = TRUE)
## # A tibble: 3 × 2
## grupo_exposicion n
## <chr> <int>
## 1 Municipio sin programa (control) 6091
## 2 Municipio con programa, sin componente 3097
## 3 Expuesta (componente intensivo) 2812
# Tabla descriptiva general
datos %>%
mutate(municipio_tratado = factor(municipio_tratado,
labels = c("Control", "Tratado"))) %>%
select(municipio_tratado, edad, aseguramiento, quintil_SES_1bajo_5alto,
ruralidad_muni, indice_pobreza_0_100, oferta_servicios_salud_0_1,
nivel_ips, tipo_ips) %>%
tbl_summary(
by = municipio_tratado,
statistic = list(all_continuous() ~ "{mean} ({sd})"),
missing = "no"
) %>%
add_p() %>%
modify_caption("**Tabla 1. Características basales por condición de tratamiento municipal**")
| Characteristic | Control N = 6,0911 |
Tratado N = 5,9091 |
p-value2 |
|---|---|---|---|
| edad | 37 (9) | 36 (9) | 0.5 |
| aseguramiento | 0.4 | ||
| Contributivo | 2,077 (34%) | 2,087 (35%) | |
| Especial | 496 (8.1%) | 485 (8.2%) | |
| No asegurada | 442 (7.3%) | 392 (6.6%) | |
| Subsidiado | 3,076 (51%) | 2,945 (50%) | |
| quintil_SES_1bajo_5alto | 0.009 | ||
| 1 | 1,845 (30%) | 1,923 (33%) | |
| 2 | 2,122 (35%) | 2,020 (34%) | |
| 3 | 1,671 (27%) | 1,494 (25%) | |
| 4 | 453 (7.4%) | 472 (8.0%) | |
| ruralidad_muni | <0.001 | ||
| Rural | 1,216 (20%) | 1,148 (19%) | |
| Intermedio | 2,211 (36%) | 1,690 (29%) | |
| Urbano | 2,664 (44%) | 3,071 (52%) | |
| indice_pobreza_0_100 | 59 (11) | 60 (13) | 0.7 |
| oferta_servicios_salud_0_1 | 0.55 (0.16) | 0.47 (0.15) | <0.001 |
| nivel_ips | <0.001 | ||
| I | 2,981 (49%) | 3,363 (57%) | |
| II | 1,981 (33%) | 1,972 (33%) | |
| III | 1,129 (19%) | 574 (9.7%) | |
| tipo_ips | <0.001 | ||
| Mixta | 1,101 (18%) | 1,072 (18%) | |
| Privada | 1,148 (19%) | 1,681 (28%) | |
| Publica | 3,842 (63%) | 3,156 (53%) | |
| 1 Mean (SD); n (%) | |||
| 2 Wilcoxon rank sum test; Pearson’s Chi-squared test | |||
# Tabla descriptiva de los desenlaces
datos %>%
mutate(municipio_tratado = factor(municipio_tratado,
labels = c("Control", "Tratado"))) %>%
select(municipio_tratado, tamizaje_realizado, resultado_positivo,
seguimiento_completo, cin2_mas, cancer_cuello_invasor) %>%
tbl_summary(
by = municipio_tratado,
statistic = list(all_dichotomous() ~ "{n} / {N} ({p}%)"),
missing = "no"
) %>%
add_p() %>%
modify_caption("**Tabla 2. Desenlaces de la cadena de atención por condición de tratamiento municipal (comparación cruda, sin ajuste)**")
| Characteristic | Control N = 6,0911 |
Tratado N = 5,9091 |
p-value2 |
|---|---|---|---|
| tamizaje_realizado | 4,215 / 6,091 (69%) | 4,019 / 5,909 (68%) | 0.2 |
| resultado_positivo | 431 / 6,091 (7.1%) | 411 / 5,909 (7.0%) | 0.8 |
| seguimiento_completo | 307 / 6,091 (5.0%) | 289 / 5,909 (4.9%) | 0.7 |
| cin2_mas | 80 / 6,091 (1.3%) | 78 / 5,909 (1.3%) | >0.9 |
| cancer_cuello_invasor | 24 / 6,091 (0.4%) | 22 / 5,909 (0.4%) | 0.8 |
| 1 n / N (%) | |||
| 2 Pearson’s Chi-squared test | |||
# Gráfico de la evolución trimestral de la tasa de tamización realizada,
datos %>%
group_by(periodo_trimestre, municipio_tratado) %>%
summarise(tasa_tamizaje = mean(tamizaje_realizado), .groups = "drop") %>%
mutate(municipio_tratado = factor(municipio_tratado,
labels = c("Control", "Tratado"))) %>%
ggplot(aes(x = periodo_trimestre, y = tasa_tamizaje,
color = municipio_tratado, group = municipio_tratado)) +
geom_line(linewidth = 1) +
geom_point() +
scale_y_continuous(labels = percent_format()) +
labs(
title = "Evolución trimestral de la tasa de tamización realizada",
subtitle = "Comparación entre municipios tratados y municipios control",
x = "Periodo (año-trimestre)", y = "Tasa de tamización",
color = "Condición municipal"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Se estima el modelo de Diferencias en Diferencias mediante regresión con
# efectos fijos por municipio
modelo_tamizaje <- feols(
tamizaje_realizado ~ post_programa | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos
)
summary(modelo_tamizaje)
## OLS estimation, Dep. Var.: tamizaje_realizado
## Observations: 12,000
## Fixed-effects: id_municipio: 55, periodo_trimestre: 16
## Standard-errors: Clustered (id_municipio)
## Estimate Std. Error t value Pr(>|t|)
## post_programa 0.026852 0.014593 1.84008 0.071255 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.462251 Adj. R2: 0.001915
## Within R2: 2.393e-4
# Se reestima el mismo modelo incorporando covariables individuales y
# contextuales
modelo_tamizaje_ajustado <- feols(
tamizaje_realizado ~ post_programa + edad + quintil_SES_1bajo_5alto +
indice_pobreza_0_100 + oferta_servicios_salud_0_1
| id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos
)
summary(modelo_tamizaje_ajustado)
## OLS estimation, Dep. Var.: tamizaje_realizado
## Observations: 12,000
## Fixed-effects: id_municipio: 55, periodo_trimestre: 16
## Standard-errors: Clustered (id_municipio)
## Estimate Std. Error t value Pr(>|t|)
## post_programa 0.026961 0.014510 1.85810 0.0686091 .
## edad 0.000880 0.000452 1.94580 0.0568911 .
## quintil_SES_1bajo_5alto 0.015796 0.005056 3.12426 0.0028659 **
## ... 2 variables were removed because of collinearity
## (indice_pobreza_0_100 and oferta_servicios_salud_0_1)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.462037 Adj. R2: 0.002672
## Within R2: 0.001165
# Tabla comparativa de los dos modelos (sin y con covariables)
etable(modelo_tamizaje, modelo_tamizaje_ajustado,
title = "Efecto del programa sobre la tamización realizada (DiD)")
## modelo_tamizaje modelo_tamizaje_..
## Dependent Var.: tamizaje_realizado tamizaje_realizado
##
## post_programa 0.0268. (0.0146) 0.0270. (0.0145)
## edad 0.0009. (0.0005)
## quintil_SES_1bajo_5alto 0.0158** (0.0051)
## Fixed-Effects: ------------------ ------------------
## id_municipio Yes Yes
## periodo_trimestre Yes Yes
## _______________________ __________________ __________________
## S.E.: Clustered by: id_municipio by: id_municipio
## Observations 12,000 12,000
## R2 0.00774 0.00866
## Within R2 0.00024 0.00117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Se replica la misma especificación del modelo principal para cada uno de
# los desenlaces secundarios definidos en la cadena de atención
modelo_seguimiento <- feols(
seguimiento_completo ~ post_programa | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos %>% filter(resultado_positivo == 1)
)
modelo_cin2 <- feols(
cin2_mas ~ post_programa | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos %>% filter(seguimiento_completo == 1)
)
modelo_cancer <- feols(
cancer_cuello_invasor ~ post_programa | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos
)
etable(modelo_seguimiento, modelo_cin2, modelo_cancer,
headers = c("Seguimiento completo", "CIN2+", "Cáncer invasor"),
title = "Efecto del programa sobre desenlaces secundarios de la cadena de atención")
## modelo_seguimiento modelo_cin2 modelo_cancer
## Seguimiento completo CIN2+ Cáncer invasor
## Dependent Var.: seguimiento_completo cin2_mas cancer_cuello_invasor
##
## post_programa -0.0336 (0.0563) 0.0283 (0.0649) 0.0028 (0.0027)
## Fixed-Effects: -------------------- --------------- ---------------------
## id_municipio Yes Yes Yes
## periodo_trimestre Yes Yes Yes
## _________________ ____________________ _______________ _____________________
## S.E.: Clustered by: id_municipio by: id_munici.. by: id_municipio
## Observations 841 593 12,000
## R2 0.09125 0.09748 0.00504
## Within R2 0.00039 0.00029 0.00015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Se explora si el efecto del programa sobre el desenlace principal varía
# según la ruralidad del municipio, el nivel de complejidad de la IPS y el
# quintil socioeconómico.
modelo_heterog_ruralidad <- feols(
tamizaje_realizado ~ post_programa * ruralidad_muni
| id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos
)
modelo_heterog_nivel_ips <- feols(
tamizaje_realizado ~ post_programa * nivel_ips
| id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos
)
etable(modelo_heterog_ruralidad, modelo_heterog_nivel_ips,
title = "Heterogeneidad del efecto del programa según contexto territorial e institucional")
## modelo_heterog_r.. modelo_heterog_n..
## Dependent Var.: tamizaje_realizado tamizaje_realizado
##
## post_programa 0.0113 (0.0326) 0.0313. (0.0164)
## post_programa x ruralidad_muniIntermedio 0.0330 (0.0327)
## post_programa x ruralidad_muniUrbano 0.0120 (0.0353)
## nivel_ipsII -0.0037 (0.0124)
## nivel_ipsIII 0.0158 (0.0134)
## post_programa x nivel_ipsII -0.0162 (0.0221)
## post_programa x nivel_ipsIII 0.0069 (0.0367)
## Fixed-Effects: ------------------ ------------------
## id_municipio Yes Yes
## periodo_trimestre Yes Yes
## ________________________________________ __________________ __________________
## S.E.: Clustered by: id_municipio by: id_municipio
## Observations 12,000 12,000
## R2 0.00781 0.00796
## Within R2 0.00031 0.00047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Para evaluar visualmente el supuesto de tendencias paralelas, se estima un
# modelo de "event study"
datos_event <- datos %>%
filter(municipio_tratado == 1) %>%
mutate(
tiempo_relativo_cat = factor(tiempo_relativo)
)
datos_event$tiempo_relativo_cat <- relevel(datos_event$tiempo_relativo_cat,
ref = "-1")
# Modelo de event study propiamente dicho, restringido a los municipios
# tratados, comparando cada periodo relativo contra el periodo -1.
modelo_event_study <- feols(
tamizaje_realizado ~ tiempo_relativo_cat | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos_event
)
summary(modelo_event_study)
## OLS estimation, Dep. Var.: tamizaje_realizado
## Observations: 5,909
## Fixed-effects: id_municipio: 27, periodo_trimestre: 16
## Standard-errors: Clustered (id_municipio)
## Estimate Std. Error t value Pr(>|t|)
## tiempo_relativo_cat-12 -0.162256 0.055391 -2.929280 0.006983 **
## tiempo_relativo_cat-11 -0.189373 0.069374 -2.729730 0.011222 *
## tiempo_relativo_cat-10 -0.127443 0.064033 -1.990276 0.057171 .
## tiempo_relativo_cat-9 -0.095456 0.045578 -2.094368 0.046121 *
## tiempo_relativo_cat-8 -0.031979 0.043246 -0.739468 0.466247
## tiempo_relativo_cat-7 -0.027538 0.044947 -0.612683 0.545408
## tiempo_relativo_cat-6 -0.045319 0.042587 -1.064160 0.297040
## tiempo_relativo_cat-5 -0.089051 0.036829 -2.417952 0.022916 *
## tiempo_relativo_cat-4 -0.066140 0.035067 -1.886127 0.070496 .
## tiempo_relativo_cat-3 -0.000419 0.037566 -0.011157 0.991183
## tiempo_relativo_cat-2 -0.012307 0.030385 -0.405033 0.688767
## tiempo_relativo_cat0 -0.015461 0.026793 -0.577055 0.568866
## tiempo_relativo_cat1 -0.028791 0.029564 -0.973855 0.339103
## tiempo_relativo_cat2 0.021813 0.037880 0.575829 0.569682
## tiempo_relativo_cat3 -0.038984 0.027940 -1.395276 0.174737
## tiempo_relativo_cat4 0.018904 0.026510 0.713071 0.482153
## tiempo_relativo_cat5 -0.002595 0.030882 -0.084031 0.933675
## tiempo_relativo_cat6 -0.016268 0.047324 -0.343765 0.733788
## tiempo_relativo_cat7 0.056669 0.050813 1.115235 0.274952
## tiempo_relativo_cat8 0.045850 0.046893 0.977752 0.337207
## tiempo_relativo_cat9 0.012670 0.040413 0.313510 0.756394
## tiempo_relativo_cat10 0.020160 0.073041 0.276015 0.784719
## ... 1 variable was removed because of collinearity
## (tiempo_relativo_cat11)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.46319 Adj. R2: 0.003168
## Within R2: 0.003199
# Se extraen los coeficientes e intervalos de confianza del modelo de event
# study y se grafican en función del tiempo relativo al inicio del programa,
coef_event <- broom::tidy(modelo_event_study, conf.int = TRUE) %>%
mutate(
tiempo_relativo = as.numeric(str_remove(term, "tiempo_relativo_cat"))
) %>%
bind_rows(tibble(tiempo_relativo = -1, estimate = 0,
conf.low = 0, conf.high = 0)) # periodo de referencia
ggplot(coef_event, aes(x = tiempo_relativo, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
geom_vline(xintercept = -0.5, linetype = "dotted", color = "firebrick") +
geom_point(size = 2) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
labs(
title = "Event study: efecto del programa por trimestre relativo al inicio",
subtitle = "Periodo de referencia: t = -1 (un trimestre antes del inicio del programa)",
x = "Trimestres relativos al inicio del programa",
y = "Efecto estimado sobre la tamización realizada"
) +
theme_minimal()
# Estimador de Callaway & Sant'Anna
panel_municipio <- datos %>%
group_by(id_municipio, periodo_num, periodo_inicio_num) %>%
summarise(tasa_tamizaje = mean(tamizaje_realizado), .groups = "drop")
resultado_cs <- att_gt(
yname = "tasa_tamizaje",
tname = "periodo_num",
idname = "id_municipio",
gname = "periodo_inicio_num",
data = panel_municipio,
control_group = "nevertreated",
clustervars = "id_municipio"
)
summary(resultado_cs)
##
## Call:
## att_gt(yname = "tasa_tamizaje", tname = "periodo_num", idname = "id_municipio",
## gname = "periodo_inicio_num", data = panel_municipio, control_group = "nevertreated",
## clustervars = "id_municipio")
##
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
##
## Group-Time Average Treatment Effects:
## Group Time ATT(g,t) Std. Error [95% Simult. Conf. Band]
## 5 2 0.0196 0.0290 -0.0647 0.1039
## 5 3 -0.3986 0.0367 -0.5053 -0.2918 *
## 5 4 0.5312 0.0457 0.3983 0.6642 *
## 5 5 -0.3314 0.0417 -0.4528 -0.2101 *
## 5 6 -0.4734 0.0521 -0.6249 -0.3219 *
## 5 7 -0.0271 0.0429 -0.1518 0.0976
## 5 8 -0.5409 0.0408 -0.6595 -0.4222 *
## 5 9 -0.0124 0.0485 -0.1536 0.1287
## 5 10 -0.0200 0.0550 -0.1801 0.1401
## 5 11 -0.0056 0.0414 -0.1261 0.1148
## 5 12 -0.3367 0.0439 -0.4645 -0.2088 *
## 5 13 0.0260 0.0370 -0.0818 0.1338
## 5 14 0.0483 0.0396 -0.0670 0.1637
## 5 15 -1.0201 0.0454 -1.1522 -0.8881 *
## 5 16 -0.2915 0.0393 -0.4060 -0.1771 *
## 6 2 -0.1110 0.0401 -0.2275 0.0056
## 6 3 0.0810 0.0982 -0.2048 0.3668
## 6 4 -0.0504 0.1449 -0.4722 0.3714
## 6 5 0.0844 0.0758 -0.1361 0.3050
## 6 6 -0.0429 0.0918 -0.3100 0.2243
## 6 7 -0.1703 0.0542 -0.3282 -0.0125 *
## 6 8 0.0113 0.0646 -0.1766 0.1993
## 6 9 -0.0339 0.0756 -0.2540 0.1861
## 6 10 -0.0398 0.0620 -0.2202 0.1405
## 6 11 0.1146 0.0917 -0.1523 0.3815
## 6 12 -0.0606 0.1399 -0.4677 0.3465
## 6 13 -0.0002 0.1060 -0.3087 0.3083
## 6 14 0.2593 0.0641 0.0729 0.4458 *
## 6 15 -0.0913 0.0474 -0.2292 0.0465
## 6 16 -0.0754 0.1641 -0.5530 0.4023
## 7 2 0.0324 0.1400 -0.3750 0.4399
## 7 3 -0.2819 0.0765 -0.5045 -0.0593 *
## 7 4 0.3964 0.1920 -0.1623 0.9551
## 7 5 -0.0688 0.1313 -0.4508 0.3132
## 7 6 -0.0769 0.1678 -0.5651 0.4114
## 7 7 -0.0708 0.1066 -0.3811 0.2394
## 7 8 -0.0821 0.0763 -0.3041 0.1399
## 7 9 -0.0420 0.1467 -0.4690 0.3850
## 7 10 -0.1587 0.0715 -0.3668 0.0493
## 7 11 0.0058 0.2119 -0.6109 0.6226
## 7 12 0.0843 0.1781 -0.4339 0.6025
## 7 13 -0.0994 0.1764 -0.6127 0.4140
## 7 14 0.0503 0.2055 -0.5478 0.6483
## 7 15 -0.1768 0.1306 -0.5568 0.2032
## 7 16 -0.1610 0.1349 -0.5536 0.2315
## 8 2 0.1455 0.1678 -0.3428 0.6337
## 8 3 -0.4338 0.1688 -0.9249 0.0574
## 8 4 0.2786 0.0491 0.1357 0.4215 *
## 8 5 0.0168 0.0828 -0.2242 0.2578
## 8 6 -0.1200 0.1363 -0.5166 0.2765
## 8 7 -0.1141 0.2119 -0.7309 0.5027
## 8 8 0.2902 0.1141 -0.0419 0.6223
## 8 9 0.2853 0.0808 0.0502 0.5204 *
## 8 10 0.2454 0.1645 -0.2335 0.7242
## 8 11 0.3413 0.2083 -0.2648 0.9474
## 8 12 0.1981 0.1313 -0.1840 0.5801
## 8 13 0.3031 0.1304 -0.0762 0.6824
## 8 14 0.3781 0.1403 -0.0303 0.7866
## 8 15 0.0492 0.3135 -0.8630 0.9614
## 8 16 0.1498 0.1459 -0.2748 0.5745
## 9 2 -0.1915 0.1277 -0.5632 0.1802
## 9 3 0.0900 0.0941 -0.1838 0.3638
## 9 4 0.0055 0.0744 -0.2108 0.2219
## 9 5 -0.0108 0.0581 -0.1800 0.1584
## 9 6 -0.0394 0.2377 -0.7312 0.6523
## 9 7 0.0888 0.1238 -0.2714 0.4491
## 9 8 0.0382 0.0411 -0.0813 0.1578
## 9 9 -0.0382 0.0611 -0.2162 0.1397
## 9 10 -0.0120 0.0480 -0.1517 0.1278
## 9 11 -0.2076 0.0589 -0.3790 -0.0362 *
## 9 12 0.0859 0.0511 -0.0629 0.2348
## 9 13 -0.0037 0.0592 -0.1759 0.1686
## 9 14 -0.0505 0.1115 -0.3748 0.2738
## 9 15 -0.0604 0.0833 -0.3029 0.1821
## 9 16 -0.0285 0.0520 -0.1797 0.1228
## 10 2 -0.0162 0.0412 -0.1361 0.1037
## 10 3 0.0490 0.0895 -0.2114 0.3093
## 10 4 -0.2322 0.0746 -0.4494 -0.0150 *
## 10 5 0.0852 0.0417 -0.0361 0.2066
## 10 6 -0.0169 0.0618 -0.1967 0.1630
## 10 7 0.1409 0.0457 0.0080 0.2738 *
## 10 8 0.0710 0.1537 -0.3762 0.5182
## 10 9 0.0275 0.1217 -0.3266 0.3817
## 10 10 -0.2191 0.0679 -0.4166 -0.0216 *
## 10 11 -0.0027 0.0502 -0.1488 0.1434
## 10 12 -0.1219 0.2429 -0.8287 0.5849
## 10 13 -0.2651 0.1141 -0.5971 0.0668
## 10 14 0.0145 0.0461 -0.1198 0.1488
## 10 15 -0.2280 0.0698 -0.4311 -0.0248 *
## 10 16 -0.0897 0.0480 -0.2293 0.0498
## 11 2 0.0508 0.0544 -0.1074 0.2089
## 11 3 -0.0299 0.0560 -0.1928 0.1329
## 11 4 0.0139 0.0680 -0.1841 0.2120
## 11 5 -0.0215 0.0564 -0.1855 0.1425
## 11 6 -0.0852 0.0599 -0.2595 0.0890
## 11 7 0.1320 0.0628 -0.0509 0.3148
## 11 8 -0.0618 0.0803 -0.2956 0.1720
## 11 9 0.0601 0.1121 -0.2660 0.3862
## 11 10 -0.0171 0.0829 -0.2583 0.2240
## 11 11 0.0225 0.0695 -0.1798 0.2249
## 11 12 0.0543 0.0971 -0.2283 0.3369
## 11 13 0.1430 0.0778 -0.0833 0.3693
## 11 14 0.0970 0.0644 -0.0904 0.2843
## 11 15 0.0155 0.0488 -0.1265 0.1575
## 11 16 -0.0465 0.0817 -0.2842 0.1912
## 12 2 -0.0862 0.1673 -0.5731 0.4008
## 12 3 -0.0323 0.0865 -0.2841 0.2196
## 12 4 0.0993 0.0876 -0.1556 0.3542
## 12 5 0.0321 0.0444 -0.0969 0.1612
## 12 6 0.0027 0.0981 -0.2829 0.2883
## 12 7 0.0187 0.1133 -0.3109 0.3484
## 12 8 -0.1821 0.0706 -0.3876 0.0234
## 12 9 0.1708 0.0521 0.0191 0.3225 *
## 12 10 -0.0710 0.1601 -0.5368 0.3948
## 12 11 0.0749 0.0768 -0.1486 0.2985
## 12 12 0.0932 0.0743 -0.1229 0.3093
## 12 13 0.0316 0.0733 -0.1816 0.2448
## 12 14 0.0387 0.1237 -0.3214 0.3987
## 12 15 -0.0783 0.1141 -0.4102 0.2536
## 12 16 0.1317 0.1186 -0.2135 0.4769
## 13 2 -0.1113 0.0290 -0.1956 -0.0270 *
## 13 3 0.1504 0.0367 0.0436 0.2571 *
## 13 4 -0.1132 0.0457 -0.2462 0.0197
## 13 5 0.1034 0.0417 -0.0180 0.2247
## 13 6 -0.0050 0.0482 -0.1453 0.1352
## 13 7 0.0017 0.0462 -0.1327 0.1361
## 13 8 -0.1461 0.0358 -0.2502 -0.0420 *
## 13 9 0.1351 0.0488 -0.0069 0.2771
## 13 10 0.1007 0.0491 -0.0422 0.2437
## 13 11 -0.1106 0.0528 -0.2642 0.0430
## 13 12 0.0287 0.0503 -0.1176 0.1750
## 13 13 0.0452 0.0475 -0.0932 0.1836
## 13 14 0.0274 0.0506 -0.1198 0.1746
## 13 15 0.0999 0.0455 -0.0325 0.2323
## 13 16 -0.0289 0.0434 -0.1552 0.0974
## ---
## Signif. codes: `*' confidence band does not cover 0
##
## Control Group: Never Treated, Anticipation Periods: 0
## Estimation Method: Doubly Robust
# Agregación del efecto a un único parámetro resumen
agregado_simple <- aggte(resultado_cs, type = "simple")
agregado_dinamico <- aggte(resultado_cs, type = "dynamic")
summary(agregado_simple)
##
## Call:
## aggte(MP = resultado_cs, type = "simple")
##
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
##
##
## ATT Std. Error [ 95% Conf. Int.]
## 0.0061 0.0374 -0.0672 0.0794
##
##
## ---
## Signif. codes: `*' confidence band does not cover 0
##
## Control Group: Never Treated, Anticipation Periods: 0
## Estimation Method: Doubly Robust
summary(agregado_dinamico)
##
## Call:
## aggte(MP = resultado_cs, type = "dynamic")
##
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna. "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015>
##
##
## Overall summary of ATT's based on event-study/dynamic aggregation:
## ATT Std. Error [ 95% Conf. Int.]
## -0.0382 0.0386 -0.1139 0.0375
##
##
## Dynamic Effects:
## Event time Estimate Std. Error [95% Simult. Conf. Band]
## -11 -0.1113 0.0296 -0.1919 -0.0308 *
## -10 -0.0270 0.1413 -0.4122 0.3581
## -9 0.0132 0.0443 -0.1075 0.1339
## -8 0.0122 0.0383 -0.0922 0.1167
## -7 -0.0180 0.0455 -0.1421 0.1062
## -6 0.0053 0.0433 -0.1128 0.1235
## -5 -0.0774 0.0506 -0.2154 0.0606
## -4 -0.0003 0.0388 -0.1059 0.1054
## -3 0.0714 0.0466 -0.0556 0.1984
## -2 -0.0245 0.0457 -0.1490 0.1000
## -1 0.0222 0.0439 -0.0974 0.1418
## 0 0.0031 0.0391 -0.1034 0.1097
## 1 -0.0031 0.0410 -0.1148 0.1086
## 2 0.0362 0.0451 -0.0867 0.1592
## 3 0.0005 0.0502 -0.1364 0.1374
## 4 0.0370 0.0360 -0.0612 0.1351
## 5 0.0290 0.0514 -0.1111 0.1692
## 6 0.0142 0.0726 -0.1837 0.2122
## 7 -0.0089 0.0862 -0.2438 0.2261
## 8 0.0893 0.0875 -0.1493 0.3279
## 9 -0.1000 0.0563 -0.2535 0.0535
## 10 -0.2643 0.2157 -0.8525 0.3238
## 11 -0.2915 0.0379 -0.3947 -0.1883 *
## ---
## Signif. codes: `*' confidence band does not cover 0
##
## Control Group: Never Treated, Anticipation Periods: 0
## Estimation Method: Doubly Robust
ggdid(agregado_dinamico)
# Como prueba de robustez adicional, se asignan fechas de inicio ficticias
# (aleatorias) a los municipios control
set.seed(123)
fechas_placebo <- datos %>%
filter(municipio_tratado == 0) %>%
distinct(id_municipio) %>%
mutate(
periodo_inicio_placebo = sample(unique(datos$periodo_num),
size = n(), replace = TRUE)
)
datos_placebo <- datos %>%
filter(municipio_tratado == 0) %>%
left_join(fechas_placebo, by = "id_municipio") %>%
mutate(post_placebo = if_else(periodo_num >= periodo_inicio_placebo, 1, 0))
modelo_placebo <- feols(
tamizaje_realizado ~ post_placebo | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos_placebo
)
summary(modelo_placebo)
## OLS estimation, Dep. Var.: tamizaje_realizado
## Observations: 6,091
## Fixed-effects: id_municipio: 28, periodo_trimestre: 16
## Standard-errors: Clustered (id_municipio)
## Estimate Std. Error t value Pr(>|t|)
## post_placebo -0.00879 0.022293 -0.394294 0.69646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.459762 Adj. R2: 0.001171
## Within R2: 2.949e-5
# Se reestima el modelo principal excluyendo los municipios extremos según con
# condición socioeconómica
limites_pobreza <- quantile(datos$indice_pobreza_0_100, probs = c(0.05, 0.95))
datos_sin_extremos <- datos %>%
filter(indice_pobreza_0_100 >= limites_pobreza[1],
indice_pobreza_0_100 <= limites_pobreza[2])
modelo_sensibilidad <- feols(
tamizaje_realizado ~ post_programa | id_municipio + periodo_trimestre,
cluster = ~id_municipio,
data = datos_sin_extremos
)
etable(modelo_tamizaje, modelo_sensibilidad,
headers = c("Modelo principal", "Excluyendo municipios extremos"),
title = "Sensibilidad del efecto principal a la exclusión de municipios con pobreza extrema")
## modelo_tamizaje modelo_sensibili..
## Modelo principal Excluyendo municipios extremos
## Dependent Var.: tamizaje_realizado tamizaje_realizado
##
## post_programa 0.0268. (0.0146) 0.0188 (0.0155)
## Fixed-Effects: ------------------ ------------------
## id_municipio Yes Yes
## periodo_trimestre Yes Yes
## _________________ __________________ __________________
## S.E.: Clustered by: id_municipio by: id_municipio
## Observations 12,000 11,045
## R2 0.00774 0.00794
## Within R2 0.00024 0.00012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1