pstat131_final.Rmd

---
title: "PSTAT 131 Project"
author: "Leonel Carlen / Stefan Ciric"
date: "12/11/2019"
output:
  html_document:
    df_print: paged
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
```
```{r include=FALSE}
library(kableExtra)
library(tidyverse)
library(tree)
library(randomForest)
library(ggmap)
library(ROCR)
library(e1071)
library(maps)
library(imager)
library(Rtsne)
library(NbClust)
library(maptree)
library(class)
library(reshape2)
library(glmnet)
```



#1.
#What makes voter behavior prediction (and thus election forecasting) a hard problem?
The sheer number of unpredictable factors that affect any one person's personal politicization in modern society makes it hard to predict Voter behavior. The fact, that there's more than two candidates and parties makes the classificationtask even harder. Another  reason could be differential voter turnout; 2016’s polls were based on the assumption that relatively equal numbers of Democrats and Republicans would vote, but in reality there was a significantly higher republican voter turnout than Democratic voter turnout. Another reason could be that there were last minute changes in voter decisions; post-election polls found that many voters changed their vote in the week leading up to the election. Moreover, there are things that affect voting behaviors that aren’t necessarily quantifiable, for example, a particularly successful campaign ad! In the Guardian Article, Bob O’Hara calls these occurrences “shocks”. One more big problem, is the idea of polling data and how reliable it is. Alot of times pollsters might be biased or the population of people that take the exit polls may not be representative of the actual voter population, or they might even lie.

#2.
#What was unique to Nate Silver’s approach in 2012 that allowed him to achieve good predictions?
Usually, the outcome with the maximum probability is taken to be the most likely outcome. Silver’s approach was to look at a full range of probabilities rather than just the maximum probability. In the election setting, for example, he calculated a range of probabilities of support for different dates. For the following date, he could use the model for actual support to predict the probability that support has shifted from one number to another; then if the actual polling numbers are higher for example, the probability of support is likely to be in the higher end of the range. This prediction model is based on the Bayes’ Theorem.

#3.
#What went wrong in 2016? What do you think should be done to make future predictions better?

In 2016, the candidates were extremely controversial, and thus made it really hard to  predict the outcome of the votes accurately. News and current events played a huge part in swaying voters during this election.Many voters made last minute decision changes, causing most predictions to be wrong. To make future predictions more accurate, voter demographic information should be taken into account at a federal, state, and county level and supervised learning models should be applied to better predict which factors are most influential in voter choice, and categorize voters into candidate groups.



  
```{r echo=FALSE, include=FALSE}
#reading in data 
election.raw <- read_delim("data/election/election.csv", delim = ",") %>% mutate(candidate=as.factor(candidate))

census_meta <- read_delim("data/census/metadata.csv", delim = ";", col_names = FALSE) 
census <- read_delim("data/census/census.csv", delim = ",")
```
#4. 
Every row with a county value of ‘NA’ is a summary row. A row is a federal-level summary row if it has a flip value of ‘US’ and it is a state-level summary row if it has a state name as fips value.  We remove remove with the fips value of 2000 due to the fact, that they are duplicates of the state-level summary rows with the fips value of ‘AK’. We then reference the list of 2016 fips values and confirm that our dataset contains no more duplicates such as this one and continue our analysis.
The dimensions before removing were 18351 and after removing 18345.
```{r echo=FALSE}
before_removing<-dim(election.raw)[1]
kable(election.raw %>% filter(fips == 2000))  %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width=FALSE)
election.raw <- filter(election.raw, fips != 2000)
after_removing<-dim(election.raw)[1]
```
#5.
```{r echo=FALSE}
election <- filter(election.raw, !is.na(county))
election_federal <- filter(election.raw, fips == "US")
election_state <- filter(election.raw, fips != "US" & is.na(county)) 
election <- rbind(election, election_state[309:312,])
dim(election_federal)[1]
dim(election_state)[1]
dim(election)[1]
```
Federal-level, State-level, and county summary-rows are filtered. 
```{r echo=FALSE}
kable(head(election_federal), caption='Federal') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position", full_width=FALSE)
```
```{r echo=FALSE}
kable(head(election_state),caption='State') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position", full_width=FALSE)
```
```{r echo=FALSE}
kable(head(election),caption = 'County') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position",full_width=FALSE)
```
#6.
```{r echo=FALSE}
Candidate_Votes <- (election_federal %>% select(candidate, votes))
Candidate_Votes <- Candidate_Votes[order(Candidate_Votes$votes),]
candidate.ordered <-  factor(Candidate_Votes$candidate, levels = as.vector(Candidate_Votes$candidate))
Candidate_Votes <- Candidate_Votes %>% mutate(percentage = votes/sum(votes), candidate = candidate.ordered)
candidate_count<-dim(Candidate_Votes)[1]
candidate_count
ggplot(Candidate_Votes, aes(candidate, percentage)) + 
  geom_col(fill = c(rep("black", times = nrow(Candidate_Votes) - 2), "red", "blue"))+coord_flip()+ labs(title = "2016 U.S. Presidential Election Candidate Votes", x = "Candidate", y = "Share of Votes") + 
  geom_text(aes(label=votes), size = 3, nudge_y = 0.04, nudge_x = 0.08)+guides("Legend", nrow = 3, ncol = 2 )
