#
# plotWaVBMEffect.r
#
# Plot county-level turnout across time in Washington
# State, noting the year that each county switched to 
# all-mail elections.
#
# Seth Hill, March 2011
#

rm(list = ls())
if (.Platform$OS.type == "windows") { # Set the working directory to the location of the script.
  frame_files <- lapply(sys.frames(), function(x) x$ofile)
  frame_files <- Filter(Negate(is.null), frame_files)
  PATH <- dirname(frame_files[[length(frame_files)]])
  setwd(PATH); rm(PATH,frame_files)
}

x.election.years <- T # Set to TRUE to have x-axis in election years, rather than in time-years, i.e.
                      # "elections before after switch" versus "years before/after switch."
plot.years <- seq(1996,2010,2)
title.figures <- FALSE # Toggle on/off printing figure titles.

# Data downloaded from http://www.sos.wa.gov/elections/absentee_stats.aspx?y=2000.
dat <- read.csv("WA-County-VotesCast.csv",as.is=T)

# Identify year of VBM switch.
dat <- dat[order(dat$County,dat$Year),]
dat[,'newlyVBM'] <- as.numeric(c(0,diff(dat$VBMOnly)) == 1)
cat("Year of VBM Switch By County:\n")
yearOfSwitch <- dat[dat$newlyVBM == 1,c('County','Year')]
names(yearOfSwitch) <- c("County","YearOfSwitch")
print(yearOfSwitch)
cat("Before merge:",dim(dat),"\n")
dat <- merge(dat,yearOfSwitch,by="County",all.x=T,all.y=F)
cat("After merge:",dim(dat),"\n")
dat[,'yearsSinceSwitch'] <- with(dat,Year - YearOfSwitch)
#dat[is.na(dat$yearsSinceSwitch),'yearsSinceSwitch'] <- 0
#dat[dat$yearsSinceSwitch < 0,'yearsSinceSwitch'] <- 0
dat <- dat[order(dat$County,dat$Year),]
# Note: There are a few discrepancies between these years and those documented on 
# the pdf at http://wei.secstate.wa.gov/osos/en/ElectionStatistics/VBM/Documents/
# 2000-2010%20WA%20history%20of%20move%20to%20vbm.pdf. My readng of the turnout 
# numbers indicate this is accurate. Any inaccuracies surround 2005 versus 2006,
# so should not much affect things.

# Line plots by county, separately for Even-Years and Odd-Years.
plots <- list("Odd-Year"=list(name="Odd-Year",flag=(!dat$Year %in% plot.years & dat$Year %% 2 == 1)),
              "Even-Year"=list(name="Even-Year",flag=dat$Year %in% plot.years))
              
rddLm <- function(dat,y,cutoff,forcingvar,power=1,covars=NULL,
                  interaction=NULL,...) {
  # Calculates the regression discontinuity effect via a linear model.
  #
  # Args:
  #  dat: Data frame or named data matrix holding variables.
  #  y: Response variable name.
  #  cutoff: Indicator variable for before/after cutoff.
  #  forcingvar: Variable that generates the cutoff.
  #  power: Polynomial power for forcing variable. Default is linear.
  #  covars: Variable names of covariates in model; Default NULL is
  #          no covars.
  #  interaction: Variable to interact with cutoff to modify treatment
  #          effect. Default NULL is no interaction.
  #  ...: Additional arguments to pass to lm.
  #
  # Returns list with elements:
  #  itt: Intent to treat effect (coefficient on cutoff).
  #  itt.se: Intent to treat standard error (se on cutoff).
  #  lm.res: lm result object.
  #

  # Build up model formula.
  form <- sprintf("%s ~ %s*%s",y,cutoff,forcingvar)
  if (power > 1) { # Build up polynomial of forcing variable.
    for (pow in 2:power) {
      form <- sprintf("%s + %s*I(%s^%s)",form,cutoff,forcingvar,pow)
    }
  }
  if (!is.null(interaction)) {
    form <- sprintf("%s + %s:%s",form,cutoff,interaction)
  }
  if (!is.null(covars)) {
    for (cv in covars) {
      # Should the discontinuity be interacted with covariates?
      #form <- sprintf("%s + %s*%s)",form,cutoff,cv)
      form <- sprintf("%s + %s",form,cv)
    }
  }

  # Model
  lm.out <- lm(as.formula(form),data=dat,...)
  mycoefs <- summary(lm.out)$coefficients

  return(list("itt" = mycoefs[cutoff,'Estimate'],
       "itt.se" = mycoefs[cutoff,'Std. Error'],
       "lm.res" = lm.out))
}

