Month: April 2019

library(ggplot2)
library(dplyr)
library(ggthemes)
library(scales)
library(colorspace)
library(tidyr)
library(sf)


setwd("C:/Users/Jonathan Tannen/Dropbox/sixty_six/posts/turnout_2019_primary/")

load("../../data/processed_data/df_major_2017_12_01.Rda")

source("../../admin_scripts/util.R")

turnout <- df_major %>%
  filter(OFFICE %in% c(
    "PRESIDENT OF THE UNITED STATES",
    "GOVERNOR",
    "MAYOR",
    "DISTRICT ATTORNEY"
  )) %>%
  group_by(WARD16, DIV16, year, election) %>%
  summarise(VOTES = sum(VOTES))


df_2018 <- read.csv("../../data/raw_election_data/2018_general.csv")
names(df_2018)<- c(
  "WARD16", "DIV16", "TYPE", "OFFICE", 
  "CANDIDATE", "PARTY", "VOTES"
)
df_2018$WARD16 <- sprintf("%02d", df_2018$WARD16)
df_2018$DIV16 <- sprintf("%02d", df_2018$DIV16)

df_2018 <- df_2018 %>% 
  filter(OFFICE == "GOVERNOR AND LIEUTENANT GOVERNOR") %>%
  group_by(WARD16, DIV16) %>%
  summarise(VOTES = sum(VOTES)) 

df_2018$election <- "general"
df_2018$year <- "2018"

df_2018_primary <- read.csv("../../data/raw_election_data/2018_primary.csv")
# head(df_2018_primary)
names(df_2018_primary)<- c(
  "WARD16", "DIV16", "TYPE", "OFFICE", 
  "CANDIDATE", "PARTY", "VOTES"
)
df_2018_primary$WARD16 <- sprintf("%02d", df_2018_primary$WARD16)
df_2018_primary$DIV16 <- sprintf("%02d", df_2018_primary$DIV16)
df_2018_primary <- df_2018_primary %>% 
  filter(PARTY == "DEMOCRATIC") %>%
  mutate(OFFICE = gsub("(.*)-DEM", "\\1", OFFICE)) df_2018_primary <- df_2018_primary %>% filter(OFFICE == "GOVERNOR") %>% group_by(WARD16, DIV16) %>% summarise(VOTES = sum(VOTES)) df_2018_primary$election <- "primary" df_2018_primary$year <- "2018" turnout <- bind_rows(turnout, df_2018) turnout <- bind_rows(turnout, df_2018_primary) turnout_wide <- turnout %>% unite(key, election, year) %>% spread(key = key, value = VOTES) cycles <- data.frame( year = 2002:2021, cycle = rep(c("Governor","Mayor","President","District Attorney"), 5), senate = rep(c(FALSE, FALSE, TRUE, FALSE, TRUE, FALSE), 4)[1:20] ) turnout_total <- turnout %>% group_by(year, election) %>% summarise(VOTES = sum(VOTES)) turnout_total <- turnout_total %>% left_join( cycles %>% mutate(year = as.character(year)), by = "year" ) 

annotation_df <- list(
  primary = tribble(
    ~year, ~VOTES, ~cycle, ~hjust,
    13, 75e3, "District Attorney", 0,
    14, 375e3, "President", 1,
    14, 180e3, "Governor", 0,
    14, 270e3, "Mayor", 0.5
  ),
  general = tribble(
    ~year, ~VOTES, ~cycle, ~hjust,
    13.5, 120e3, "District Attorney", 0,
    12, 675e3, "President", 0,
    14.2, 400e3, "Governor", 0,
    14.3, 248e3, "Mayor", 0
  )
)


senate_label_pos <- list(
  primary = tribble(
    ~year, ~VOTES, ~senate, ~label,
    2.3, 25e3, TRUE, "Senate",
    2.3, 50e3, FALSE, "Non-Senate"
  ),
  general = tribble(
    ~year, ~VOTES, ~senate, ~label,
    2.3, 37.5e3, TRUE, "Senate",
    2.3, 75e3, FALSE, "Non-Senate"
  )
)


