Taller 5: Variables Instrumentales II

Carlos Daboin Contreras (cdaboincontreras@udesa.edu.ar)

2023-10-01

Segunda semana de IV

Librerias

library(tidyverse)    # manipulacion y visualizacion de datos
library(haven)        # leer datos en formato dta.
library(modelsummary) # para resultados
library(gt)           # para tablas
library(AER)          # ivreg
library(patchwork) #sirve para pegar graficas de lado a lado.

Angrist & Krueger (1991)

Does Compulsory School Attendance Affect Schooling and Earnnings?

La pregunta

¿Cúal es el efecto de la obligatoriedad de la escuela sobre las tasas de graduación y los ingresos futúros de los jóvenes?

  • Como funcionarios e investigadores no puedemos obligar a la gente a finalizar el secundario, mucho menos impedirleslo.
  • Eso implica que no podemos diseñar asignaciones experimentales en esta área.
  • Usar datos observacionales (i.e, la EPH) pudiera traernos problemas. ¿Por qué?
  • ¿Que tal si pudieramos encontrar grupos de personas iguales (en promedio) en varias dimesiones con la exepción de la educación que recibieron?

Angrist & Krueger (1991)

Revisan la historia y se dan cuenta de algo:

  • La educación secundaria no era obligatoria para los mayores de 16 añps en las décadas de 1930 y 1940. Los mayores de 16 podían retirarse si querían.

  • Aquellos estudiantes de último año de secundarion que cumplieron años antes de recibirse (en Q1 y Q2) tenían mayores chances de retirarse que el resto.

  • El trimestre (Q) de nacimiento puede ser un instrumento útil para capturar la variablidad exógena en el nivel educativo: No hay diferencias entre los nacido en Q1-Q2 y los nacidos en Q3-Q4.

Quarter of birth and education data

## Recuerden cambiar esta línea de acuerdo a su working directory y al lugar donde guardaron el archivo QOB.dta
qob_data <- read_dta("../data/QOB.dta")

# Que variables tenemos
names(qob_data)
 [1] "age"     "ageq"    "ageqsq"  "educ"    "lwage"   "married" "census" 
 [8] "qob"     "race"    "smsa"    "yob"    
# Estadisticas descriptivas generales
summary(qob_data)
      age             ageq           ageqsq            educ      
 Min.   :30.00   Min.   :30.25   Min.   : 915.1   Min.   : 0.00  
 1st Qu.:33.00   1st Qu.:33.50   1st Qu.:1122.2   1st Qu.:12.00  
 Median :38.00   Median :38.00   Median :1444.0   Median :12.00  
 Mean   :38.49   Mean   :38.86   Mean   :1544.2   Mean   :13.25  
 3rd Qu.:43.00   3rd Qu.:43.75   3rd Qu.:1914.1   3rd Qu.:16.00  
 Max.   :50.00   Max.   :50.00   Max.   :2500.0   Max.   :20.00  
     lwage           married           census        qob       
 Min.   :-2.342   Min.   :0.0000   Min.   :80   Min.   :1.000  
 1st Qu.: 5.583   1st Qu.:1.0000   1st Qu.:80   1st Qu.:2.000  
 Median : 5.922   Median :1.0000   Median :80   Median :3.000  
 Mean   : 5.839   Mean   :0.8277   Mean   :80   Mean   :2.525  
 3rd Qu.: 6.186   3rd Qu.:1.0000   3rd Qu.:80   3rd Qu.:4.000  
 Max.   :11.225   Max.   :1.0000   Max.   :80   Max.   :4.000  
      race              smsa             yob       
 Min.   :0.00000   Min.   :0.0000   Min.   :30.00  
 1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:36.00  
 Median :0.00000   Median :0.0000   Median :42.00  
 Mean   :0.08166   Mean   :0.1858   Mean   :40.76  
 3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:46.00  
 Max.   :1.00000   Max.   :1.0000   Max.   :49.00

Quarter of birth and education data

  • Los datos son tidy y están a nivel de individuo.

  • Los anos de nacimiento van de 1930 hasta 1949.

  • No hay variable género (sólo varones).

  • No hay mucha diferencia entre las variables de edad y las variables de edad medidas en trimestres.[^3]

QoB como instrumento

stats<-qob_data %>% 
  filter(yob >=30 & yob <=49) %>%
  mutate(decade=case_when(yob>=30 & yob<=39~"30's",
                          TRUE ~ "40's"),
         wage=exp(lwage)) %>%
  # promedios
  group_by(decade,qob) %>% 
  summarise(mean_educ=mean(educ),
            mean_wage=mean(wage))

