HW06 - Regression Discontinuity

Carpenter & Dobkin (2009) data

library(haven)
library(tidyverse)
library(dplyr)
library(stargazer)

#Loading Data
AEJfigs <- read_dta("D:/Downloads/AEJfigs.dta")

# Creating required variable
over21 = as.factor(ifelse(AEJfigs$agecell >= 21, 1, 0))

AEJfigs <- AEJfigs %>%
    mutate(over21)
AEJfigs_2 <- AEJfigs[complete.cases(AEJfigs),]

Replicating Figure 3

plot(AEJfigs$agecell,AEJfigs$all,
     pch=18,
     xlim=c(19, 23),ylim = c(0,120),
     xlab="Age",ylab="Deaths per 100,000 person-years",
     main="Figure 3. Age Profile for Death Rates")

points(AEJfigs$agecell,AEJfigs$external,pch=17,cex = 0.8)
points(AEJfigs$agecell,AEJfigs$internal,pch=22,cex = 0.8)
points(AEJfigs$agecell,AEJfigs$all, pch= 16, col="grey") # plotting data points

L21 <- nrow(AEJfigs[AEJfigs$over21== 0,]) 
lines(AEJfigs$agecell[1:L21], AEJfigs$allfitted [1:L21], lty=3, col = "black") # Plotting fitted lines data points
lines(AEJfigs$agecell[(L21 + 1): nrow(AEJfigs)],
      AEJfigs$allfitted [(L21 + 1):nrow(AEJfigs)], lty=3, col = "black")

lines(AEJfigs$agecell[1:L21], AEJfigs$externalfitted[1:L21], lty=1)
lines(AEJfigs$agecell[(L21 + 1): nrow(AEJfigs)], AEJfigs$externalfitted[(L21 + 1): nrow(AEJfigs)], lty=1)

lines(AEJfigs$agecell[1:L21], AEJfigs$internalfitted [1:L21], lty=5)
lines(AEJfigs$agecell[(L21 + 1): nrow(AEJfigs)], AEJfigs$internalfitted [(L21 + 1): nrow(AEJfigs)], lty=5)
     
legend("right",
       lty=c(0,0,0,3,5,1),
       pch=c(16,22,17,NA,NA,NA),
       col = c("grey","black","black","black","black","black"),
       legend = c("All","Internal","External",
                  "All fitted", "Internal fitted", "External fitted"),
       ncol=2, cex = 0.4) # Editing legend

Simplest regression-discontinuity design

simple_RDD <- lm(all ~ agecell + over21, data=AEJfigs_2) # running regression

stargazer(simple_RDD, type = "text", title = "Simplest regression-discontinuity design")

## 
## Simplest regression-discontinuity design
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 all            
## -----------------------------------------------
## agecell                       -0.975           
##                               (0.632)          
##                                                
## over211                      7.663***          
##                               (1.440)          
##                                                
## Constant                    112.310***         
##                              (12.668)          
##                                                
## -----------------------------------------------
## Observations                    48             
## R2                             0.595           
## Adjusted R2                    0.577           
## Residual Std. Error       2.493 (df = 45)      
## F Statistic           32.995*** (df = 2; 45)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Plotting fitted results

fitted_RDD <- fitted(simple_RDD) # fitting the regression line

AEJfigs_2 <- AEJfigs_2 %>%
  mutate(fitted_RDD)                # Merging the var

#Plotting the graph

plot(AEJfigs_2$agecell,AEJfigs_2$all,
     pch=18,
     xlim=c(19, 23),ylim = c(85,110),
     xlab="Age",ylab="Deaths per 100,000 person-years",
     main="Figure: Simple RDD")
lines(AEJfigs_2$agecell[1:L21], fitted_RDD[1:L21], lty=1)
lines(AEJfigs_2$agecell[(L21 + 1): nrow(AEJfigs_2)],
      fitted_RDD[(L21 + 1): nrow(AEJfigs_2)], lty=1)
legend("topleft", 
       lty=c(0,1),
       pch=c(18,NA),
       legend = c("All"
                  ,"Fitted"),
       ncol=2, cex = 0.6)