turnout_plot <- function(use_election){
  ggplot(
    turnout_total %>% filter(election == use_election), 
    aes(
      x = year, 
      y = VOTES,
      color = cycle,
      group = interaction(cycle, election)
    )
  ) + 
    geom_point(size = 3, aes(shape = senate)) +
    geom_line() +
    expand_limits(y = 0) +
    scale_y_continuous("Votes Cast", labels = comma) +
    theme_sixtysix() +
    theme(axis.title.x = element_blank()) +
    geom_point(
      data = senate_label_pos[[use_election]],
      x = 2,
      aes(shape = senate),
      color = "grey20",
      group = NA,
      size = 3
    )+
    geom_text(
      data = senate_label_pos[[use_election]],
      aes(label=label),
      color = "grey20",
      group = NA,
      size = 4,
      hjust = 0
    )+
    geom_text(
      data = annotation_df[[use_election]],
      aes(label=cycle, hjust=hjust, color=cycle),
      group = NA,
      size = 4,
      fontface="bold"
    )+
    scale_shape_discrete(guide = FALSE)+
    scale_color_discrete(guide = FALSE)+
    ggtitle(sprintf(
      "Turnout in Philadelphia %s",
      ifelse(use_election == "general", "Generals", "Democratic Primaries")
    ))
}

turnout_plot("primary")

turnout_plot("general")

da_results <- df_major %>%
  filter(
    election == "primary" & PARTY == "DEMOCRATIC" & OFFICE == "DISTRICT ATTORNEY"
  ) %>%
  group_by(year, WARD16, DIV16, CANDIDATE) %>%
  summarise(VOTES = sum(VOTES)) %>%
  group_by(year, WARD16, DIV16) %>%
  mutate(
    total_votes = sum(VOTES),
    pct_vote = VOTES / total_votes
  )

da_results$candidate_name <- format_name(da_results$CANDIDATE)

turnout_wide <- turnout_wide %>%
  mutate(
    typical_turnout_da = (primary_2013 + primary_2009 + primary_2005)/3,
    typical_turnout_governor = (general_2014 + general_2010 + general_2006 + general_2002)/4
  )

krasner_results <- da_results %>% 
  filter(candidate_name == "Lawrence S Krasner") %>%
  left_join(turnout_wide)

library(sf)

divs <- st_read("../../data/gis/2016/2016_Ward_Divisions.shp", quiet = TRUE)
divs <- divs %>% st_transform(2272)
wards <- st_read("../../data/gis/2016/2016_Wards.shp",  quiet = TRUE)
wards <- wards %>% st_transform(2272)

divs$area <- as.numeric(st_area(divs$geometry)) / (5280^2)
wards$area <- wards$AREA_SFT / (5280^2)

divs <- st_simplify(divs, 500)
divs <- divs %>% mutate(
  WARD16 = sprintf("%02d", WARD),
  DIV16 = sprintf("%02d", DIVSN)
) %>% select(-WARD, -DIVSN)
wards$WARD16 = sprintf("%02d", asnum(wards$WARD))

krasner_results_wards <- krasner_results %>%
  group_by(WARD16) %>%
  summarise(
    turnout_2017 = sum(primary_2017),
    typical_turnout_da = sum(typical_turnout_da),
    pct_vote = weighted.mean(pct_vote, w=total_votes)
  )

turnout_wide_wards <- turnout_wide %>%
  group_by(WARD16) %>%
  summarise_at(
    .funs = funs(sum(., na.rm = TRUE)), 
    vars(
      starts_with("primary_"), 
      starts_with("general_"), 
      starts_with("typical_turnout_")
    )
  )

krasner_turnout_per_mile <- ggplot(
  # divs %>%
  #   left_join(turnout_wide)
  wards %>% 
    left_join(turnout_wide_wards)
) +
  geom_sf(
    aes(
      fill = pmin(primary_2017 / area, 10e3)
    ), 
    color = NA
  ) +
  scale_fill_viridis_c(
    "Votes per mile",
    labels = scales::comma
  ) +
  theme_map_sixtysix() +
  ggtitle("Votes per mile in the 2017 primary")

krasner_turnout_change <- ggplot(
  # divs %>%
  #   left_join(turnout_wide)
  wards %>%
    left_join(turnout_wide_wards)
) +
  geom_sf(
    aes(fill = pmin(primary_2017 / typical_turnout_da, 4)), 
    color = NA
  ) +
  expand_limits(fill = c(1,3)) +
  scale_fill_viridis_c(
    "Turnout in 2017/\n Typical DA Turnout",
    # breaks = 0:5,
    labels = function(x) paste0(x, "x")
  ) +
  theme_map_sixtysix() +
  ggtitle("Surged nearly 3x in Krasner's base")