table1<-gt(stats) %>% 
  gt::cols_label(
    decade="Decada",
    qob="Trimestre de nacimiento",
    mean_educ="Anos de educacion promedio",
    mean_wage="Salario promedio"
  )
Trimestre de nacimiento Anos de educacion promedio Salario promedio
30's
1 12.68807 436.1368
2 12.74472 438.5450
3 12.80541 441.2872
4 12.83943 441.8054
40's
1 13.50922 391.5326
2 13.56706 391.4990
3 13.58315 391.8774
4 13.62326 388.2357

QoB como instrumento

  • Nacidos en ultimos trimestres tienen mayor promedio de escolaridad (años).
  • Tabién tienen mayores salarios promedio.
  • Las diferencias en el ingreso parecen desaparecer en la década de 1940.
  • Conviene remover los 40s.
Trimestre de nacimiento Anos de educacion promedio Salario promedio
30's
1 12.68807 436.1368
2 12.74472 438.5450
3 12.80541 441.2872
4 12.83943 441.8054
40's
1 13.50922 391.5326
2 13.56706 391.4990
3 13.58315 391.8774
4 13.62326 388.2357

Efecto QoB Visualmente

# Summarize the average wages by age:
data_age <- qob_data %>%
  group_by(qob, yob) %>%
  summarise(lnw = mean(lwage), m_educ = mean(educ)) %>%
  mutate(q4 = (qob == 4)) %>% 
  arrange( yob, qob)

p1<-ggplot(data_age, aes(x = yob + (qob - 1) / 4, y = m_educ)) +
      geom_line() +
      geom_point(mapping = aes(color = ifelse(q4==1,"4th Quarter","Rest"))) +
      geom_text(
        mapping = aes(label=qob),
        size=2, nudge_y = 0.05
      )+
      scale_x_continuous("Year of birth") +
      scale_y_continuous("Years of schooling") +
      labs(color="Quarter")
    
p2<-ggplot(data_age, aes(x = yob + (qob - 1) / 4, y = lnw)) +
      geom_line() +
      geom_point(mapping = aes(color = ifelse(q4==1,"4th Quarter","Rest"))) +
        geom_text(
        mapping = aes(label=qob),
        size=2, nudge_y = 0.01
      )+
      scale_x_continuous("Year of birth" ) +
      scale_y_continuous("Log weekly wages") +
      labs(color="Quarter")

p1+p2+patchwork::plot_layout(guides = "collect")

Efecto QoB Visualmente

Estimaciones via MCO

qob_data_30<-qob_data %>%
  filter(yob >=30 & yob <=39)

ols<-lm(lwage~educ + ageq+ageqsq+race+married+smsa, 
        data=qob_data_30) 
broom::tidy(ols) %>% select(term,estimate,std.error,statistic)
# A tibble: 7 × 4
  term          estimate std.error statistic
  <chr>            <dbl>     <dbl>     <dbl>
1 (Intercept)  4.97       0.293       17.0  
2 educ         0.0642     0.000338   190.   
3 ageq        -0.00499    0.0130      -0.383
4 ageqsq       0.0000913  0.000144     0.632
5 race        -0.269      0.00403    -66.8  
6 married      0.243      0.00318     76.4  
7 smsa        -0.195      0.00282    -69.3  
broom::glance(ols) %>% select(r.squared, adj.r.squared)
# A tibble: 1 × 2
  r.squared adj.r.squared
      <dbl>         <dbl>
1     0.157         0.157

Estimaciones via 2SLS (IV)

Primero creamos variables dummy para cada trimestre de nacimiento :

qob_data_30<-qob_data_30 %>% 
  mutate(z1=ifelse(qob==1,1,0),
         z2=ifelse(qob==2,1,0),
         z3=ifelse(qob==3,1,0)) # solo necesito 3 ¿recuerdan por que?

Antes de seguir recuerda. Instrumentos deben ser

  • Relevantes (chequeable).

  • Válidos (inchequeable).

Test de instrumentos débiles

Se usa para verificar relevancia del instrumento.

Se rechaza la \(H_0\) de que los trimestres tienen un impacto conjunto nulo sobre la variabilidad de la educación.

ols_test<-lm(educ~z1+z2+z3, data = qob_data_30) 

broom::glance(ols_test) %>% select(statistic, p.value)
# A tibble: 1 × 2
  statistic  p.value
      <dbl>    <dbl>
