Respuestas a la tarea 3

Pregunta 1

Hace unas semanas se presentaron los resultados de una evaluación del impacto del programa Jóvenes Construyendo el Futuro (JCF), realizada usando métodos de matching. Las tablas 1 y 2 del reporte muestran el ATE estimado en el ingreso trimestral entre los jóvenes que no asisten a la escuela y no están empleados y el ATE en la probabilidad de encontrar un trabajo entre los jóvenes en general, respectivamente. En este ejercicio extenderemos los resultados encontrados.

Los datos en datos_jcf_analisis.csv están listos para analizarse. Estos se construyeron a partir de la ENIGH 2020, que incluyó un módulo especial para el programa JCF y que pueden descargar de aquí. Para limpiar los datos use el script en limpiar_jcf_analisis.r. Para ejecutar este script debería descargar de la página de INEGI antes citada los archivos ingresos_jcf.csv, ingresos.csv, poblacion.csv, viviendas.csv, hogares.csv y concentradohogar.csv. Esto sería necesario si quisiera agregar nuevas variables al análisis, pero bien puede trabajar con los datos en datos_jcf_analisis.csv que yo están limpios.

El propensity score (PS) usado en la evaluación usa los siguientes regresores: mujer (dummy de sexo), indigena (dummy de pertenencia a una etnia), rural (dummy del ámbito rural), escoacum (años de escolaridad), casadounion (dummy para casados o en unión libre), jefehog (dummy para jefes del hogar), haymenores (dummy para la presencia de menores de edad en el hogar), proggob (dummy para beneficiarios de programas de gobierno), y tot_integ (número de miembros del hogar), así como dummies de estado, cve_ent.

  1. [10 puntos] Estime ahora el TOT (TT o ATT) del programa en el ingreso trimestral, ingtot_tri. Para estimar el impacto en el ingreso trimestral se comparan a los beneficiarios de JCF con los jóvenes que no asisten a la escuela y no están empleados. Los beneficiarios tienen jcf2==1 y los jóvenes que no asisten a la escuela y no están empleados tienen jcf2==0. Realice la inferencia estadística usando el método de simulación de King, Tomz y Wittenberg (2000) visto en clase. ¿De qué tamaño es el TOT estimado y es este efecto estadísticamente significativo?

    Limpiamos datos originales de la ENIGH:

    
    setwd("C:/Users/rojas/Dropbox/Evaluación de Programas Sociales/2021/tareas/Tarea 3")
    
    #Ingresos de JCF
    ing.jcf <- read_csv(
      "./ingresos_jcf.csv",
      locale = locale(encoding = "latin1")) %>% 
      rename(ingjcf_1=ing_6, #Arreglar nombres para que los ingresos estén en orden cronológico
         ingjcf_2=ing_5,
         ingjcf_3=ing_4,
         ingjcf_4=ing_3,
         ingjcf_5=ing_2,
         ingjcf_6=ing_1,
         ingjcf_tri=ing_tri,
         mesjcf_1=mes_6,
         mesjcf_2=mes_5,
         mesjcf_3=mes_4,
         mesjcf_4=mes_3,
         mesjcf_5=mes_2,
         mesjcf_6=mes_1) %>% 
      dplyr::select(-c(clave))
    
    #Ingreso laboral
    ing.lab <- read_csv(
      "./ingresos.csv",
      locale = locale(encoding = "latin1")) %>% 
      filter(clave %in% c("P001","P002","P003","P004",
                      "P005","P006","P007","P010",
                      "P011","P012","P013","P014",
                      "P018","P019","P020","P021","P022")) 
    
    
    ing.lab.fuentes <- ing.lab %>%
      dplyr::select(folioviv,foliohog,numren,ing_tri,
         ing_1=ing_6, #Arreglar nombres para que los ingresos estén en orden cronológico
           ing_2=ing_5,
           ing_3=ing_4,
           ing_4=ing_3,
           ing_5=ing_2,
           ing_6=ing_1) %>% 
      group_by(folioviv,foliohog,numren) %>%
      fsum
    
    ing.lab.meses <- ing.lab %>% 
      dplyr::select(folioviv,foliohog,numren,
         mes_1=mes_6,
         mes_2=mes_5,
         mes_3=mes_4,
         mes_4=mes_3,
         mes_5=mes_2,
         mes_6=mes_1) %>% 
      mutate(mes_1=as.numeric(mes_1),
         mes_2=as.numeric(mes_2),
         mes_3=as.numeric(mes_3),
         mes_4=as.numeric(mes_4),
         mes_5=as.numeric(mes_5),
         mes_6=as.numeric(mes_6)) %>% 
      group_by(folioviv,foliohog,numren) %>%
      fmean
    
    ing.todos <- ing.lab.fuentes %>% 
      left_join(ing.lab.meses, by=c("folioviv","foliohog","numren")) %>% 
      full_join(ing.jcf, by=c("folioviv","foliohog","numren")) %>% 
      mutate(ing_1=ifelse(is.na(ing_1),0,ing_1),
         ing_2=ifelse(is.na(ing_2),0,ing_2),
         ing_3=ifelse(is.na(ing_3),0,ing_3),
         ing_4=ifelse(is.na(ing_4),0,ing_4),
         ing_5=ifelse(is.na(ing_5),0,ing_5),
         ing_6=ifelse(is.na(ing_6),0,ing_6),
         ingjcf_1=ifelse(is.na(ingjcf_1),0,ingjcf_1),
         ingjcf_2=ifelse(is.na(ingjcf_2),0,ingjcf_2),
         ingjcf_3=ifelse(is.na(ingjcf_3),0,ingjcf_3),
         ingjcf_4=ifelse(is.na(ingjcf_4),0,ingjcf_4),
         ingjcf_5=ifelse(is.na(ingjcf_5),0,ingjcf_5),
         ingjcf_6=ifelse(is.na(ingjcf_6),0,ingjcf_6),
         ct_futuro=ifelse(is.na(ct_futuro),0,ct_futuro),
         ing_1=ifelse(ct_futuro!=9,ing_1-ingjcf_1,ing_1),
         ing_2=ifelse(ct_futuro!=9,ing_2-ingjcf_2,ing_2),
         ing_3=ifelse(ct_futuro!=9,ing_3-ingjcf_3,ing_3),
         ing_4=ifelse(ct_futuro!=9,ing_4-ingjcf_4,ing_4),
         ing_5=ifelse(ct_futuro!=9,ing_5-ingjcf_5,ing_5),
         ing_6=ifelse(ct_futuro!=9,ing_6-ingjcf_6,ing_6),
         ing_tri=ifelse(is.na(ing_tri),0,ing_tri),
         ingjcf_tri=ifelse(is.na(ingjcf_tri),0,ingjcf_tri),
         ing_tri=ifelse(ct_futuro!="9",ing_tri-ingjcf_tri,ing_tri),
         ingtot_tri=ingjcf_tri+ing_tri)
    
    
    poblacion <- read_csv(
      "./poblacion.csv",
      locale = locale(encoding = "latin1"),
      col_types = list(min_8 = col_double())) %>% 
      mutate(mujer=ifelse(sexo==2,1,0),
         indigena=ifelse(etnia==1,1,0),
         analfabeta=ifelse(alfabetism==2,1,0),
         joven=ifelse(edad>=18 & edad <=29,1,0),
         jnc=ifelse(edad>=18 & edad <=29 & asis_esc==2 & trabajo_mp==2,1,0),
         jtrabaja=ifelse(edad>=18 & edad <=29 & trabajo_mp==1,1,0),
         escoacum=case_when(nivelaprob==0 | nivelaprob==1 ~ 0,
                            nivelaprob==2 ~ gradoaprob + 1,
                            nivelaprob==3 ~ gradoaprob + 6,
                            nivelaprob==4 ~ gradoaprob + 9,
                            nivelaprob==5 & antec_esc==1 ~ gradoaprob + 6,
                            nivelaprob==5 & antec_esc==2 ~ gradoaprob + 9,
                            nivelaprob==5 & antec_esc==3 ~ gradoaprob + 12,
                            nivelaprob==6 & antec_esc==1 ~ gradoaprob + 6,
                            nivelaprob==6 & antec_esc==2 ~ gradoaprob + 9,
                            nivelaprob==6 & antec_esc==3 ~ gradoaprob + 12,
                            nivelaprob==7 ~ gradoaprob + 12,
                            nivelaprob==8 ~ gradoaprob + 16,
                            nivelaprob==9 ~ gradoaprob + 18),
         casadounion=ifelse(edo_conyug==1 | edo_conyug==2,1,0),
         mintrab=ifelse(is.na(hor_1) | is.na(min_1),0,(hor_1*60)+min_1),
         minest=ifelse(is.na(hor_2) | is.na(min_2),0,(hor_2*60)+min_2),
         mincuid=ifelse(is.na(hor_4) | is.na(min_4),0,(hor_4*60)+min_4),
         minqueh=ifelse(is.na(hor_6) | is.na(min_6),0,(hor_6*60)+min_6),
         minrecr=ifelse(is.na(hor_8) | is.na(min_8),0,(as.numeric(hor_8)*60)+min_8),
         afilss=ifelse(pop_insabi==1 | atemed==1, 1, 0),
         jefehog=ifelse(parentesco==101,1,0))
    
    
    viviendas <- read_csv("./viviendas.csv",
                      col_types = list(mat_pisos = col_character()),
                      locale = locale(encoding = "latin1")) %>% 
      clean_names() %>% 
      mutate(pisotierra=ifelse(mat_pisos==1, 1 , 0),
           aguaent=ifelse(disp_agua==1 | disp_agua==2,1,0),
           drenfosa=ifelse(drenaje==1 | drenaje==2,1,0),
           sinelec=ifelse(disp_elect==5,1,0))
    
    hogares <- read_csv(
      "./hogares.csv",
      locale = locale(encoding = "latin1")) %>% 
      mutate(sincomida=ifelse(acc_alim2==1,1,0),
          internet=ifelse(conex_inte==1,1,0))
    
    concentrado <- read_csv(
      "./concentradohogar.csv",
      locale = locale(encoding = "latin1")) %>% 
      mutate(cve_ent=substr(ubica_geo,1,2),
         cve_mun=substr(ubica_geo,3,6),
         rural=ifelse(tam_loc==4,1,0),
         haymenores=ifelse(menores!=0,1,0))
    
    #Un solo conjunto de datos
    
    data.jcf <- poblacion %>% 
      dplyr::select(-ct_futuro) %>% 
      left_join(viviendas, by=c("folioviv")) %>% 
      left_join(hogares, by=c("folioviv", "foliohog")) %>% 
      left_join(concentrado, by=c("folioviv", "foliohog")) %>% 
      left_join(ing.todos, by=c("folioviv","foliohog","numren")) %>% 
      mutate(ing_1=ifelse(is.na(ing_1),0,ing_1),
         ing_2=ifelse(is.na(ing_2),0,ing_2),
         ing_3=ifelse(is.na(ing_3),0,ing_3),
         ing_4=ifelse(is.na(ing_4),0,ing_4),
         ing_5=ifelse(is.na(ing_5),0,ing_5),
         ing_6=ifelse(is.na(ing_6),0,ing_6),
         ingjcf_1=ifelse(is.na(ingjcf_1),0,ingjcf_1),
         ingjcf_2=ifelse(is.na(ingjcf_2),0,ingjcf_2),
         ingjcf_3=ifelse(is.na(ingjcf_3),0,ingjcf_3),
         ingjcf_4=ifelse(is.na(ingjcf_4),0,ingjcf_4),
         ingjcf_5=ifelse(is.na(ingjcf_5),0,ingjcf_5),
         ingjcf_6=ifelse(is.na(ingjcf_6),0,ingjcf_6),
         ingtot_tri=ifelse(is.na(ingtot_tri),0,ingtot_tri),
         ing_tri=ifelse(is.na(ing_tri),0,ing_tri),
         ingjcf_tri=ifelse(is.na(ingjcf_tri),0,ingjcf_tri))
    
      #Definición de grupos de comparación
    data.jcf <- data.jcf %>% 
      mutate(jcf=ifelse(ct_futuro %in% c(1,2,3,8,9),1,0)) %>% 
      filter((edad>=18 & edad <=29) | jcf==1) %>% 
      mutate(ingbene=ifelse(ct_futuro==9,bene_gob-ingjcf_tri,NA),
        ingbene=ifelse(jcf==0,bene_gob,ingbene),
        ingbene=ifelse(ingbene<0 | is.na(ingbene),0,ingbene),
        proggob=ifelse(ingbene>0,1,0)) %>% 
      mutate(jcf2=ifelse(jnc==1,0,NA),
         jcf2=ifelse(jcf==1,1,jcf2)) %>% 
      mutate(jcf3=ifelse(jtrabaja==1 & is.na(jcf2),0,NA),
         jcf3=ifelse(jcf==1,1,jcf3)) %>% 
      mutate(trabajo1=ifelse(ing_1>0 & !is.na(ing_1),1,0),
         trabajo2=ifelse(ing_2>0 & !is.na(ing_2),1,0),
         trabajo3=ifelse(ing_3>0 & !is.na(ing_3),1,0),
         trabajo4=ifelse(ing_4>0 & !is.na(ing_4),1,0),
         trabajo5=ifelse(ing_5>0 & !is.na(ing_5),1,0),
         trabajo6=ifelse(ing_6>0 & !is.na(ing_6),1,0)) %>% 
      mutate(benef1=ifelse(ingjcf_1>0 & !is.na(ingjcf_1),1,0),
         benef2=ifelse(ingjcf_2>0 & !is.na(ingjcf_2),1,0),
         benef3=ifelse(ingjcf_3>0 & !is.na(ingjcf_3),1,0),
         benef4=ifelse(ingjcf_4>0 & !is.na(ingjcf_4),1,0),
         benef5=ifelse(ingjcf_5>0 & !is.na(ingjcf_5),1,0),
         benef6=ifelse(ingjcf_6>0 & !is.na(ingjcf_6),1,0))
    
    #Transiciones de empleo
    
    data.jcf <- data.jcf %>%
      mutate(jcf_a_emp=0,
         jcf_a_emp=ifelse(trabajo2==1 & benef1==1,1,jcf_a_emp),
         jcf_a_emp=ifelse(trabajo3==1 & (benef1==1 | benef2==1),1,jcf_a_emp),
         jcf_a_emp=ifelse(trabajo4==1 & (benef1==1 | benef2==1 | benef3==1),1,jcf_a_emp),
         jcf_a_emp=ifelse(trabajo5==1 & (benef1==1 | benef2==1 | benef3==1 | benef4==1),1,jcf_a_emp),
         jcf_a_emp=ifelse(trabajo6==1 & (benef1==1 | benef2==1 | benef3==1 | benef4==1 | benef5==1),1,jcf_a_emp)) %>% 
      mutate(desemp_a_emp=0,
         desemp_a_emp=ifelse(trabajo1==0 & (trabajo2==1 | trabajo3==1 | trabajo4==1 | trabajo5==1 | trabajo6==1) & benef1==0 & jcf==0,1,desemp_a_emp),
         desemp_a_emp=ifelse(trabajo2==0 & (trabajo3==1 | trabajo4==1 | trabajo5==1 | trabajo6==1) & benef2==0 & jcf==0,1,desemp_a_emp),
         desemp_a_emp=ifelse(trabajo3==0 & (trabajo4==1 | trabajo5==1 | trabajo6==1) & benef3==0 & jcf==0,1,desemp_a_emp),
         desemp_a_emp=ifelse(trabajo4==0 & (trabajo5==1 | trabajo6==1) & benef4==0 & jcf==0,1,desemp_a_emp),
         desemp_a_emp=ifelse(trabajo5==0 & (trabajo6==1) & benef5==0 & jcf==0,1,desemp_a_emp)) %>% 
      mutate(siempre_jcf=ifelse(benef1==1 & benef2==1 & benef3==1 & benef4==1 & benef5==1 & benef6==1,1,0),
         siempre_emp=ifelse((trabajo1==1 & trabajo2==1 & trabajo3==1 & trabajo4==1 & trabajo5==1 & trabajo6==1) & jcf==0,1,0),
         encontro=ifelse(jcf_a_emp==1 | desemp_a_emp==1,1,0)) %>% 
      mutate(transicion=ifelse(benef6==1 | siempre_jcf==1| siempre_emp==1,0,1))

    Estadística descriptiva:

    descr(dplyr::select(data.jcf, jcf, jcf2, jcf3, ingtot_tri,
                    mujer, indigena, cve_ent, rural,
                    escoacum, casadounion, jefehog, haymenores,
                    proggob, tot_integ, factor.x),
      round.digits = 2,
      headings = FALSE, # quitar encabezados
      stats = "common") %>%  # estadísticas más usadas
      tb()
    ## Error in pryr::where(obj_name) : length(name) == 1 is not TRUE
    ## # A tibble: 14 x 8
    ##    variable         mean        sd   min   med     max n.valid pct.valid
    ##    <chr>           <dbl>     <dbl> <dbl> <dbl>   <dbl>   <dbl>     <dbl>
    ##  1 casadounion    0.374      0.484     0    0       1    60660     100  
    ##  2 escoacum      11.3        3.13      0   12      23    60659     100. 
    ##  3 factor.x     407.       411.       10  282    5284    60660     100  
    ##  4 haymenores     0.542      0.498     0    1       1    60660     100  
    ##  5 indigena       0.301      0.459     0    0       1    60660     100  
    ##  6 ingtot_tri  9486.     14550.        0 3688. 727411.   60660     100  
    ##  7 jcf            0.0142     0.118     0    0       1    60660     100  
    ##  8 jcf2           0.0567     0.231     0    0       1    15152      25.0
    ##  9 jcf3           0.0223     0.148     0    0       1    38446      63.4
    ## 10 jefehog        0.125      0.331     0    0       1    60660     100  
    ## 11 mujer          0.504      0.500     0    1       1    60660     100  
    ## 12 proggob        0.269      0.444     0    0       1    60660     100  
    ## 13 rural          0.372      0.483     0    0       1    60660     100  
    ## 14 tot_integ      4.72       2.02      1    4      25    60660     100

    Uso matchit para construir las muestras emparejadas:

    sub.data <- data.jcf %>%
        dplyr::select(ingtot_tri, jcf2, mujer, indigena, cve_ent, rural, escoacum, casadounion,
            jefehog, haymenores, proggob, tot_integ, factor.x)
    
    sub.data <- sub.data[complete.cases(sub.data), ]
    
    m.out.a <- matchit(formula = jcf2 ~ mujer + indigena + factor(cve_ent) + rural +
        escoacum + casadounion + jefehog + haymenores + proggob + tot_integ, method = "nearest",
        distance = "glm", replace = FALSE, data = sub.data)

    La simulación de King, Tomz y Wittenberg (2000):

    set.seed(1711)
    
    z.out.a <- zelig(formula = ingtot_tri ~ jcf2, data = match.data(m.out.a), model = "ls")
    ## How to cite this model in Zelig:
    ##   R Core Team. 2007.
    ##   ls: Least Squares Regression for Continuous Dependent Variables
    ##   in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
    ##   "Zelig: Everyone's Statistical Software," https://zeligproject.org/
    
    # Simulo las dos situaciones
    x.out.a <- setx(z.out.a, jcf2 = 0)
    x1.out.a <- setx1(z.out.a, jcf2 = 1)
    
    # Corremos la simulación
    sim.out.a <- sim(z.out.a, x = x.out.a, x1 = x1.out.a)
    
    # Vemos los resultados
    summary(sim.out.a)
    ## 
    ##  sim x :
    ##  -----
    ## ev
    ##       mean       sd      50%     2.5%    97.5%
    ## 1 1628.963 235.5454 1630.271 1156.995 2064.117
    ## pv
    ##          mean       sd     50%      2.5%    97.5%
    ## [1,] 1369.212 7060.959 1437.84 -12210.35 14859.28
    ## 
    ##  sim x1 :
    ##  -----
    ## ev
    ##      mean       sd      50%     2.5%    97.5%
    ## 1 10102.6 240.5413 10103.57 9626.654 10564.29
    ## pv
    ##          mean       sd      50%      2.5%    97.5%
    ## [1,] 9847.497 7157.148 9936.032 -3541.964 23555.74
    ## fd
    ##       mean      sd      50%     2.5%    97.5%
    ## 1 8473.637 326.528 8467.072 7850.093 9078.618

    El TOT estimado es 8473.637 pesos adicionales (e.e. 326.528) para los beneficiarios de JCF, comparados con otros jóvenes que no asisten a la escuela ni están empleados. Este efecto es muy parecido al ATE estimado en el reporte, aunque el error estándar es considerablemente más pequeño.

  2. [5 puntos] En el matching de la parte a., evalúe qué tan bueno es el procedimiento en balancear las características observadas una vez realizado el matching. Cree un love plot y realice pruebas formales para contrastar las diferencias en características observables antes y después del matching.

    Veamos un resumen del emparejamiento:

    # Con esto elimino las dummies de estado de la salida
    m.out.a[["X"]][["factor(cve_ent)"]] <- NULL
    
    summary(m.out.a, standardize = T)
    ## 
    ## Call:
    ## matchit(formula = jcf2 ~ mujer + indigena + factor(cve_ent) + 
    ##     rural + escoacum + casadounion + jefehog + haymenores + proggob + 
    ##     tot_integ, data = sub.data, method = "nearest", distance = "glm", 
    ##     replace = FALSE)
    ## 
    ## Summary of Balance for All Data:
    ##             Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
    ## distance           0.1687        0.0500          0.8203     4.3918    0.3349   0.4973
    ## mujer              0.5867        0.7810         -0.3945          .    0.1943   0.1943
    ## indigena           0.4389        0.2997          0.2806          .    0.1392   0.1392
    ## rural              0.4948        0.4161          0.1572          .    0.0786   0.0786
    ## escoacum          11.7590       10.4076          0.4698     0.7742    0.0606   0.1965
    ## casadounion        0.3946        0.5367         -0.2906          .    0.1421   0.1421
    ## jefehog            0.1048        0.0486          0.1833          .    0.0561   0.0561
    ## haymenores         0.5786        0.6717         -0.1886          .    0.0931   0.0931
    ## proggob            0.3446        0.2473          0.2048          .    0.0973   0.0973
    ## tot_integ          4.7544        4.9276         -0.0853     0.9603    0.0092   0.0394
    ## 
    ## 
    ## Summary of Balance for Matched Data:
    ##             Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
    ## distance           0.1687        0.1661          0.0179     1.1056    0.0002   0.0198
    ## mujer              0.5867        0.6298         -0.0875          .    0.0431   0.0431
    ## indigena           0.4389        0.4412         -0.0047          .    0.0023   0.0023
    ## rural              0.4948        0.4843          0.0210          .    0.0105   0.0105
    ## escoacum          11.7590       11.9744         -0.0749     0.8116    0.0137   0.0861
    ## casadounion        0.3946        0.3679          0.0548          .    0.0268   0.0268
    ## jefehog            0.1048        0.0920          0.0418          .    0.0128   0.0128
    ## haymenores         0.5786        0.5774          0.0024          .    0.0012   0.0012
    ## proggob            0.3446        0.3411          0.0073          .    0.0035   0.0035
    ## tot_integ          4.7544        4.7043          0.0246     1.1304    0.0057   0.0291
    ##             Std. Pair Dist.
    ## distance             0.0184
    ## mujer                0.6549
    ## indigena             0.7554
    ## rural                0.8126
    ## escoacum             0.8843
    ## casadounion          0.6883
    ## jefehog              0.5132
    ## haymenores           0.7521
    ## proggob              0.7226
    ## tot_integ            0.8294
    ## 
    ## Percent Balance Improvement:
    ##             Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
    ## distance               97.8       93.2      99.9     96.0
    ## mujer                  77.8          .      77.8     77.8
    ## indigena               98.3          .      98.3     98.3
    ## rural                  86.7          .      86.7     86.7
    ## escoacum               84.1       18.4      77.4     56.2
    ## casadounion            81.2          .      81.2     81.2
    ## jefehog                77.2          .      77.2     77.2
    ## haymenores             98.8          .      98.8     98.8
    ## proggob                96.4          .      96.4     96.4
    ## tot_integ              71.1     -202.9      38.0     26.1
    ## 
    ## Sample Sizes:
    ##           Control Treated
    ## All         14293     859
    ## Matched       859     859
    ## Unmatched   13434       0
    ## Discarded       0       0

    No hay un solo criterio para decidir cuándo el emparejamiento fue exitoso. Una propuesta es observar las diferencias promedio estandarizadas (SMD). Esto nos permite eliminar el tema de la escala y concentrarnos en las diferencias entre los grupos tratados y no tratados, antes y después del emparejamiento. Podemos definitir la SMD como:

    \[SMD_X=\frac{\bar{X}_T-\bar{X}_{NT}}{\sqrt{S^2_T+S^2_{NT}}}\]

    También vale la pena no perder de vista la razón de varianzas (VR). Se espera que este ratio no sea muy distinto de 1 después de hacer el emparejamiento:

    \[VR=\frac{S^2_T}{S^2_{NT}}\]

    Una regla de dedo de 0.1 en el SMD para juzgar el balance. Siguiendo esta regla, todas las variables están balanceadas después del emparejamiento. Por ejemplo, la escolaridad acumulada tenía un SDM de 0.4698 en la muestra en bruto, pero con el emparejamiento el SDM se vuelve de solo 0.0749.

    La librería cobalt tiene una función que hace algo similar, de una forma más compacta:

    # Con esto elimino las dummies de estado de la salida
    bal.tab(m.out.a, m.threshold = 0.1, un = T)
    ## Call
    ##  matchit(formula = jcf2 ~ mujer + indigena + factor(cve_ent) + 
    ##     rural + escoacum + casadounion + jefehog + haymenores + proggob + 
    ##     tot_integ, data = sub.data, method = "nearest", distance = "glm", 
    ##     replace = FALSE)
    ## 
    ## Balance Measures
    ##                 Type Diff.Un Diff.Adj    M.Threshold
    ## distance    Distance  0.8203   0.0179 Balanced, <0.1
    ## mujer         Binary -0.1943  -0.0431 Balanced, <0.1
    ## indigena      Binary  0.1392  -0.0023 Balanced, <0.1
    ## rural         Binary  0.0786   0.0105 Balanced, <0.1
    ## escoacum     Contin.  0.4698  -0.0749 Balanced, <0.1
    ## casadounion   Binary -0.1421   0.0268 Balanced, <0.1
    ## jefehog       Binary  0.0561   0.0128 Balanced, <0.1
    ## haymenores    Binary -0.0931   0.0012 Balanced, <0.1
    ## proggob       Binary  0.0973   0.0035 Balanced, <0.1
    ## tot_integ    Contin. -0.0853   0.0246 Balanced, <0.1
    ## 
    ## Balance tally for mean differences
    ##                    count
    ## Balanced, <0.1        10
    ## Not Balanced, >0.1     0
    ## 
    ## Variable with the greatest mean difference
    ##  Variable Diff.Adj    M.Threshold
    ##  escoacum  -0.0749 Balanced, <0.1
    ## 
    ## Sample sizes
    ##           Control Treated
    ## All         14293     859
    ## Matched       859     859
    ## Unmatched   13434       0

    Finalmente, podemos representar esta información en un loveplot:

    m.out.a[["X"]][["factor(cve_ent)"]] <- NULL
    
    love.plot(bal.tab(m.out.a), threshold = 0.1)

    Concluimos que existe un balance apropiado de los covariables usando el modelo del PS propuesto y el tratamiento definido por la variable jcf2.

  3. [10 puntos] Para la probabilidad de encontrar empleo, encontro, se comparan a los beneficiarios de JCF con los jóvenes en general. Los beneficiarios tienen jcf==1, mientras que el resto de los jóvenes tienen jcf==0. Realice la estimación del TOT y la inferencia, de manera análoga a lo realizado en la parte a.

    Seguimos un procedimiento similar al anterior, pero teniendo en cuenta que la variable de tratamiento es jcf:

    sub.data <- data.jcf %>%
        filter(transicion == 1) %>%
        dplyr::select(encontro, jcf, mujer, indigena, cve_ent, rural, escoacum, casadounion, jefehog,
            haymenores, proggob, tot_integ, factor.x)
    
    sub.data <- sub.data[complete.cases(sub.data), ]
    
    m.out.b <- matchit(formula = jcf ~ mujer + indigena + factor(cve_ent) + rural + escoacum + casadounion +
        jefehog + haymenores + proggob + tot_integ, method = "nearest", ratio = 1, distance = "glm",
        replace = FALSE, data = sub.data)
    
    # Y entonces hacemos nuestra comparación
    z.out.b <- zelig(formula = encontro ~ jcf, data = match.data(m.out.b), model = "ls")
    ## How to cite this model in Zelig:
    ##   R Core Team. 2007.
    ##   ls: Least Squares Regression for Continuous Dependent Variables
    ##   in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
    ##   "Zelig: Everyone's Statistical Software," https://zeligproject.org/
    
    # Simularemos el valor esperado de las diferencias cuando t==0
    x.out.b <- setx(z.out.b, jcf = 0)
    
    # Con respecto a cuando t==1
    x1.out.b <- setx1(z.out.b, jcf = 1)
    
    # Corremos la simulación
    sim.out.b <- sim(z.out.b, x = x.out.b, x1 = x1.out.b)
    
    # Vemos los resultados
    summary(sim.out.b)
    ## 
    ##  sim x :
    ##  -----
    ## ev
    ##        mean        sd       50%     2.5%     97.5%
    ## 1 0.2582677 0.0235609 0.2582197 0.211883 0.3026533
    ## pv
    ##           mean       sd       50%       2.5%    97.5%
    ## [1,] 0.2485625 0.446642 0.2553932 -0.6632649 1.080309
    ## 
    ##  sim x1 :
    ##  -----
    ## ev
    ##        mean         sd       50%      2.5%    97.5%
    ## 1 0.4902346 0.02375526 0.4914075 0.4412101 0.535353
    ## pv
    ##           mean        sd       50%       2.5%    97.5%
    ## [1,] 0.4825499 0.4571934 0.4920552 -0.4678052 1.361199
    ## fd
    ##        mean         sd      50%      2.5%     97.5%
    ## 1 0.2319669 0.03370108 0.230831 0.1667102 0.2972108

    El TOT estimado es de un incremento de 23.45% (e.e 3.33) en la probabilidad de encontrar un empleo para los beneficiarios del programa, comparados con el resto de jóvenes. La magnitud del impacto estimado es muy parecida al ATE reportado en el informe, aunque la precisión es menor.

  4. [5 puntos] Evalúe qué tan bueno es el procedimiento de la parte c. en balancear las características observadas una vez realizado el matching. Cree un love plot y realice pruebas formales para contrastar las diferencias en características observables antes y después del matching.

    Veamos las SMDs y el love plot:

    m.out.b[["X"]][["factor(cve_ent)"]] <- NULL
    
    bal.tab(m.out.b, m.threshold = 0.1, un = T)
    ## Call
    ##  matchit(formula = jcf ~ mujer + indigena + factor(cve_ent) + 
    ##     rural + escoacum + casadounion + jefehog + haymenores + proggob + 
    ##     tot_integ, data = sub.data, method = "nearest", distance = "glm", 
    ##     replace = FALSE, ratio = 1)
    ## 
    ## Balance Measures
    ##                 Type Diff.Un Diff.Adj    M.Threshold
    ## distance    Distance  0.6471   0.0036 Balanced, <0.1
    ## mujer         Binary  0.0285  -0.0273 Balanced, <0.1
    ## indigena      Binary  0.1177  -0.0298 Balanced, <0.1
    ## rural         Binary  0.0757   0.0174 Balanced, <0.1
    ## escoacum     Contin.  0.2335   0.0604 Balanced, <0.1
    ## casadounion   Binary  0.0732   0.0099 Balanced, <0.1
    ## jefehog       Binary  0.0484   0.0174 Balanced, <0.1
    ## haymenores    Binary  0.0505  -0.0074 Balanced, <0.1
    ## proggob       Binary  0.0926   0.0323 Balanced, <0.1
    ## tot_integ    Contin. -0.0450  -0.0087 Balanced, <0.1
    ## 
    ## Balance tally for mean differences
    ##                    count
    ## Balanced, <0.1        10
    ## Not Balanced, >0.1     0
    ## 
    ## Variable with the greatest mean difference
    ##  Variable Diff.Adj    M.Threshold
    ##  escoacum   0.0604 Balanced, <0.1
    ## 
    ## Sample sizes
    ##           Control Treated
    ## All         35530     403
    ## Matched       403     403
    ## Unmatched   35127       0
    
    love.plot(bal.tab(m.out.b), threshold = 0.1)

    De nuevo parece haber un buen balance en los covariables.

  5. [5 puntos] Estime ahora el TOT en el ingreso trimestral, como en la parte a., pero usando un caliper de 0.1 y 3 vecinos a ser emparejados. ¿Cómo cambian sus resultados respecto a los de la parte a.?

    El procedimiento es:

    sub.data <- data.jcf %>%
        dplyr::select(ingtot_tri, jcf2, mujer, indigena, cve_ent, rural, escoacum, casadounion,
            jefehog, haymenores, proggob, tot_integ, factor.x)
    
    sub.data <- sub.data[complete.cases(sub.data), ]
    
    m.out.e <- matchit(formula = jcf2 ~ mujer + indigena + factor(cve_ent) + rural +
        escoacum + casadounion + jefehog + haymenores + proggob + tot_integ, method = "nearest",
        caliper = 0.1, ratio = 3, distance = "glm", replace = FALSE, data = sub.data)
    
    z.out.e <- zelig(formula = ingtot_tri ~ jcf2, data = match.data(m.out.a), model = "ls")
    ## How to cite this model in Zelig:
    ##   R Core Team. 2007.
    ##   ls: Least Squares Regression for Continuous Dependent Variables
    ##   in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
    ##   "Zelig: Everyone's Statistical Software," https://zeligproject.org/
    
    # Simulo las dos situaciones
    x.out.e <- setx(z.out.e, jcf2 = 0)
    x1.out.e <- setx1(z.out.e, jcf2 = 1)
    
    # Corremos la simulación
    sim.out.e <- sim(z.out.e, x = x.out.e, x1 = x1.out.e)
    
    # Vemos los resultados
    summary(sim.out.e)
    ## 
    ##  sim x :
    ##  -----
    ## ev
    ##       mean       sd      50%     2.5%    97.5%
    ## 1 1624.597 233.7063 1616.959 1173.134 2089.763
    ## pv
    ##          mean       sd      50%      2.5%    97.5%
    ## [1,] 1654.698 6918.868 1403.177 -11718.46 15172.09
    ## 
    ##  sim x1 :
    ##  -----
    ## ev
    ##       mean       sd      50%     2.5%    97.5%
    ## 1 10097.69 235.5765 10092.84 9635.253 10564.24
    ## pv
    ##          mean       sd      50%      2.5%    97.5%
    ## [1,] 10125.25 7049.094 10151.67 -3460.996 23491.42
    ## fd
    ##       mean       sd      50%     2.5%    97.5%
    ## 1 8473.093 329.8256 8470.301 7824.977 9119.494

    El efecto estimado es de 8466.431 pesos (e.e. 342.3131), apenas distinto de los 8473.637 encontrados en la parte a. El resultado estimado es robusto a usar varios vecinos y un caliper.

  6. [15 puntos] Proponga una estrategia de PSM para evaluar el efecto del programa en la probabilidad de encontrar empleo, superior a la de las partes c. y d. Para ello, escoja un modelo para el PS y un algoritmo de emparejamiento. El modelo del PS puede modificarse de diversas maneras: añadiendo polinomios de las variables continuas, incluyendo interacciones, agregando variables disponibles en los datos, o construyendo nuevas variables a partir de los datos en bruto, modificando el script limpiar_jcf_analisis.r. Reporte qué tan bueno es su procedimiento para construir grupos balanceados y compare los resultados que obtiene con los de las partes c. y d.

    Aquí cometí un error al escribir la pregunta, porque debería decir “… superior a la de las partes c. y d.” Mi idea era que mejoraran el modelo y el balance que se encontró para buscar el impacto en la empleabilidad. No quise deshacer lo que ya hubieran realizado hasta el momento en que me di cuenta. Consideraré este hecho cuando revise sus tareas, dando todo el crédito para lo que hayan hecho y que tenga sentido.

    Intento encontrar un mejor modelo para evaluar el efecto en la empleabilidad haciendo un par de modificaciones. La primera consiste en introducir un término cuadrático para las variables de escolaridad y tamaño del hogar. La segunda es especificar un mecanismo de emparejamiento genético.

    El algoritmo genético hace los emparejamientos no solo en términos del PS, sino que también pone atención a los covariables mismos. El algoritmo otorga pesos a los covariables y al PS, por lo que los métodos que vimos en clase son un caso particular en los que el peso a los covariables es cero y solo nos enfocamos en el PS. Este método fue introducido por Diamond y Sekhon (2013) y busca optimizar el balance de los covariables.

    La idea del método es especificar una función de distancia entre los covariables del grupo tratado y no tratado, \(GMD(W)\) por genetic matching distance, que otorga pesos \(W\) a cada covariable y al PS. El algoritmo comienza dando unos pesos iniciales \(W\) iguales a todos los covariables y al PS (los genes parentales). Luego, genera \(p_w\) mutaciones aleatoria a los pesos, dando lugar a una nueva generación de pesos \(W\) de tamaño \(p_w\) (pop.size en el código). Se calcula la distancia usando cada una de los candidatos a pesos, es decir, se calcula la función de pérdida \(p_w\) veces. Se escoge la \(W\) que minimiza la función de pérdida (selección natural) para una generación. Se repite el proceso comenzando con los pesos \(W\) que minimizan \(GMD(W)\) en la generación previa, hasta que se alcance el número máximo de iteraciones o el balance en los covariables deseado.

    Estimé este procedimiento usando un \(p_w=100\), lo cual tardó más o menos una hora en correr. Se recomienda usar un \(p_w\) varias veces más grandes que el que yo usé.

    
    sub.data <- data.jcf %>%
        filter(transicion == 1) %>%
        dplyr::select(encontro, jcf, mujer, indigena, cve_ent, rural, escoacum, casadounion,
            jefehog, haymenores, proggob, tot_integ, factor.x)
    
    sub.data <- sub.data[complete.cases(sub.data), ]
    
    m.out.f <- matchit(formula = jcf ~ mujer + indigena + factor(cve_ent) + rural + escoacum +
        escoacum^2 + casadounion + jefehog + haymenores + proggob + tot_integ^2, method = "genetic",
        pop.size = 100, distance = "glm", replace = FALSE, data = sub.data)
    
    
    # Corroboro balance
    m.out.f[["X"]][["factor(cve_ent)"]] <- NULL
    
    bal.tab(m.out.f, m.threshold = 0.1, un = T)
    ## Call
    ##  matchit(formula = jcf ~ mujer + indigena + factor(cve_ent) + 
    ##     rural + escoacum + escoacum^2 + casadounion + jefehog + haymenores + 
    ##     proggob + tot_integ^2, data = sub.data, method = "genetic", 
    ##     distance = "glm", replace = FALSE, pop.size = 100)
    ## 
    ## Balance Measures
    ##                 Type Diff.Un Diff.Adj    M.Threshold
    ## distance    Distance  0.6471   0.0173 Balanced, <0.1
    ## mujer         Binary  0.0285  -0.0025 Balanced, <0.1
    ## indigena      Binary  0.1177  -0.0074 Balanced, <0.1
    ## rural         Binary  0.0757   0.0050 Balanced, <0.1
    ## escoacum     Contin.  0.2335   0.0105 Balanced, <0.1
    ## casadounion   Binary  0.0732   0.0074 Balanced, <0.1
    ## jefehog       Binary  0.0484   0.0025 Balanced, <0.1
    ## haymenores    Binary  0.0505   0.0000 Balanced, <0.1
    ## proggob       Binary  0.0926   0.0025 Balanced, <0.1
    ## tot_integ    Contin. -0.0450  -0.0210 Balanced, <0.1
    ## 
    ## Balance tally for mean differences
    ##                    count
    ## Balanced, <0.1        10
    ## Not Balanced, >0.1     0
    ## 
    ## Variable with the greatest mean difference
    ##   Variable Diff.Adj    M.Threshold
    ##  tot_integ   -0.021 Balanced, <0.1
    ## 
    ## Sample sizes
    ##           Control Treated
    ## All         35530     403
    ## Matched       403     403
    ## Unmatched   35127       0
    
    love.plot(bal.tab(m.out.f), threshold = 0.1)

    
    
    z.out.f <- zelig(formula = encontro ~ jcf, data = match.data(m.out.f), model = "ls")
    ## How to cite this model in Zelig:
    ##   R Core Team. 2007.
    ##   ls: Least Squares Regression for Continuous Dependent Variables
    ##   in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
    ##   "Zelig: Everyone's Statistical Software," https://zeligproject.org/
    
    # Simulación
    x.out.f <- setx(z.out.f, jcf = 0)
    x1.out.f <- setx1(z.out.f, jcf = 1)
    sim.out.f <- sim(z.out.f, x = x.out.f, x1 = x1.out.f)
    summary(sim.out.f)
    ## 
    ##  sim x :
    ##  -----
    ## ev
    ##        mean         sd     50%      2.5%    97.5%
    ## 1 0.2629635 0.02324159 0.26217 0.2167832 0.308694
    ## pv
    ##          mean        sd       50%       2.5%    97.5%
    ## [1,] 0.253472 0.4713359 0.2748007 -0.6675696 1.130598
    ## 
    ##  sim x1 :
    ##  -----
    ## ev
    ##        mean        sd       50%      2.5%     97.5%
    ## 1 0.4919948 0.0237728 0.4913797 0.4469201 0.5395469
    ## pv
    ##           mean        sd       50%       2.5%    97.5%
    ## [1,] 0.5153231 0.4736652 0.5132267 -0.3829502 1.479568
    ## fd
    ##        mean         sd       50%      2.5%     97.5%
    ## 1 0.2290312 0.03411491 0.2285096 0.1642413 0.2969145

    Notemos que el balance mejora, como es claro en el love plot: los puntos azules están muy cerca del cero. El efecto estimado, sin embargo, no es muy distinto al encontrado en la parte c.

    Una de las motivaciones de hacer esta pregunta es que, dado que tenemos dos definiciones diferentes de tratamiento, jcf y jcf2, nada nos garantiza que el balance se adecuado si mantenemos fijo el modelo del PS. Esto lo fui pensando al leer la tesina de su compañera Daniela Díaz, quien evaluó el efecto de tener distintos productos financieros en medidas de vulnerabilidad. Ella modificó el modelo PS para cada uno de los tratamientos, definidos como tener o no cierto producto financiero.