gridExtra::grid.arrange(
  krasner_turnout_per_mile,
  krasner_turnout_change, 
  nrow=1  
)

votes_per_mile_2018 <- ggplot(
  # divs %>%
    # left_join(turnout_wide)
    wards %>% left_join(turnout_wide_wards)
) +
  geom_sf(
    aes(
      fill = pmin(general_2018 / area, 20e3)
    ), 
    color = NA
  ) +
  scale_fill_viridis_c(
    "Votes per mile",
    labels = scales::comma
  ) +
  theme_map_sixtysix() +
  ggtitle("Votes per mile in the 2018 general")

turnout_change_2018 <- ggplot(
  # divs %>%
  #   left_join(turnout_wide)
  wards %>%
    left_join(turnout_wide_wards)
) +
  geom_sf(
    aes(fill = pmin(general_2018 / typical_turnout_governor, 4)), 
    color = NA
  ) +
  expand_limits(fill = c(1, 3)) +
  scale_fill_viridis_c(
    "Turnout in 2018/\n Typical Governor Turnout",
    labels = function(x) paste0(x, "x")
  ) +
  theme_map_sixtysix() +
  ggtitle(
    'Surged "only" 2x in Krasner\'s base'
  )

gridExtra::grid.arrange(
  votes_per_mile_2018,
  turnout_change_2018, 
  nrow=1  
)

krasner_pct_map <- ggplot(
  divs %>% left_join(krasner_results)
  # wards %>% 
  #   left_join(krasner_results_wards)
) +
  geom_sf(aes(fill = 100 * pct_vote), color = NA) +
  scale_fill_viridis_c("Percent\n of vote") +
  theme_map_sixtysix() +
  ggtitle("Percent of vote for Krasner")

print(krasner_pct_map)

ggplot(
  krasner_results,
  aes(
    y = primary_2017 / typical_turnout_da,
    x = 100 *pct_vote
  )
) +
  geom_point(
    aes(
      size = typical_turnout_da
    ),
    alpha = 0.3,
    color = strong_green,
    pch = 16
  ) +
  geom_smooth(
    aes(weight = typical_turnout_da), 
    color = "grey10"
  )+
  scale_size_area("Division's average votes\nin 2005, 2009, 2013") +
  scale_y_continuous(
    labels = function(x) return(paste0(x, "x")),
    limits = c(0, 8),
    breaks = seq(0, 10, 2)
  ) +
  geom_hline(yintercept = 1, linetype="dashed") +
  annotate(
    "text", 
    x=55, y=0.95,
    label="2017 turnout = typical turnout",
    hjust=0, vjust=1
  ) +
  labs(
    title = "Krasner's popularity also drove turnout",
    subtitle = "Democratic Primary turnout was 1.5x in low-support divisions, but 3x in high",
    y = "Votes in 2017 / Average votes in 2005, 2009, 2013",
    x = "Percent of vote for Krasner"
  )+
  theme_sixtysix()

load_council_races <- function(year){
  df_year <- read.csv(
    paste0("../../data/raw_election_data/",year,"_primary.csv")
  )

  df_year <- df_year %>%
    mutate(
      WARD16 = sprintf("%02d", asnum(WARD)),
      DIV16 = sprintf("%02d", asnum(DIVISION))
    ) %>%
    group_by(WARD16, DIV16, OFFICE, CANDIDATE, PARTY) %>%
    summarise(VOTES = sum(VOTES)) %>%
    group_by()


  district_regex <- "DISTRICT COUNCIL(-|\\s)([0-9]+)[A-Z]+ DIST(RICT)?-D(EM)?"
  council_districts <- df_year %>%
    filter(grepl(district_regex, OFFICE)) %>%
    mutate(
      district = asnum(gsub(district_regex, "\\2", OFFICE))
    ) %>%
    mutate(
      candidate_name = format_name(CANDIDATE),
      last_name = get_last_name(candidate_name),
      year = year
    )
  return(council_districts)
}

council_2015 <- load_council_races(2015)
council_2011 <- load_council_races(2011)
council_2007 <- load_council_races(2007)
council_2003 <- load_council_races(2003)

council_df <- bind_rows(
  council_2015, council_2011, council_2007, council_2003
)

