Instrumental Variables

Welcome

# These are the libraries we will use today. Make sure to install them in your console in case you have not done so previously.

library(tidyverse) # To load the collection of packages in the tidyverse
library(dplyr) # To wrangle our data (this package is also loaded by library(tidyverse))
library(readr) # To load the .csv data (this package is also loaded by library(tidyverse))
library(ggplot2) # To create plots (this package is also loaded by library(tidyverse))
library(tidyr) # To use pivot_wider() (this package is also loaded by library(tidyverse))
library(ggdag) # To dagify and plot our DAG objects in R
#library(summarytools) # for ctable()
library(broom) # To format regression output
library(janitor) # To examine our data with tabyl()
library(stargazer) # To format model output
library(knitr) # To create HTML tables with kable()
library(kableExtra) # To format the HTML tables
library(AER) # To run 2SLS with ivreg()

Measuring the effect of mosquito net use on malaria infections

Imagine that the organization you work for is laying out a project to distribute mosquito nets to help combat malaria transmitions.

The funding agency requires a impact evaluation report from your organization. You are in charge of running the evaluation of this program.

You realize that there are potential unobserved counfounders that could bias the observed differences in malaria risk for mosquito net users and non-users. You also think about the ethical considerations of fully randomizing who receives the nets, so you remember your Statistics II lecture on IVs and set up an encouragement design.

Your benificiaries are scattered across ten villages. You decide to randomly select five villages to send SMS reminders every night encouraging them to use the mosquito nets. (This example is using simulated data)

Exploring our data

evaluation_df <- readr::read_csv("https://raw.githubusercontent.com/seramirezruiz/stats-ii-lab/master/Session%205/data/evaluation_data.csv") # loading simulated data frame of the intervention

You receive the results of your intervention from the M&E officers. There are 1000 inhabitants across the ten villages. This is what the data look like:

• village_name: A character string with the name of the village
• sms: A binary marker for the SMS encouragement
• net_use: A binary marker for mosquito net use
• malaria: A binary marker for malaria infection

Compliance types

You may remember, this table from the lecture:

We can crosstabulate our data with janitor::tabyl() and the additional features of the janitor::adorn_* functions.

Why janitor::tabyl()? Because as prospective policy analysts we will do a lot of counting.

As the vignette of the package even puts it:

Analysts do a lot of counting. Indeed, it’s been said that 'data science is mostly counting things.' But the base R function for counting, table(), leaves much to be desired:

• It doesn’t accept data.frame inputs (and thus doesn’t play nicely with the %>% pipe)
• It doesn’t output data.frames
• Its results are hard to format. Compare the look and formatting choices of an R table to a Microsoft Excel PivotTable.

This is how it works. Say we are interested in exploring the number of persons in each of the observed strata, we would do:

evaluation_df %>% # your data frame
  janitor::tabyl(net_use, sms) %>% # the two dimensions for the table (D, Z)
  janitor::adorn_totals(c("row", "col")) %>% # add totals for rows and cols
  knitr::kable() %>% # turn into a kable table for nice rendering in HTML
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(c("", "sms" = 2, "")) #add header for sms

If you want to learn more about the syntax of tabyl(), make sure to check the vignette

Exercise

evaluation_df %>% # your data frame
  janitor::tabyl(malaria, sms) %>% # the two dimensions for the table (Y, Z)
  janitor::adorn_totals(c("row", "col")) %>% # add totals for rows and cols
  knitr::kable() %>% # turn into a kable table for nice rendering in HTML
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(c("", "sms" = 2, "")) #add header for sms

ggplot(evaluation_df, aes(x = factor(sms), 
                          y = factor(net_use), 
                          color = factor(malaria))) +
  geom_jitter() +
  theme_minimal() +
  scale_x_discrete(labels = c("SMS = 0", "SMS = 1")) +
  scale_y_discrete(limits = c("1","0"), labels = c("NET = 1", "NET = 0")) +
  scale_color_manual(values = c("#CDCDCD", "#CC0055"),
                       labels = c("Not infected", "Infected")) +
  labs(x = "Encouragement",
       y = "Treatment",
       color = "") +
  theme(legend.position = "bottom")

Exploring our set-up