1      34.0 5.74e-22

Rule of thumb: Si el F-statistic de la primera etapa es mayor a 10 los instrumentos son relevantes.

Estimaciones 2SLS

### Output first stage
first_stage<-lm(educ~z1+z2+z3+
                  ageq+ageqsq+race+married+smsa, 
                data = qob_data_30)

tidy(first_stage)
# A tibble: 9 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept) 13.3      1.51         8.83  1.08e-18
2 z1          -0.0930   0.0160      -5.81  6.33e- 9
3 z2          -0.0561   0.0160      -3.49  4.77e- 4
4 z3          -0.0236   0.0157      -1.50  1.32e- 1
5 ageq         0.0508   0.0672       0.756 4.50e- 1
6 ageqsq      -0.00127  0.000744    -1.70  8.88e- 2
7 race        -1.71     0.0206     -83.1   0       
8 married      0.157    0.0164       9.62  6.62e-22
9 smsa        -1.17     0.0144     -81.2   0       
### Output second stage
qob_data_30$educ_hat_over<-first_stage$fitted.values

second_stage<-lm(lwage~educ_hat_over+ ageq + ageqsq + race + married + smsa, 
                 data = qob_data_30)

tidy(second_stage)
# A tibble: 7 × 5
  term           estimate std.error statistic   p.value
  <chr>             <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    3.98      0.538        7.40  1.39e- 13
2 educ_hat_over  0.138     0.0330       4.18  2.94e-  5
3 ageq          -0.00856   0.0138      -0.619 5.36e-  1
4 ageqsq         0.000184  0.000158     1.17  2.44e-  1
5 race          -0.143     0.0566      -2.53  1.14e-  2
6 married        0.231     0.00618     37.3   1.18e-304
7 smsa          -0.109     0.0387      -2.82  4.73e-  3

Estimaciones IV en un paso

iv_fit <- AER::ivreg(lwage ~ educ+ageq + ageqsq + race + married + smsa | 
                  z1 + z2 + z3 + ageq + ageqsq + race + married + smsa ,
                data = qob_data_30) 

summary(iv_fit)

Call:
AER::ivreg(formula = lwage ~ educ + ageq + ageqsq + race + married + 
    smsa | z1 + z2 + z3 + ageq + ageqsq + race + married + smsa, 
    data = qob_data_30)

Residuals:
     Min       1Q   Median       3Q      Max 
-9.29526 -0.28520  0.06063  0.36577  4.96270 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.9825151  0.5467481   7.284 3.25e-13 ***
educ         0.1379401  0.0335265   4.114 3.88e-05 ***
ageq        -0.0085581  0.0140335  -0.610  0.54197    
ageqsq       0.0001838  0.0001601   1.148  0.25099    
race        -0.1432307  0.0574831  -2.492  0.01271 *  
married      0.2309541  0.0062800  36.776  < 2e-16 ***
smsa        -0.1093169  0.0393044  -2.781  0.00541 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6668 on 329502 degrees of freedom
Multiple R-Squared: 0.03511,    Adjusted R-squared: 0.03509 
Wald test:  3668 on 6 and 329502 DF,  p-value: < 2.2e-16

Resultados

tabla2<-
  msummary(
  list(
    'OLS1'=ols,
    'First Stage'=first_stage, 
    'Second Stage'=second_stage,
    'IV'= iv_fit
    ),
  stars = c('*' = .1, '**' = .05, '***' = .01),
  coef_map = c( 
    '(Intercept)' = 'Constant',
    'educ' ='Educ',
    'educ_hat_over'='Pred Educ',
    'ageq'='Age',
    'ageqsq'='Age sq',
    'race'='Race',
    'married'='Married',
    'smsa'='SMSA'
  ),
  gof_omit = "AIC|BIC|Log|RMSE")
OLS1 First Stage Second Stage  IV
Constant 4.968*** 13.327*** 3.983*** 3.983***
(0.293) (1.510) (0.538) (0.547)
Educ 0.064*** 0.138***
(0.000) (0.034)
Pred Educ 0.138***
(0.033)
Age −0.005 0.051 −0.009 −0.009
(0.013) (0.067) (0.014) (0.014)
Age sq 0.000 −0.001* 0.000 0.000
(0.000) (0.001) (0.000) (0.000)
Race −0.269*** −1.709*** −0.143** −0.143**
(0.004) (0.021) (0.057) (0.057)
Married 0.243*** 0.157*** 0.231*** 0.231***
(0.003) (0.016) (0.006) (0.006)
SMSA −0.195*** −1.169*** −0.109*** −0.109***
(0.003) (0.014) (0.039) (0.039)
Num.Obs. 329509 329509 329509 329509
R2 0.157 0.042 0.064 0.035
R2 Adj. 0.157 0.042 0.064 0.035
F 10211.543 1785.255 3782.605
* p < 0.1, ** p < 0.05, *** p < 0.01