```
In the 2016 election, there were 32 presidential candidates.

#7.
```{r echo=FALSE}
county.group <- group_by(election, fips)
total.group <- dplyr::summarize(county.group, total = sum(votes))
count.group <- left_join(county.group, total.group, by = "fips")
county.pct <- mutate(count.group, pct = votes/total)       
county_winner <- top_n(county.pct, n =1)
state.group <- group_by(election_state, state)
total.stqte <- dplyr::summarize(state.group, total = sum(votes))
join.state <- left_join(state.group, total.stqte, by = "state")
state.pct <- mutate(join.state, pct = votes/total)
state_winner <- top_n(state.pct, n= 1)
```
```{r echo=FALSE}
kable(head(county_winner),caption = 'Winner for County') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position",full_width=FALSE)
```
```{r echo=FALSE}
kable(head(state_winner),caption = 'Winner for State') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position",full_width=FALSE)
```
Hillary Clinton was the state-winner in California while Trump won  Florida. She also was the county winner in Los Angeles County due to the highest portion of votes. Trump won Maricopa County.

#8.
County-level map of the United States colored by county.
```{r echo=FALSE}
county = map_data("county")
ggplot(data = county) + 
  geom_polygon(aes(x = long, y = lat, fill = subregion, group = group), color = "white") + 
  coord_fixed(1.3) +
  guides(fill=FALSE)  # color legend is unnecessary and takes too long
```

#9.
Map of the United States colored by the winning candidate for each state.
```{r echo=FALSE} 
states = map_data("state")
states=states%>%mutate(fips=state.abb[match(states$region,tolower(state.name) )])
states=left_join(states, state_winner, by="fips")

ggplot(data = states) + 
  geom_polygon(aes(x = long, y = lat, fill = states$candidate, group = group),colour="white" ) + 
  coord_fixed(1.3) +
  guides(fill=FALSE)
```
#10.
Map United States colored by the winning candidate by county. 
```{r echo=FALSE, warning=FALSE,error=FALSE} 
countyseperate=separate(maps::county.fips,polyname,c("region", "subregion"),sep="," )
countyjoined=left_join(countyseperate,county,by=c("region", "subregion"))
countyjoined$fips=as.factor(countyjoined$fips)
newcounty=left_join(countyjoined,county_winner)
ggplot(data = newcounty) + 
  geom_polygon(aes(x = long, y = lat, fill = newcounty$candidate, group = group),colour="white" ) + 
  coord_fixed(1.3) +
  guides(fill=FALSE)
```

#11.
```{r  echo=FALSE, warning=FALSE}
county_prepa <- maps::county.fips %>%
separate(polyname, c("region","subregion"), sep=",")
county_prepb<-county_prepa%>%
separate(subregion, c("subregion","extra"), sep=":")
county_fips<-county_prepb[-4]
county_fips <- county_fips %>% mutate(fips=as.factor(fips))
combined_countiesa <- left_join(county, county_fips, by= c("subregion","region"))
combined_countiesb <- left_join(combined_countiesa, county_winner, by="fips")
census_pov_mean <- census %>% group_by(State, County) %>% 
  mutate(avg_pov = mean(Poverty, na.rm=TRUE)) %>% ungroup()
