1 Procedimiento

  1. Carga, inspección y limpieza de los datos.
  2. Construcción de las variables clave del diseño (exposición, tiempo relativo al tratamiento, periodo pre/post).
  3. Análisis descriptivo y verificación de balance entre grupos.
  4. Estimación del efecto causal mediante DiD escalonado sobre el desenlace principal y los desenlaces secundarios de la cadena de atención.
  5. Verificación de supuestos (tendencias paralelas) y análisis de sensibilidad/robustez.

2 Carga de paquetes

# 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)

3 Carga y limpieza de los datos

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,…

4 Construcción de variables del diseño DiD escalonado

# 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

5 Análisis descriptivo y balance entre grupos

# 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**")
Tabla 1. Características basales por condición de tratamiento municipal
Characteristic Control
N = 6,091
1
Tratado
N = 5,909
1
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)**")
Tabla 2. Desenlaces de la cadena de atención por condición de tratamiento municipal (comparación cruda, sin ajuste)
Characteristic Control
N = 6,091
1
Tratado
N = 5,909
1
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))


6 Estimación del efecto del programa: Diferencias en Diferencias

6.1 Modelo principal: tamización realizada (desenlace principal)

# 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

6.2 Desenlaces secundarios de la cadena de atención

# 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

6.3 Análisis de heterogeneidad del efecto

# 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

7 Verificación de supuestos: tendencias paralelas (event study)

# 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()


8 Análisis de sensibilidad y robustez

8.1 Estimador robusto a heterogeneidad en el momento de adopción (Callaway & Sant’Anna)

# 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)

8.2 Prueba de placebo: fechas de inicio ficticias en municipios control

# 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

8.3 Sensibilidad a la exclusión de municipios extremos en pobreza

# 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