Conclusiones

  • Los estimados de retorno de la educación por el método IV duplican los obtenidos vía MCO (14% vs 6.4%).

  • Los coeficientes relevantes son significativos al 1%

2SLS en situaciones de imperfect compliance

Vouchers para acceder a mosquiteros

  • Salud quiere promover el uso de mosquiteros durante el verano para prevenir una suba de casos de dengue.

  • Se quiere evaluar el desempeño del programa en Salto Encantado, Misiones.

  • Se envía un voucher digital a los jefes de hogar que puede ser canjeado por un kit de mosquiteros y perticidas para el hogar. Se hace por sorteo.

  • Equipo pasa a censar el pueblo al final del programa. Chequean si usaron el vaucher y su estado de salud.

Datos

Estos son los datos recolectados (son símulados).

bed_nets <- read_csv("../data/bed_nets_observed_epp.csv") %>% 
  # Las variables son categoricas. Asignamos categorias 
  # ordinales con fct_relevel, donde "No bed net" sea la primera
  # eso hara que salga como control en la regresion.
  mutate(bed_net = fct_relevel(bed_net, "No bed net"))

head(bed_nets)
# A tibble: 6 × 3
  lotery  bed_net    health
  <chr>   <fct>       <dbl>
1 Control No bed net   24.6
2 Control No bed net   36.4
3 Voucher Bed net      49.1
4 Voucher Bed net      47.0
5 Control Bed net      58.6
6 Voucher Bed net      63.9

¿Podemos simplemente calcular diferencia de medias?

Al haber sido sorteado, grupos de control y tratamiento son comparables.

Por ello el efecto del programa fue:

tidy(lm(health~lotery,data=bed_nets))
# A tibble: 2 × 5
  term          estimate std.error statistic  p.value
  <chr>            <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)      40.9      0.444     92.1  0       
2 loteryVoucher     5.99     0.630      9.51 5.36e-21

Pues No! No hay compliance perfecto.

Compliance imperfecto

¿Todos los que fueron incentivados tenian mosquitero al momento de la revisión?

  • Si hay compliance perfecto podemos calcular la diferencia entre el score de salud promedio entre grupos de control y tratamiento.

  • Si no lo hay, tenemos que buscar el efecto del tratamiento en los compliers. Este efecto se conoce como el Local Average Treatment Effect (LATE).

Evidentemente no lo hay. Muchos que no recibieron voucher tenian mosquitero (always takers evidentes), y muchos que recibieron voucher no tenian (never takers evidentes).

bed_nets %>% 
  group_by(lotery,bed_net) %>% 
  count()
# A tibble: 4 × 3
# Groups:   lotery, bed_net [4]
  lotery  bed_net        n
  <chr>   <fct>      <int>
1 Control No bed net   808
2 Control Bed net      196
3 Voucher No bed net   388
4 Voucher Bed net      608

Lo que vemos del compliance I

set.seed(1234)
ggplot(bed_nets, aes(y = health, x = bed_net)) + 
  geom_point(aes(shape = bed_net), 
             alpha=0.5,
             position = position_jitter(height = NULL, width = 0.25)) + 
  facet_wrap(vars(lotery)) + 
  theme_bw()+
  labs(shape = "Compliance",
       x = NULL, y = "Health status")

Lo que vemos del compliance I

Lo que vemos del compliance II

Lo que no vemos

Imaginemos que tenemos la capacidad de leer las intenciones de todas las personas que participaron del programa piloto.

IV al rescate

  • Usa el voucher para aislar la variabildiad exógena en la adquisición de mosquiteros.

  • Luego usar esta variabilidad para encontrar diferencias en la salud promedio de los hogares con y sin mosquitero.

Requisitos:

  • Exclusión: Suponemos que los vouchers no tienen ninguna incidencia sobre la salud de los hogares por fuera de su rol en la adquisión de mosquiteros.

  • El instrumento es exógeno: La distribución de vouchers fue aleatoria. ¿Por que?

  • Debe ser relevante: Se verifica si hay correlación fuerte entre voucher y tenencia de mosquiteros.