Pregunta 2

Suponga que se convierte en asesor de la instancia gubernamental encargada de la seguridad alimentaria. Al gobierno le interesa que la seguridad alimentaria de las familias productoras de maíz para autoconsumo no se vea afectada negativamente por la presencia de cierta plaga y dará una transferencia per cápita a todos los pequeños productores de maíz cuyos cultivos se considere están afectados por dicha plaga. Para determinar qué hogares reciben la transferencia se decide usar un índice de prevalencia de la plaga y se selecciona un umbral por arriba del cual está demostrado que los rendimientos del cultivo del maíz se ven seriamente afectados. Esta inspección se llevará a cabo por autoridades federales y el umbral es conocido solo por estas autoridades. Cuando se determine que la prevalencia está por encima del umbral, el monto del programa será transferido de manera inmediata, electrónicamente.

  1. [5 puntos] ¿Qué aspectos del programa permitirían emplear un diseño de regresión discontinua para evaluar la efectividad de este sobre la seguridad alimentaria y cómo mostraría su validez empíricamente?

    En este caso podemos usar el método de regresión discontinua por las siguientes razones:

    i. La variable de selección es continua.

    ii. Es estatus de tratamiento es una función determinística de la posición de la variable de selección respecto al umbral.

    iii. La probabilidad de recibir el tratamiento es discontinua en el umbral.

    iv. Los productores no pueden manipular la prevalencia de la plaga para posicionarse estratégicamente por encima del umbral.

  2. [5 puntos] ¿Cómo emplearía el diseño de este programa para evaluar su efectividad con un modelo de regresión discontinua nítida? Elabore una gráfica donde explique una situación en la que el programa muestra ser efectivo. Describa cómo usaría una regresión para hacer inferencia respecto a la efectividad del programa.

    La forma gráfica de inspeccionar la presencia de una regresión consiste en graficar la variable de resultados en función de la variable de asignación. En este caso, esperaríamos que las familias que están por encima del umbral tengan una diferencia notable en términos de seguridad alimentaria si la transferencia empleada se usa para comprar alimentos. No era estrictamente necesario simular un proceso para obtener una representación gráfica, pero aquí lo hice así. Quizás esto pueda ser de utilidad para futuras aplicaciones:

    set.seed(1711)
    
    plaga <- runif(1000, -1, 1)
    y <- 3 + 2 * plaga + 5 * (plaga >= 0) + rnorm(1000, mean = 0, sd = 0.2)
    
    data.sharp <- data.frame(y, plaga, c = 0)
    
    
    data.sharp %>%
        ggplot() + geom_point(aes(x = plaga, y = y), size = 0.5, alpha = 0.5) + geom_abline(intercept = 3,
        slope = 2, linetype = "dashed") + geom_abline(intercept = 8, slope = 2, linetype = "dashed") +
        geom_vline(xintercept = 0)

    Un caso que no es una discontinuidad difusa es el representado por la siguiente situación:

    set.seed(1711)
    
    plaga <- runif(1000, -1, 1)
    y <- ((4 + 0.5 * plaga) * (plaga >= 0)) + ((4 + 3 * plaga) * (plaga < 0)) + rnorm(1000,
        mean = 0, sd = 0.2)
    
    data.kink <- data.frame(y, plaga, c = 0)
    
    data.kink %>%
        ggplot() + geom_point(aes(x = plaga, y = y), size = 0.5, alpha = 0.5) + geom_abline(intercept = 4,
        slope = 0.5, linetype = "dashed") + geom_abline(intercept = 4, slope = 3, linetype = "dashed") +
        geom_vline(xintercept = 0)

    A esto se le conoce como regresión con pliegues (o regression kinks). La idea es que en \(x_0\) ocurre un cambio en la pendiente de la función de esperanza de \(y\). Card et al. (2016) describen la teoría y la práctica de estos diseños que no alcanzamos a cubrir en clase.

    Paramétricamente, la forma más sencilla de identificar el efecto de la discontinuidad es especificando una regresión como sigue: \[y_i=\alpha+\tau D_i+ \beta x_i+\varepsilon_i\] donde \(x_i\) es la variable de selección y \(D_i\) es una variable indicadora que toma el valor de uno cuando el índice de prevalencia de la plaga rebasa el umbral. Controlar por \(x_i\) captura la relación que tiene la prevalencia de la plaga en la seguridad alimentaria, por ejemplo, vía los rendimientos. Se recomiendan al menos dos tipos de procedimientos más completos para comprobar la robustez de los efectos encontrados.

    El primero es incluir un polinomio lo suficientemente flexible de \(x_i\): \[y_i=\alpha+\tau D_i+ Bf(x_i)+\varepsilon_i\]

    El segundo consiste en permitir que la pendiente sea diferente antes y después de la discontinuidad: \[y_i=\alpha+\tau D_i+ +\beta_0(x_i-x_0)+\beta_1(x_i-x_0)D_i+\varepsilon_i\]

    Más aún, es posible combinar estas dos posibilidades para dar lugar a modelos más flexibles. Se espera que las conclusiones sean robustas al uso de modelos extremadamente complejos.

  3. [5 puntos] ¿Qué factores podrían invalidar el uso de este método para evaluar el programa?

    La principal preocupación es la posibilidad de manipulación de la prevalencia de la plaga para que la medición lo clasifique como receptor del programa. Podemos pensar en situaciones donde esto pudiera suceder con un individuo altamente sofisticado que pudiera manipular la presencia de la plaga de forma estratégica. Pensando que esto es costoso, el individuo estratégicamente debería seleccionar un punto justo por encima del umbral. Aunque difícil de suceder esta posibilidad podría investigarse empíricamente, por ejemplo, verificando que no haya “amontonamientos” justo por encima de la discontinuidad.

    Si existiera corrupción y muchos no elegibles recibieran la transferencia o si las familias no gastaran la transferencia en alimentos que mejoren su seguridad alimentaria el diseño también estaría comprometido.

  4. Suponga que otro de los asesores juzga como demasiado paternalista la transferencia y propone que, en su lugar, se otorgue un cupón válido para canjearse por bultos de un plaguicida. Asumiendo que en una encuesta posterior usted podría conocer la cantidad precisa de plaguicida aplicado, ¿cómo emplearía un diseño de regresión discontinua difusa para evaluar el efecto del uso del plaguicida sobre la seguridad alimentaria? En particular, describa:

    1. [5 puntos] ¿Cómo estimaría la forma reducida? ¿Cuál es el coeficiente relevante y cuál es su interpretación?

      El problema puede ser visto entonces como un diseño de regresión discontinua difusa. La discontinuidad define la intensidad del tratamiento, en este caso dada por la cantidad de plaguicida efectivamente aplicado. La forma reducida se estima con una regresión de la variable de resultados sobre el instrumento. Al igual que cuando se estudió la interpretación del LATE, este coeficiente da la correlación entre la seguridad alimentaria y el estado del tratamiento, pero no toma en cuenta que la seguridad alimentaria también depende de la cantidad de plaguicida usado, una decisión endógena.

    2. [5 puntos] ¿Cómo estimaría la primera y la segunda etapa? ¿Cuáles son los coeficientes relevantes y cuál es su interpretación?

      La primera etapa consiste en estimar la relación entre la variable endógena y el instrumento. En este caso, el instrumento es una variable indicadora que toma valor de 1 si la prevalencia de la plaga rebasa el umbral. La decisión endógena es la cantidad de plaguicida empleado. Se estima por una regresión de la variable endógena en función del instrumento.

      La segunda etapa consiste en estimar el efecto sobre la seguridad alimentaria de la cantidad plaguicida que predice el instrumento. Conceptualmente es como si se corriera una regresión de la variable de seguridad alimentaria en función de los valores ajustados en la primera etapa de la cantidad de plaguicida empleado. En la práctica, nunca se estiman dos regresiones separadas, sino que se usa la definición del estimador de mínimos cuadrados en dos etapas. El coeficiente es el efecto del uso de plaguicida en la seguridad alimentaria.

    3. [5 puntos] ¿Cuáles son los supuestos necesarios para estimar este modelo usando mínimos cuadrados en dos etapas?

      Los supuestos econométricos para la estimación del modelo de regresión discontinua difusa son los mismos que para cualquier otro problema de variables instrumentales: 1) Exclusión: el instrumento no pertenece a la ecuación estructural; y 2) Relevancia de la primera etapa: el instrumento está correlacionado con la variable endógena.