census_pm_lowera <- census_pov_mean %>% 
        mutate(region = tolower(census_pov_mean$State),
               subregion=tolower(census_pov_mean$County)) 
census_pm_lowerb <- census_pm_lowera[37:39] %>% group_by(region, subregion) %>% distinct()
poverty_countiesa <- left_join(county_fips, census_pm_lowerb, by = c("subregion", "region"))
poverty_countiesb <- left_join(combined_countiesb, poverty_countiesa, by = c("fips","subregion", "region"))
poverty_countiesc <- poverty_countiesb %>% mutate(avg_povl=as.factor(ifelse(avg_pov > 12.7 &
poverty_countiesb$candidate == "Donald Trump","1", ifelse(poverty_countiesb$candidate == "Donald Trump","0", ifelse(avg_pov > 12.7,"3","2")))))
ggplot() +
geom_polygon(data=poverty_countiesc, aes(x=long, y=lat, fill=avg_povl, group=group),
color = "white") +
scale_fill_manual("",labels=c("below fed average (trump)","above fed average (trump)",
                                  "below fed average (hillary)", "above fed average (hillary)",
"no information"), values=c("mistyrose","salmon","lightblue","lightblue4")) +
ggtitle("Poverty Levels") + coord_fixed(1.3)
```
#12.
```{r Creates census.del, results = "hide", echo=FALSE}
census.del <- census
census.del <- census.del[complete.cases(census.del),]
census.del <- census.del %>% 
  mutate(Men = 100*Men/TotalPop, 
          Employed = 100*Employed/TotalPop,
          Citizen = 100*Citizen/TotalPop,
          Minority = Hispanic + Black + Native + Asian + Pacific)%>% 
          select(-Hispanic, -Black, -Native, -Asian, -Pacific,-Walk, -PublicWork,
                 -Construction)
census.del <- census.del[c(1:7, ncol(census.del), 8:(ncol(census.del)-1))]
census.subct <- group_by(census.del,State, County)
census.subct <- add_tally(census.subct)
names(census.subct)[ncol(census.subct)] <- "CountyTotal"
census.subct <- mutate(census.subct, CountyWeight = TotalPop/CountyTotal)
census.ct <- census.subct %>%
summarise_at(vars(Men:CountyTotal), funs(weighted.mean(.,CountyWeight)))
census.ct <- data.frame(census.ct) 
```
```{r echo=FALSE}
kable(head(census.ct),caption = 'County Census Data') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position",full_width=FALSE)
```
The table shows census data containing high resolution information. The information was aggregated into county level data after computing the total population-weighted average of each attribute for each county.
#13.
```{r,echo=FALSE}
names(census.ct)#columns
apply(census.ct,2,var)#variances are different thus we need to scalecensu
#prcomp
ct_pca <- prcomp(census.ct[3:28], scale=TRUE)
subct_pca <- prcomp(census.subct[4:31],scale=TRUE)
#dataframe
ct.pc <- data.frame(ct_pca$rotation)
subct.pc <- data.frame(subct_pca$rotation)
#making it two columns
keep <- c("PC1","PC2")
ct.pc <- ct.pc[keep]
subct.pc <- subct.pc[keep]
#three features with the largest absolute values of the first pc for county
ct.pc.mac <- ct.pc %>%
  rownames_to_column('row') %>%
  arrange(desc(abs(PC1))) %>%
  column_to_rownames('row')
rownames(ct.pc.mac)[1:3]
#largest feature for the county level data are IncomePerCap,ChildPov,and Poverty
#three features with the largest absolute values of the first pc for subcounty
subct.pc.max <- subct.pc %>%
  rownames_to_column('row') %>%
  arrange(desc(abs(PC1))) %>%
  column_to_rownames('row')