council_totals <- council_df %>%
  group_by(year, candidate_name, district) %>%
  summarise(votes = sum(VOTES)) %>%
  arrange(year, district, desc(votes))

council_races <- council_totals %>%
  group_by(year, district) %>%
  summarise(
    winner = candidate_name[which.max(votes)],
    pct_winner = max(votes) / sum(votes),
    is_competitive = pct_winner < 0.9
  ) %>%
  group_by()

council_turnout <- turnout %>% 
  rename(total_votes = VOTES) %>%
  group_by() %>%
  filter(
    election == "primary" &
      year %in% seq(2003, 2015, 4)
  ) %>%
  mutate(year = asnum(year)) %>%
  left_join(
    council_df %>%
      select(WARD16, DIV16, district) %>%
      unique()
  ) %>%
  left_join(council_races)

## FYI: 2015 did not have a Dem Primary in the 10th

council_turnout$is_competitive <- with(
  council_turnout,
  replace(is_competitive, is.na(is_competitive), FALSE)
)

council_turnout <- council_turnout %>%
  mutate(WARD_DIVSN = paste0(WARD16, DIV16))

council_turnout <- council_turnout %>%
  mutate(competitive_mayor = year %in% c(2015, 2007))


fit_competitive <- lm(
  log(total_votes + 1) ~ 
    as.character(year) +
    WARD_DIVSN +
    is_competitive * competitive_mayor,
  data = council_turnout
) 

coef_council_is_competitiveTRUE <- coef(fit_competitive)['is_competitiveTRUE']

coef_council_mayor_is_competitive_interaction <- coef(fit_competitive)['is_competitiveTRUE:competitive_mayorTRUE']

# ncoef <- length(fit_competitive$coefficients)
# fit_competitive$coefficients %>% tail(4)
# vcov <-vcov(fit_competitive)
# vcov[(nrow(vcov)-2):nrow(vcov), (nrow(vcov)-2):nrow(vcov)] %>% diag %>% sqrt

projected_turnout <- turnout_wide %>%
  left_join(
    council_df %>%
      filter(year == 2015) %>%
      select(WARD16, DIV16, district) %>%
      unique()
  ) 

projected_turnout <- projected_turnout %>%
  mutate(
    is_contested_2019 = district %in% c(1,2,3,4,7)
  ) %>%
  left_join(
    council_races %>%
      select(year, district, is_competitive) %>%
      filter(year %in% c(2011, 2003)) %>%
      mutate(year = paste0("was_contested_", year)) %>%
      spread(year, is_competitive, fill=FALSE)
  ) 

projected_turnout <- projected_turnout %>%
  mutate(
    baseline_turnout_avg = 0.5 * (
      primary_2003 / ifelse(
        was_contested_2003, exp(coef_council_is_competitiveTRUE), 1 
      ) + 
        primary_2011 / ifelse(
          was_contested_2011, exp(coef_council_is_competitiveTRUE), 1
        )
    ),
    baseline_turnout_2011 =
      primary_2011 / ifelse(
        was_contested_2011, exp(coef_council_is_competitiveTRUE), 1
      )
  )

projected_turnout <- projected_turnout %>%
  mutate(
    competitive_scaling = ifelse(
      is_contested_2019, exp(coef_council_is_competitiveTRUE), 1
    )
  )

projected_turnout <- projected_turnout %>%
  left_join(
    turnout_wide %>% 
      select(
        WARD16, DIV16, primary_2017, 
        general_2018, typical_turnout_da, typical_turnout_governor
      ) %>%
      mutate(
        scale_2017 = primary_2017 / typical_turnout_da,
        scale_2018 = general_2018 / typical_turnout_governor
      )
  ) %>%
  mutate(
    high_projection = baseline_turnout_2011 *
      competitive_scaling *
      # (scale_2018 + scale_2017)/2
      scale_2018,
    low_projection = baseline_turnout_avg *
      competitive_scaling
  )

baseline_turnout <- sum(
  projected_turnout$low_projection,
  na.rm = TRUE
)
surged_turnout <- sum(
  projected_turnout$high_projection, 
  na.rm = TRUE
)

turnout_2019 <- tribble(
  ~year, ~election, ~cycle, ~senate,
  '2019', "primary", "Mayor", FALSE
) %>% full_join(
  data.frame(
    year = '2019',
    sim = c("Low", "High"),
    VOTES = c(baseline_turnout, surged_turnout)
  )
)