#
# Make plots.
#
for (p in plots) {
  cat(p$name,"\n")
  d2 <- dat[p$flag,]
  d2 <- d2[order(d2$County,d2$Year),]
  d2[,'firstVBM'] <- as.numeric(c(0,diff(d2$VBMOnly)) == 1)
  
  # Grab residuals from regressing turnout on county and year fixed effects.
  lm.out <- lm(TurnoutPct ~ as.factor(Year) + as.factor(County), data=d2)
  d2[,'resid'] <- residuals(lm.out)
  
  if (x.election.years) {
    # Recode yearsSinceSwitch variable to actually measure "elections of
    # type since switch."
    d2 <- d2[order(d2$County,d2$Year),]
    elecsSinceSwitch <- function(vbm.vector) {
      # Calculate the distance in units between each element of a vector
      # and the single element of that vector equal to 1.
      if (any(is.na(vbm.vector))) return(rep(NA,length(vbm.vector)))
      if (!any(vbm.vector == 1))  return(rep(NA,length(vbm.vector)))
      ind.one <- which(vbm.vector == 1)
      return( 1:length(vbm.vector) - ind.one )
    }
    d2[,'ess'] <- unlist(tapply(d2$firstVBM,d2$County,elecsSinceSwitch))
    d2[,'yearsSinceSwitch'] <- d2$ess
  }

  # Plot of residuals on years before/after switch, with overlying smooth.
  pdf(sprintf("WAVBMEffect-%s.pdf",p$name),width=10)
  par(mar=c(5.1,4.1,ifelse(title.figures,2.1,.6),1.1))
  with(d2,plot(x=yearsSinceSwitch,y=resid,type='n',ann=F,axes=F,ylim=c(-.1,.1)))
  abline(v=-.5,lwd=2,lty=2)
  text(x=-.5,y=.1,"<-- Polling Place Elections",pos=2)
  text(x=-.5,y=.1,"All-Mail Elections -->",pos=4)
  axis(1,at=seq(-14,8,2)) ; axis(2,las=2)
  title(xlab=paste(p$name,ifelse(x.election.years,'Elections','Years'),
                  'Before/After Switch to All-Mail -- Lines Connect Same County.'),
        ylab='Residual From Year and County Fixed Effects Regression Model')
  if (title.figures) {
    title(main=sprintf("Effect of Vote By Mail on\n%s Election Turnout in Washington Counties",p$name))
  }
  for (county in unique(d2$County)) {
    flag <- d2$County == county
    lines(x=d2[flag,'yearsSinceSwitch'],y=d2[flag,'resid'],col='lightgray',lwd=1)
  }
  with(d2,points(x=yearsSinceSwitch,y=resid,pch=19,cex=.7))
  
  # Add smooth lines
  d3 <- na.omit(d2[,c('yearsSinceSwitch','resid')])
  library(RColorBrewer)
  mycols <- brewer.pal(6,"Set1")[c(3,4)]

  # RDD estimate of the discontinuity effect.
  d3$cutoff <- as.numeric(d3$yearsSinceSwitch >= 0)
  rdd.res <- rddLm(dat=d3,y="resid",cutoff="cutoff",forcingvar="yearsSinceSwitch",
                  power=2,covars=NULL,interaction=NULL)
  newdata <- data.frame(yearsSinceSwitch=sort(unique(d3[,'yearsSinceSwitch'])))
  newdata$cutoff <- as.numeric(newdata$yearsSinceSwitch >= 0)
  pred <- predict(rdd.res$lm.res,newdata,se.fit=T,interval="confidence")$fit
  cat(sprintf(" Predicted effect of discontinuity: %1.3f (%1.3f).\n",
    100*rdd.res$itt,100*rdd.res$itt.se))
  # Plot the effect for the cubic polynomial.
  flag.cutt <- newdata$cutoff == 1
  lines(x=newdata[!flag.cutt,1],y=pred[!flag.cutt,1],col=mycols[1],lwd=4)
  lines(x=newdata[!flag.cutt,1],y=pred[!flag.cutt,2],col=mycols[1],lwd=2,lty=2)
  lines(x=newdata[!flag.cutt,1],y=pred[!flag.cutt,3],col=mycols[1],lwd=2,lty=2)
  lines(x=newdata[flag.cutt,1],y=pred[flag.cutt,1],col=mycols[2],lwd=4)
  lines(x=newdata[flag.cutt,1],y=pred[flag.cutt,2],col=mycols[2],lwd=2,lty=2)
  lines(x=newdata[flag.cutt,1],y=pred[flag.cutt,3],col=mycols[2],lwd=2,lty=2)
  title(sub="Outer Lines are 95 Percent Confidence Interval for Fit",font.sub=3)
  legend('bottomright',col=c(mycols),lty=c(1,1),lwd=c(3,3),bty='n',
    legend=c('RDD Fit Before Change','RDD Fit After Change'))
  dev.off()
  
  # Plot residuals against year of election, with different pch for all-mail and 
  # polling place elections.
  with(d2,plot(x=Year,y=resid,pch=ifelse(VBMOnly==1,1,2),col=ifelse(VBMOnly==1,mycols[2],mycols[1])))
}

# Summary of turnout variability.
cat("Maximum and minimum even-year turnout:\n")
print(dat[plots[['Even-Year']][['flag']],][which.max(dat[plots[['Odd-Year']][['flag']],'TurnoutPct']),c('County','Year','TurnoutPct')])
print(dat[plots[['Even-Year']][['flag']],][which.min(dat[plots[['Odd-Year']][['flag']],'TurnoutPct']),c('County','Year','TurnoutPct')])

cat("Maximum and minimum odd-year turnout:\n")
print(dat[plots[['Odd-Year']][['flag']],][which.max(dat[plots[['Odd-Year']][['flag']],'TurnoutPct']),c('County','Year','TurnoutPct')])
print(dat[plots[['Odd-Year']][['flag']],][which.min(dat[plots[['Odd-Year']][['flag']],'TurnoutPct']),c('County','Year','TurnoutPct')])