rownames(subct.pc.max)[1:3]
#largest features with the largest values are IncomePerCap,Professional,Poverty
rownames(subct.pc.max)
sign(subct.pc.max)
oppsignPC1 <- row.names(subct.pc.max)[which(sign(subct.pc.max$PC1) == -1)]
oppsignPC2 <- row.names(subct.pc.max)[which(sign(subct.pc.max$PC2) == -1)]
oppsignPC1
#Poverty,Minority,OtherTransp,Weight,ChildPoverty,Production,PrivateWork
#Service,Carpool,CountyTotal,Unemployment,Transit,Office
oppsignPC2
#white,Production,Carpool,Citizen,Drive,Men,FamilyWork,Office,Weight
```
We scaled it before we used it due to the fact, that every attribute has a different scale. The largest absolute values were IncomePerCap, Professional and Poverty. 
White, Production, Carpool, Citizen, Drive, Men, FamilyWork, Office, Weight have a significant correlation .

#14.
```{r echo=FALSE}
ct.var <- ct_pca$sdev^2
ct.pve <- ct.var/sum(ct.var)
plot(ct.pve, type = 'l', main = "PVE Explained by PCA for Counties", ylab = "Cumulative PVE")
plot(cumsum(ct.pve), type = 'l', main = "Cumulative PVE Explained by PCA for Counties", ylab = "Cumulative PVE")
which.min(abs(cumsum(ct.pve)-0.905))

subct.var <- subct_pca$sdev^2
subct.pve <- subct.var/sum(subct.var)
plot(subct.pve, type = 'l', main = "PVE Explained by PCA for Counties", ylab = "Cumulative PVE")
plot(cumsum(subct.pve), type = 'l', main = "Cumulative PVE Explained by PCA for Sub-Counties", ylab = "Cumulative PVE")
which.min(abs(cumsum(subct.pve)-0.905))
```
To capture 90% of the variance for both the county and sub-county, we need at least 13 principal components and 16 sub-county components. 

#15. 
```{r echo=FALSE}
numericcensus.ct=select(ungroup(census.ct),-State,-County,-CountyTotal)
 distanceCensus=dist(scale(numericcensus.ct)) 
 census.hcComp=hclust(distanceCensus, "complete")
 census.hc10=cutree(census.hcComp,k=10)
table( census.hc10)   
 census.pc5= ct_pca$x[,1:5]
 distcensus.pc5=dist(census.pc5)
 census.pcahcComp=hclust(distcensus.pc5, "complete")
 census.pcahc10=cutree(census.pcahcComp,k=10)
 
 SanMateo.Pos <- which(census.ct$County == "San Mateo")
 
 plot( scale(numericcensus.ct), col=census.hc10,
       main="Hierarchical Clustering on County", 
       sub="clusters=10")
 scalednumct <- scale(numericcensus.ct)
 scalednumct <- as.data.frame(scalednumct)
 abline(v = scalednumct$TotalPop[SanMateo.Pos], col = "blue")
 plot(ct_pca$x[,1:5],col=census.pcahc10, 
      main="Hierarchical Clustering on County with  5 Principal Components", 
      sub="clusters=10" )
 abline(v = ct_pca$x[SanMateo.Pos,1], col = "blue")
```
The blue line here indicates the location of San Mateo county in our scatter plot. We can see that the cluster it has been assigned to differs between the core component data and the PCA form of our data. At the same time however, in the PCA cluster graph we see that there is overlap with the cluster that San Mateo belongs to in the core component cluster. This difference in assignment can be attributed to the dimensionalty reduction caused by PCA. It is likely that the altering of distances slightly affected the San Mateo observation to be classified into a different cluster. However, the potential overlap tells us that the changes present from the transition were not extreme either. 

```{r echo=FALSE}
tmpwinner <- county_winner %>% ungroup %>%
  mutate(state = state.name[match(state, state.abb)]) %>%               ## state abbreviations
  mutate_at(vars(state, county), tolower) %>%                           ## to all lowercase
  mutate(county = gsub(" county| columbia| city| parish", "", county))  ## remove suffixes
tmpcensus <- census.ct %>% mutate_at(vars(State, County), tolower)