turnout_plot("primary") +
  geom_point(
    data=turnout_2019,
    size=3
  ) +
  geom_segment(
    data=turnout_2019,
    aes(
      xend=year,
      yend=VOTES,
      color=cycle
    ),
    x=14,
    y=turnout_total %>% 
      filter(year == 2015 & election == "primary") %>% 
      with(VOTES),
    linetype="dashed"
  ) +
  geom_text(
    data=turnout_2019,
    aes(label = sim),
    vjust = -1
  ) +
  labs(subtitle = "Projections under pre- and post-2016 assumptions.")

write.csv(
  projected_turnout %>% 
    select(
      WARD16, DIV16, district,
      high_projection, low_projection, 
      primary_2011, primary_2015, 
      typical_turnout_governor, general_2018, 
      typical_turnout_da, primary_2017
    ),
  file = "turnout_projections_2019.csv",
  row.names=FALSE
)

library(tidyverse)
library(sp)
library(rgeos)
library(rgdal)
library(sf)
source("../../admin_scripts/util.R")

df_major <- safe_load("../../data/processed_data/df_major_2017_12_01.Rda")

df_major$WARD_DIVSN <- with(df_major, paste0(WARD16, DIV16))

df_major <- df_major %>%
  filter(
    election == "primary" & CANDIDATE != "Write In" & PARTY == "DEMOCRATIC"
  )
df_major <- df_major %>%
  group_by(WARD_DIVSN, OFFICE, year) %>%
  mutate(pct_vote = VOTES / sum(VOTES))

df_major <- df_major %>%
  filter(OFFICE %in% c("COUNCIL AT LARGE", "DISTRICT ATTORNEY"))


bg_17_acs <- read.csv("../../data/census/acs_2013_2017_phila_bg_race_income.csv")
bg_17_acs <- bg_17_acs %>% 
  mutate(Geo_FIPS = as.character(Geo_FIPS)) %>%
  select(
    Geo_FIPS, pop, pop_nh_white, pop_nh_black, pop_nh_asian, pop_hisp, pop_median_income_2017
  )

sp_divs <- readOGR("../../data/gis/2016/2016_Ward_Divisions.shp", verbose = FALSE)
sp_divs <- spChFIDs(sp_divs, as.character(sp_divs$WARD_DIVSN))
sp_divs <- spTransform(sp_divs, CRS("+init=EPSG:4326"))

library(tigris)
options(tigris_use_cache = TRUE)
bg_shp <- block_groups(42, 101, year = 2015)
bg_shp <- spChFIDs(bg_shp, as.character(bg_shp$GEOID))
bg_shp <- spTransform(bg_shp, CRS(proj4string(sp_divs)))

sp_divs$bg <- over(
  gCentroid(sp_divs, byid = TRUE), 
  bg_shp
)$GEOID

sp_divs@data <- sp_divs@data %>%
  left_join(bg_17_acs, by = c("bg"="Geo_FIPS"))

df_major <- df_major %>%
  left_join(sp_divs@data) %>%
  mutate(
    pct_wht = pop_nh_white / pop,
    pct_blk = pop_nh_black / pop,
    pct_asian = pop_nh_asian/ pop,
    pct_hisp = pop_hisp / pop
  )

endorsements <- read_csv("da_2017_endorsements.csv")
endorsements$ward <- sprintf("%02d", endorsements$ward)