# 1st Stage: Fijate que la Y es una variable logica (TRUE/FALSE)
glance(lm(bed_net=="Bed net"~lotery, data=bed_nets))
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.179         0.179 0.444      436. 7.88e-88     1 -1215. 2436. 2452.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

IV estimates (2SLS)

De confiar en los supuestos anteriores, podemos confiar en la siguiente estimación: Usar mosquitero mejoró el puntaje de salud de los hogares participantes en 14.42 puntos.

model_2sls <- ivreg(health ~ bed_net | lotery, 
                    data = bed_nets)
summary(model_2sls)

Call:
ivreg(formula = health ~ bed_net | lotery, data = bed_nets)

Residuals:
      Min        1Q    Median        3Q       Max 
-35.32685  -7.88569  -0.03803   7.25495  47.45599 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     38.1229     0.5668   67.26   <2e-16 ***
bed_netBed net  14.4212     1.2527   11.51   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.63 on 1998 degrees of freedom
Multiple R-Squared: 0.3473, Adjusted R-squared: 0.347 
Wald test: 132.5 on 1 and 1998 DF,  p-value: < 2.2e-16

IV estimates manualmente

El efecto de haber sido asignado a la loteria (ITT) es el promedio ponderado de el efecto causal en los compliers, el efecto de los mosquiteros en los always takers, y el efecto causal en los never takers

\[ \begin{aligned} \mathrm{ITT} & =\pi_{\mathrm{C}} \mathrm{CACE}+\pi_{\mathrm{A}} \mathrm{ATACE}+\pi_{\mathrm{N}} \mathrm{NTACE} \end{aligned} \] No hay efecto de los mosquiteros en los always takers y never takers.

\[ \begin{aligned} \mathrm{ITT} & =\pi_{\mathrm{C}} \mathrm{CACE}+\pi_{\mathrm{A}} \times 0+\pi_{\mathrm{N}} \times 0 \\ \mathrm{ITT} & =\pi_{\mathrm{C}} \mathrm{CACE} \end{aligned} \] Entonces solo debo:

\[ \mathrm{CACE}=\frac{\mathrm{ITT}}{\pi_{\mathrm{C}}} \]

Podemos calcular el ITT

\[ \mathrm{CACE}=\frac{\mathrm{ITT}}{\pi_{\mathrm{C}}} \]

(ITT_table<-subset(bed_nets) %>% 
  summarise(mean_health=mean(health),.by =lotery ) %>% 
  mutate(ITT=mean_health-lag(mean_health)))
# A tibble: 2 × 3
  lotery  mean_health   ITT
  <chr>         <dbl> <dbl>
1 Control        40.9 NA   
2 Voucher        46.9  5.99

Podemos calcular el \(\pi_C\)

Los que reciben voucher y tienen mosquitero son combinacion de AA y compliers.

\[ \pi_C + \pi_A = \frac{N. BedNet \& Voucher}{N. Voucher} \] La proporcion de AA es la misma en grupo de Voucher y Control:

\[ \pi_A = \frac{N.BedNet\&Control}{N. Control} \] Entonces:

\[ \pi_C = \frac{N. BedNet \& Voucher}{N. Voucher}- \frac{N.BedNet\&Control}{N. Control} \]

IV estimates manualmente

\[ \mathrm{CACE}=\frac{\mathrm{ITT}}{\pi_{\mathrm{C}}} \]

(ITT<-ITT_table$ITT[2])
[1] 5.987992
bed_nets %>% 
  group_by(lotery, bed_net) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n)) %>% 
  filter(bed_net=="Bed net")
# A tibble: 2 × 4
# Groups:   lotery [2]
  lotery  bed_net     n  prop
  <chr>   <fct>   <int> <dbl>
1 Control Bed net   196 0.195
2 Voucher Bed net   608 0.610
(pi_c <- 0.6104418 - 0.1952191)
[1] 0.4152227
(LATE=ITT/pi_c)
[1] 14.42116

Otros diseños populares de IV

  • Sorteos o Loterías (ya lo vimos): Son excluyentes por definición pero deben ser relevantes.

  • Encouragement design (invito aleatoreamente a personas a usar un servicio): Debe ser relevante y excluyente. Hay que tener cuidado con el mensaje de aliento no rompa la exogeneidad.

  • Discontinuidades: Lo veremos en el bloque de Regression Discontinuity Design (RDD)

Fin

–> –> –> –> –> –> –> –> –>