lm(net_use ~ sms, data = evaluation_df) %>%
  stargazer::stargazer(type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               net_use          
## -----------------------------------------------
## sms                          0.450***          
##                               (0.028)          
##                                                
## Constant                     0.330***          
##                               (0.020)          
##                                                
## -----------------------------------------------
## Observations                   1,000           
## R2                             0.205           
## Adjusted R2                    0.204           
## Residual Std. Error      0.444 (df = 998)      
## F Statistic          257.315*** (df = 1; 998)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

naive_model <- lm(malaria ~ net_use, data = evaluation_df)
stargazer::stargazer(naive_model, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               malaria          
## -----------------------------------------------
## net_use                      -0.615***         
##                               (0.023)          
##                                                
## Constant                     0.683***          
##                               (0.017)          
##                                                
## -----------------------------------------------
## Observations                   1,000           
## R2                             0.415           
## Adjusted R2                    0.414           
## Residual Std. Error      0.363 (df = 998)      
## F Statistic          707.002*** (df = 1; 998)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

itt_model <- lm(malaria ~ sms, data = evaluation_df)
stargazer::stargazer(itt_model, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               malaria          
## -----------------------------------------------
## sms                          -0.324***         
##                               (0.028)          
##                                                
## Constant                     0.504***          
##                               (0.020)          
##                                                
## -----------------------------------------------
## Observations                   1,000           
## R2                             0.117           
## Adjusted R2                    0.116           
## Residual Std. Error      0.446 (df = 998)      
## F Statistic          131.753*** (df = 1; 998)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Let's gather out local average treatment effect (LATE)

We have several options for this:

• Wald Estimator We can calculate this by hand. Let's try that now using the values from the table we created earlier. We can also calculate the average malaria rates among those who did and did not receive an SMS (no SMS = 0.504, yes SMS = 0.18).

Table: Observations across strata (D,Z)

Table: Count of malaria infections across strata (Y,Z)

Local Average Treatment Effect (LATE) manually

\[LATE = \frac{cov(Y_i,Z_i)}{cov(D_i,Z_i)}\]

Let's take a look at our numerator \(cov(Y_i,Z_i)\), also \(ITT\)

\(E(Y|Z = 1) = \frac{90}{(410+90)} = 0.18\)
\(E(Y|Z = 0) = \frac{252}{(252+248)} = 0.504\)

Let's take a look at our denominator \(cov(D_i,Z_i)\)

\(E(D∣Z = 1) = \frac{390}{(390 + 110)} = 0.78\)
\(E(D∣Z = 0) = \frac{165}{(165 + 335)} = 0.33\)

We can calculate our LATE

\[LATE = \frac{(0.18 - 0.504)}{(0.78 - 0.33)} = -0.72\]

• Two-stage Least Squares (2SLS). We will learn how to do this with ivreg(), which is part of the AER package. It fits instrumental-variable regression through two-stage least squares. The syntax is the following:

ivreg(outcome ~ treatment | instrument, data)

In our case:

late_model <- AER::ivreg(malaria ~ net_use | sms, data = evaluation_df)
summary(late_model)

## 
## Call:
## AER::ivreg(formula = malaria ~ net_use | sms, data = evaluation_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7416 -0.0216 -0.0216  0.2584  0.9784 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.74160    0.03089   24.00   <2e-16 ***
## net_use     -0.72000    0.05159  -13.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3671 on 998 degrees of freedom
## Multiple R-Squared: 0.4025,  Adjusted R-squared: 0.4019 
## Wald test: 194.8 on 1 and 998 DF,  p-value: < 2.2e-16

Exercise

• What is the substantive reading of these results?
• What would you tell the funding partner in your evaluation report?

Mechanics behind two-stage least squares (2SLS)

What ivreg() is doing in the background is the following:

net_use_hat <- lm(net_use ~ sms, evaluation_df)$fitted.values # get fitted values from first stage (the part of x that is exogenously driven by z)
summary(lm(evaluation_df$malaria ~ net_use_hat)) # run second stage with instrumented x (careful, the standard errors are wrong; better use ivreg() from AER instead)

## 
## Call:
## lm(formula = evaluation_df$malaria ~ net_use_hat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.504 -0.180 -0.180  0.496  0.820 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.74160    0.03757   19.74   <2e-16 ***
## net_use_hat -0.72000    0.06273  -11.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4463 on 998 degrees of freedom
## Multiple R-squared:  0.1166, Adjusted R-squared:  0.1157 
## F-statistic: 131.8 on 1 and 998 DF,  p-value: < 2.2e-16

Thinking about the validity of instruments

We can also adapt the ivreg() syntax to accomodate valid conditional instruments:

AER::ivreg(Y ~ D + W | Z + W, data = your_df)