da_results <- df_major %>% 
  filter(election == "primary" & year == 2017 & OFFICE == "DISTRICT ATTORNEY") %>%
  mutate(
    last_name = gsub(
      "^.*\\s([A-Z])([A-Z]+)

quot;, "\\U\\1\\L\\2", CANDIDATE, perl = TRUE ) ) %>% group_by(WARD_DIVSN) %>% mutate(total_votes = sum(VOTES)) %>% group_by() %>% mutate(pvote = VOTES / total_votes) da_results$last_name <- with( da_results, ifelse( last_name == "Neill", "O'Neill", ifelse(last_name == "Shabazz", "El-Shabazz", last_name) ) ) da_results %>% group_by(WARD16, last_name) %>% summarise(votes = sum(VOTES)) %>% group_by(WARD16) %>% mutate( ward_votes = sum(votes), pvote = votes/ ward_votes ) %>% left_join( endorsements %>% mutate(is_endorsed = TRUE), by = c("WARD16" = "ward", "last_name" = "endorsement") ) %>% mutate( is_endorsed = replace(is_endorsed, is.na(is_endorsed), FALSE) ) %>% group_by(last_name, is_endorsed) %>% summarise( pct_vote = 100 * weighted.mean(pvote, w = ward_votes), total_votes = sum(ward_votes), n_wards = n() ) %>% group_by(last_name) %>% summarise( pct_vote_overall = weighted.mean(pct_vote, w = total_votes), wards_endorsed = ifelse(any(is_endorsed), n_wards[is_endorsed], 0), turnout_endorsed = ifelse(any(is_endorsed), total_votes[is_endorsed], 0), pct_vote_notendorsed = pct_vote[!is_endorsed], pct_vote_endorsed = ifelse(any(is_endorsed), pct_vote[is_endorsed], NA) ) %>% arrange(desc(pct_vote_overall)) %>% knitr::kable( digits = 0, format = "html", format.args = list(big.mark = ","), col.names = c( "Candidate", "Citywide % of vote", "Number of ward endorsements", "Turnout in endorsed wards", "% of vote in un-endorsed wards", "% of vote in endorsed wards" ) ) 

wards <- readOGR("../../data/gis/2016","2016_Wards", verbose=FALSE) %>%
  spTransform(CRS(proj4string(sp_divs))) 

ggwards <- fortify(spChFIDs(wards, sprintf("%02d", wards$WARD)))

bbox <- sp_divs[substr(sp_divs$WARD_DIVSN,1,2) %in% c("10", "50", "09"),] %>% 
  gUnionCascaded() %>% 
  bbox()


bbox <- rowMeans(bbox) + 1.2 * sweep(bbox, 1, rowMeans(bbox))

polygon_in_bbox <- function(p) {
  coords <- p@Polygons[[1]]@coords
  any(
    coords[,1] > bbox[1,1] &
      coords[,1] < bbox[1,2] &
      coords[,2] > bbox[2,1] &
      coords[,2] < bbox[2,2] 
  )
}
sp_divs$in_bbox <- sapply(sp_divs@polygons, polygon_in_bbox) 

ggdivs <- fortify(spChFIDs(sp_divs, as.character(sp_divs$WARD_DIVSN)))

ggdivs <- ggdivs %>%
  left_join(
    sp_divs@data %>% select(WARD_DIVSN, in_bbox),
    by = c("id" = "WARD_DIVSN")
  ) %>%
  left_join(
    da_results %>% filter(last_name %in% c("Khan", "Krasner", "El-Shabazz")),
    by = c("id" = "WARD_DIVSN")
  )

ward_centroids <- gCentroid(wards, byid=TRUE) %>% as.data.frame()
ward_centroids$ward <- wards$WARD

ggplot(
  ggdivs %>% filter(in_bbox),
  aes(x=long, y=lat)
) +
  geom_polygon(aes(fill = 100 * pvote, group=group), color = NA) +
  geom_polygon(data = ggwards, aes(group=group), fill = NA, color = "white") +
  geom_text(data = ward_centroids, aes(x=x, y=y, label=ward), color = "white") +
  facet_wrap(~last_name) +
  scale_fill_viridis_c("% of vote") +
  theme_map_sixtysix() +
  coord_map(xlim=bbox[1,], ylim=bbox[2,]) +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal"
  ) +
  ggtitle("Percent of the Vote in Northwest Wards", "2017 DA Race")

da_fit <- lm(
  pvote ~
    CANDIDATE * pct_wht + 
    CANDIDATE * pct_blk + 
    CANDIDATE * pct_hisp,
    # CANDIDATE * log(pop_median_income_2017),
  data = da_results
)

da_results$predicted <- predict(da_fit, newdata = da_results)
da_results$resid <- with(da_results, pvote - predicted)

neighbors <- st_intersection(st_as_sf(sp_divs), st_as_sf(sp_divs))
neighbors <- neighbors %>%
  filter(WARD != WARD.1)

neighbors <- neighbors %>%
  mutate(geometry_type = st_geometry_type(geometry)) %>%
  filter(!geometry_type %in% c("POINT", "MULTIPOINT"))