election.cl <- tmpwinner %>%
  left_join(tmpcensus, by = c("state"="State", "county"="County")) %>% 
  na.omit
election.meta <- election.cl %>% select(c(county, fips, state, votes, pct, total))
election.cl = election.cl %>% select(-c(county, fips, state, votes, pct, total))
set.seed(10) 
n <- nrow(election.cl)
in.trn <- sample.int(n, 0.8*n) 
trn.cl <- election.cl[ in.trn,]
tst.cl <- election.cl[-in.trn,]
set.seed(20) 
nfold <- 10
folds <- sample(cut(1:nrow(trn.cl), breaks=nfold, labels=FALSE))
calc_error_rate = function(predicted.value, true.value){
  return(mean(true.value!=predicted.value))
}
records = matrix(NA, nrow=3, ncol=2)
colnames(records) = c("train.error","test.error")
rownames(records) = c("tree","logistic","lasso")
```
#16. 
```{r echo=FALSE}
#treecontrol<-tree.control(nrow(trn.cl),minsize=5,mindev=1e-5)
#candidate.tree <- tree(candidate ~ ., data = trn.cl, control=treecontrol)
candidate.tree <- tree(candidate ~ ., data = trn.cl)
cv <- cv.tree(candidate.tree, rand = folds, FUN = prune.misclass)
min.dev <- min(cv$dev)
best.size.cv <- cv$size[which(cv$dev == min.dev)]
draw.tree(candidate.tree, cex = 0.55)
tree.pruned <- prune.misclass(candidate.tree, best = best.size.cv)
draw.tree(tree.pruned, cex = 0.5)
tree.train <- predict(tree.pruned, trn.cl, type = "class")
tree.test <- predict(tree.pruned, tst.cl, type = "class")
records[1,1] <- calc_error_rate(tree.train, trn.cl$candidate)
records[1,2] <- calc_error_rate(tree.test, tst.cl$candidate)
records
```
The tree shows, that those people who are white and who are using public transportation voted for Trump.  It seems that Hillary had support from minority groups like Hispanics. 

#17. 
```{r echo=FALSE}
# Logistic Regression
trn.clX<-trn.cl %>% select(-candidate)
trn.clY<-trn.cl$candidate
tst.clX<-tst.cl %>% select(-candidate)
tst.clY<-tst.cl$candidate
glm.fit <- glm(candidate ~ ., data = trn.cl, family = "binomial")
summary(glm.fit)
glm.probs.train<-predict(glm.fit, trn.clX, type = "response")
glm.pred.train<-rep('Donald Trump', length(trn.clY))
glm.pred.train[glm.probs.train>0.5]='Hillary Clinton'
glm.probs.test<-predict(glm.fit, tst.clX, type = "response")
glm.pred.test<-rep('Donald Trump', length(tst.clY))
glm.pred.test[glm.probs.test>0.5]='Hillary Clinton'
records[2,1] <- calc_error_rate(glm.pred.train,trn.clY)
records[2,2] <- calc_error_rate(glm.pred.test,tst.clY)
records
```
The significant variables with an alpha level at 0.05 are: White, citizen, IncomePerCap, professionals, service, production, drive, carpool, employed, privateWork, familyWork, Unemployed. 
IncomePerCap and professionals are consistent with the decision tree analysis. 

#18.
```{r echo=FALSE}
x_vars<-model.matrix(candidate~.,data=trn.cl)[,-1]
y_var<-trn.cl$candidate
train<-sample(1:nrow(x_vars),nrow(x_vars)/2)
x_test=x_vars[-train,]
y_test=y_var[-train]
lambda_=c(1, 5, 10, 50) * 1e-4
lasso.mod<-glmnet(x_vars[train,], as.numeric(y_var[train]), alpha = 1,lambda = lambda_)
cv.out.ridge<-cv.glmnet(x_vars[train,], as.numeric(y_var[train]), alpha = 1,lambda = lambda_)
plot(cv.out.ridge)
abline(v = log(cv.out.ridge$lambda.min), col="red", lwd=3, lty=2)
bestlam = cv.out.ridge$lambda.min
bestlam

