# Book: Mastering Metrics
# Chapter 4: Regression discontinuity design
# Figure 4.5: RD estimates of MLDA effects on mortality by cause of death

data = read.table("age-mva-internal.txt",header=TRUE)

attach(data)

n = nrow(data)

over21 = rep(0,n)
over21[age>21]=1


par(mfrow=c(1,1))
plot(age,mva,pch=over21+15,ylim=c(10,40),xlab="Age",ylab="Death rate (per 100,000)")
points(age,internal,col=2,pch=over21+15)
abline(v=21,lty=2)
legend("bottomright",legend=c("Motor Vehicle Accidents","Deaths from Internal Causes"),
       col=1:2,pch=16,bty="n")
abline(v=21,lty=2)

# Linear models
fit  = lm(mva~age+over21)
right  = c(1,21,1)%*%fit$coef
left   = c(1,21,0)%*%fit$coef
delta  =  right-left
sigma = summary(fit)$sigma

fit1  = lm(internal~age+over21)
right1 = c(1,21,1)%*%fit1$coef
left1  = c(1,21,0)%*%fit1$coef
delta1 =  right1-left1
sigma1 = summary(fit1)$sigma

par(mfrow=c(1,2))
plot(age,mva,pch=over21+15,ylim=c(10,40),xlab="Age",ylab="Death rate (per 100,000)")
points(age,internal,col=2,pch=over21+15)
abline(v=21,lty=2)
lines(age,fit$fit,lwd=3)
lines(age,fit$fit,lwd=3)
lines(age[over21==0],fit$fit[over21==0],lwd=3)
lines(age[over21==1],fit$fit[over21==1],lwd=3)
lines(age[over21==0],fit1$fit[over21==0],lwd=3,col=2)
lines(age[over21==1],fit1$fit[over21==1],lwd=3,col=2)
legend("bottomright",legend=c(
       paste("MVA: ",round(delta/sigma,3)," stdevs",sep=""),
       paste("DIC: ",round(delta1/sigma1,3)," stdevs",sep="")),
       col=1:2,pch=16,bty="n")
title("Motor Vehicle Accidents (MVA)\nDeaths from Internal Causes (DIC)")


# Quadratic models
age2       = age^2
over21age  = over21*age
over21age2 = over21*age2

fit  = lm(mva~age+over21+age2+over21age+over21age2)
right  = c(1,21,1,21^2,21,21^2)%*%fit$coef
left   = c(1,21,0,21^2,0,0)%*%fit$coef
delta  =  right-left
sigma = summary(fit)$sigma

fit1  = lm(internal~age+over21+age2+over21age+over21age2)
right1 = c(1,21,1,21^2,21,21^2)%*%fit1$coef
left1  = c(1,21,0,21^2,0,0)%*%fit1$coef
delta1 =  right1-left1
sigma1 = summary(fit1)$sigma

plot(age,mva,pch=over21+15,ylim=c(10,40),xlab="Age",ylab="Death rate (per 100,000)")
points(age,internal,col=2,pch=over21+15)
abline(v=21,lty=2)
lines(age[over21==0],fit$fit[over21==0],lwd=3)
lines(age[over21==1],fit$fit[over21==1],lwd=3)
lines(age[over21==0],fit1$fit[over21==0],lwd=3,col=2)
lines(age[over21==1],fit1$fit[over21==1],lwd=3,col=2)
legend("bottomright",legend=c(
       paste("MVA: ",round(delta/sigma,3)," stdevs",sep=""),
       paste("DIC: ",round(delta1/sigma1,3)," stdevs",sep="")),
       col=1:2,pch=16,bty="n")
title("Motor Vehicle Accidents (MVA)\nDeaths from Internal Causes (DIC)")