neighbors <- neighbors %>% 
  mutate(
    WARD.0 = sprintf("%02d", asnum(WARD)),
    WARD.1 = sprintf("%02d", asnum(WARD.1))
  ) %>%
  left_join(
    endorsements %>% rename(endorsement.0 = endorsement),
    by = c("WARD.0" = "ward")
  ) %>%
  left_join(
    endorsements %>% rename(endorsement.1 = endorsement),
    by = c("WARD.1" = "ward")
  ) 

neighbors <- neighbors %>% 
  left_join(
    da_results %>% 
      select(WARD_DIVSN, last_name, total_votes, pvote, resid) %>%
      rename(total_votes.0 = total_votes, pvote.0 = pvote, resid.0 = resid),
    by = c("WARD_DIVSN" = "WARD_DIVSN", "endorsement.0" = "last_name")
  ) %>% 
  left_join(
    da_results %>% 
      select(WARD_DIVSN, last_name, total_votes, pvote, resid) %>%
      rename(total_votes.1 = total_votes, pvote.1 = pvote, resid.1 = resid),
    by = c("WARD_DIVSN.1" = "WARD_DIVSN", "endorsement.0" = "last_name")
  )

library(lme4)

df0 <- neighbors %>% filter(endorsement.0 != endorsement.1)

fit_lmer <- function(neighbor_df){
  re_fit <- lmer(
    resid.0 - resid.1 ~ (1 | WARD.0),
    data = neighbor_df,
    weights = neighbor_df %>%
      with(pmin(total_votes.0, total_votes.1))
  )

  re <- ranef(re_fit)$WARD.0 
  re <- re %>%
    mutate(
      ward = row.names(re),
      effect = re_fit@beta + `(Intercept)`
    )

  return(
    list(
      fit = re_fit,
      re = re
    )
  )
}

fit_baseline <- fit_lmer(df0)

n_boot <- 200
bs_list <- vector(mode = "list", length = n_boot)
for(b in 1:n_boot){

  sample_divs = sample(unique(df0$WARD_DIVSN), replace = TRUE)

  #if(b %% floor(n_boot / 10) == 0) print(b)
  df_samp <- data.frame(WARD_DIVSN = sample_divs) %>% left_join(df0)
  bs_fit <- fit_lmer(df_samp)
  bs_list[[b]] <- bs_fit
}

fixef_ci <- quantile(
  sapply(bs_list, function(x) fixef(x$fit)),
  c(0.025, 0.975)
)
cat(paste0(
  "Average Effect of a Ward Endorsement:\n",
  sprintf(
    "%0.1f (%0.1f, %0.1f)",
    fixef(fit_baseline$fit)["(Intercept)"] * 100,
    fixef_ci[1] * 100,
    fixef_ci[2] * 100
  )
))

ranef_ci <- bind_rows(
  lapply(bs_list, function(x) x$re), 
  .id = "sim"
) %>%
  group_by(ward) %>%
  summarise(
    p025 = quantile(effect, 0.025),
    p975 = quantile(effect, 0.975)
  )

fit_baseline$re %>% 
  select(ward, effect) %>% 
  left_join(ranef_ci) %>%
  mutate(
    ci = sprintf("(%0.1f, %0.1f)", 100 * p025, 100*p975)
  ) %>%
  left_join(
    da_results %>% 
      group_by(WARD16, last_name) %>%
      summarise(
        pvote = 100 * weighted.mean(pvote, w = total_votes),
        total_votes = sum(total_votes)
      ) %>%
      inner_join(
        endorsements, 
        by = c("WARD16" = "ward", "last_name" = "endorsement")
      ),
    by = c("ward" = "WARD16")
  ) %>%
  mutate(
    pvote = round(pvote, 0),
    effect = round(100 * effect, 0),
    vote_effect = round(effect/100 * total_votes)
  ) %>%
  rename(endorsement = last_name) %>%
  select(ward, endorsement, pvote, effect, ci, total_votes, vote_effect) %>%
  arrange(desc(vote_effect)) %>% 
  DT::datatable(
    rownames=FALSE,
    colnames = c("Ward", "Endorsee", "% of Vote in Ward", "Endorsement Effect at Boundary","CI", "Ward Votes", "Vote Effect of Endorsement")
  )

ggplot(ggwards, aes(x = long, y=lat)) +
  geom_polygon(aes(group=group), fill = strong_green, color = "white") +
  geom_text(data = ward_centroids, aes(x=x, y=y, label=ward), color = "white") +
  theme_map_sixtysix() +
  coord_map() +
  ggtitle("Philadelphia's Wards")