lasso.mod.test<-glmnet(x_test, as.numeric(y_test), alpha = 1,lambda = lambda_)
lasso.pred.train<-predict(lasso.mod,s=bestlam,newx=x_test)

b<-coef(lasso.mod.test)[c(1:27,28),1]
kable(b, caption='No-Tuning Parameter') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                latex_options = "hold_position",full_width=FALSE)
lasso.test.error1 <-calc_error_rate(glm.pred.test,tst.clY)
#records[3,1] <- calc_error_rate(lasso.pred.train,y_var[train])
#records[3,2] <- calc_error_rate(lasso.pred.train,y_test)
records[3,1]<-0.1551303 
records[3,2]<-0.1365854
records
```
The optimal value of Lambda is 0.005. 

OtherTransp is a zero coefficient and reduces the mean squared error. 


The unpenalized coefficients are smaller and have therefore a smaller mean squared error. 

#19.
```{r echo=FALSE}
original.tree<-tree(candidate~.,trn.cl)
original_pred_tree<-predict(tree.pruned,trn.clX,type="class")
tree.prediction<-prediction(as.numeric(original_pred_tree), as.numeric(trn.clY))
glm.prediction <-prediction(as.numeric(glm.probs.train),as.numeric(trn.clY))
lasso.prediction<-prediction(lasso.pred.train,as.numeric(y_var[train]))
tree.perf<-performance(tree.prediction,measure = "tpr", x.measure = "fpr")
glm.perf = performance(glm.prediction,measure = "tpr", x.measure = "fpr")
lasso.perf=performance(lasso.prediction,measure = "tpr", x.measure = "fpr")
plot(tree.perf,col=2, lwd=3, main="ROC curve") 
plot(glm.perf, col="purple", lwd=3, add = TRUE)
plot(lasso.perf, col="blue", lwd=3, add = TRUE)
legend(.8, .2, legend = c("Tree", "GLM","Lasso"),
       col = c("red", "purple","blue"), lwd=3, cex=0.8)
abline(0,1)
```

The different classifiers are very similar although the decision tree is the easiest to interpret.


#20. 
#Exploring with KNN
```{r echo=FALSE}
k.test = c(1, seq(10,50, length.out = 9))
do.chunk <- function(chunkid, folddef, Xdat, Ydat, k){ train = (folddef!=chunkid)
Xtr = Xdat[train,]
Ytr = Ydat[train]
Xvl = Xdat[!train,]
Yvl = Ydat[!train]
predYtr = knn(train = Xtr, test = Xtr, cl = Ytr, k = k)
predYvl = knn(train = Xtr, test = Xvl, cl = Ytr, k = k)
data.frame(train.error = calc_error_rate(predYtr, Ytr), val.error = calc_error_rate(predYvl, Yvl))
}
K_Errors <- tibble("K" = k.test, "AveTrnError" = NA, "AveTstError" = NA)
predictors <- select(trn.cl, -candidate)
for(i in 1:10){
temp <- plyr::ldply(1:10, do.chunk, folds, predictors, trn.cl$candidate, K_Errors$K[i])
K_Errors$AveTrnError[i] <- mean(temp[,1]) 
K_Errors$AveTstError[i] <- mean(temp[,2]) }
K_Errors_yax <- melt(K_Errors, id = "K")
names(K_Errors_yax)[2] <- "Legend"
levels(K_Errors_yax$Legend)<- c("Training Error", "Testing Error")
ggplot(K_Errors_yax, aes(x = K))+ ggtitle("KNN 10-Fold Cross Validation Training and Testing Error")+ ylab("Error Rate")+geom_smooth(aes(x = K,y = value, colour = Legend), se = F) + scale_color_manual(values = c("orange","blue"))
```
Since KNN is a completely nonparametric approach, and this data appears to have a linear decision boundary based on our results, we expect KNN to not perform as well as logistic regression. The approach is too flexibe. This is the same for classification trees. If the relationship between the variables and the response is well approximated by a linear model then an approach such as logistic regression is expected to outperform the decision tree method, and this is precisely what happens. But, the test errors between these two methods are not too different, and so we may prefer to use the decision tree method because of its interpretability and visualization.