Pregunta 3

La base de datos headstar.csv contiene información de 2,810 condados de los Estados Unidos. La variable mort_age59_related_postHS indica la mortalidad infantil en cada uno de los condados. El programa Head Star otorgó fondos de su componente de salud a todos los condados con un índice de pobreza superior a 59.1968. La variable povrate60 es el índice de pobreza para cada condado. Se desea estimar el efecto del programa en la mortalidad infantil empleando un diseño de regresión discontinua.

  1. [5 puntos] Genere una gráfica donde muestre evidencia de una discontinuidad en la tasa de mortalidad para aquellos condados que no recibieron fondos del programa.

    Pomdemos usar las técnicas vistas en clase o, alternativamente, el paquete rdplot. Aquí uso rdplot con un polinomio de orden 1 y de orden 2, respectivamente:

    data.hs <- read_csv("./headstar.csv", locale = locale(encoding = "latin1"))
    
    (rdplot(y = data.hs$mort_age59_related_postHS, x = data.hs$povrate60, c = 59.1968,
        p = 1))

    ## Call: rdplot
    ## 
    ## Number of Obs.                 2783
    ## Kernel                      Uniform
    ## 
    ## Number of Obs.                 2489             294
    ## Eff. Number of Obs.            2489             294
    ## Order poly. fit (p)               1               1
    ## BW poly. fit (h)             43.988          22.373
    ## Number of bins scale          1.000           1.000
    
    
    (rdplot(y = data.hs$mort_age59_related_postHS, x = data.hs$povrate60, c = 59.1968,
        p = 2))

    ## Call: rdplot
    ## 
    ## Number of Obs.                 2783
    ## Kernel                      Uniform
    ## 
    ## Number of Obs.                 2489             294
    ## Eff. Number of Obs.            2489             294
    ## Order poly. fit (p)               2               2
    ## BW poly. fit (h)             43.988          22.373
    ## Number of bins scale          1.000           1.000
  2. [5 puntos] Estime la versión más básica de un modelo de regresión discontinua. Reporte el coeficiente estimado del efecto del tratamiento y su significancia estadística. Use una ventana de 10 puntos en el índice de pobreza antes y después del corte. Interprete su resultado.

    Podemos estimar el salto en la discontinuidad como:

    data.hs <- data.hs %>%
        mutate(ispoor = ifelse(povrate60 >= 59.1968, 1, 0))
    
    summary(lm(mort_age59_related_postHS ~ povrate60 + ispoor, data = filter(data.hs,
        povrate60 >= 59.1968 - 10 & povrate60 <= 59.1968 + 10)))
    ## 
    ## Call:
    ## lm(formula = mort_age59_related_postHS ~ povrate60 + ispoor, 
    ##     data = filter(data.hs, povrate60 >= 59.1968 - 10 & povrate60 <= 
    ##         59.1968 + 10))
    ## 
    ## Residuals:
    ##    Min     1Q Median     3Q    Max 
    ## -3.479 -2.905 -2.331  1.774 61.686 
    ## 
    ## Coefficients:
    ##             Estimate Std. Error t value Pr(>|t|)  
    ## (Intercept) -1.63212    4.43998  -0.368   0.7133  
    ## povrate60    0.08643    0.08188   1.056   0.2916  
    ## ispoor      -1.53264    0.89167  -1.719   0.0862 .
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ## 
    ## Residual standard error: 5.409 on 568 degrees of freedom
    ##   (4 observations deleted due to missingness)
    ## Multiple R-squared:  0.006279,   Adjusted R-squared:  0.00278 
    ## F-statistic: 1.795 on 2 and 568 DF,  p-value: 0.1671

    Esto es, aproximadamente una reducción de 1.53 puntos en la mortalidad infantil.

  3. [5 puntos] Estime la misma especificación que en la parte b., pero ahora con una ventana de 5 puntos en el índice de pobreza. Interprete sus resultados.

    summary(lm(mort_age59_related_postHS ~ povrate60 + ispoor, data = filter(data.hs,
        povrate60 >= 59.1968 - 5 & povrate60 <= 59.1968 + 5)))
    ## 
    ## Call:
    ## lm(formula = mort_age59_related_postHS ~ povrate60 + ispoor, 
    ##     data = filter(data.hs, povrate60 >= 59.1968 - 5 & povrate60 <= 
    ##         59.1968 + 5))
    ## 
    ## Residuals:
    ##    Min     1Q Median     3Q    Max 
    ## -3.880 -2.929 -2.249  2.125 61.481 
    ## 
    ## Coefficients:
    ##             Estimate Std. Error t value Pr(>|t|)  
    ## (Intercept) -11.0201    12.6116  -0.874    0.383  
    ## povrate60     0.2520     0.2227   1.132    0.259  
    ## ispoor       -2.2562     1.3100  -1.722    0.086 .
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ## 
    ## Residual standard error: 5.806 on 306 degrees of freedom
    ##   (1 observation deleted due to missingness)
    ## Multiple R-squared:  0.01116,    Adjusted R-squared:  0.004697 
    ## F-statistic: 1.727 on 2 and 306 DF,  p-value: 0.1796

    El efecto estimado es de 2.26 puntos menos en la mortalidad infantil, estimado con una mayor precisión que en la parte b.

  4. [5 puntos] Regrese a una ventana de 10 puntos como en la parte b., pero ahora incluya un polinomio de grado 2 para el índice de pobreza y permita un coeficiente distinto para el índice de pobreza antes y después del corte. Interprete sus resultados.

    summary(lm(mort_age59_related_postHS ~ povrate60 + I(povrate60^2) + ispoor + povrate60 *
        ispoor + I(povrate60^2) * ispoor, data = filter(data.hs, povrate60 >= 59.1968 -
        10 & povrate60 <= 59.1968 + 10)))
    ## 
    ## Call:
    ## lm(formula = mort_age59_related_postHS ~ povrate60 + I(povrate60^2) + 
    ##     ispoor + povrate60 * ispoor + I(povrate60^2) * ispoor, data = filter(data.hs, 
    ##     povrate60 >= 59.1968 - 10 & povrate60 <= 59.1968 + 10))
    ## 
    ## Residuals:
    ##    Min     1Q Median     3Q    Max 
    ## -3.767 -2.853 -2.552  1.868 61.665 
    ## 
    ## Coefficients:
    ##                         Estimate Std. Error t value Pr(>|t|)
    ## (Intercept)             90.78703  113.33949   0.801    0.423
    ## povrate60               -3.29138    4.18847  -0.786    0.432
    ## I(povrate60^2)           0.03077    0.03861   0.797    0.426
    ## ispoor                -277.33745  240.62148  -1.153    0.250
    ## povrate60:ispoor         9.06706    7.87095   1.152    0.250
    ## I(povrate60^2):ispoor   -0.07481    0.06496  -1.152    0.250
    ## 
    ## Residual standard error: 5.415 on 565 degrees of freedom
    ##   (4 observations deleted due to missingness)
    ## Multiple R-squared:  0.009406,   Adjusted R-squared:  0.0006398 
    ## F-statistic: 1.073 on 5 and 565 DF,  p-value: 0.3743

    El efecto estimado es un raro -277.33, pero no significativo. El modelo más sencillo parece darnos los resultados más sensatos. De hecho, esta pregunta estuvo basada en un artículo de Calonico et al. (2019), en el que los autores desarrollan la teoria para introducir covariables en los diseños con discontinuidades. El impacto estimado es de entre -2.51 y -2.41, parecido a lo obtenido en la parte c.

Previous