## Preliminaries

Libraries:
```{r}
library(tidyverse)
library(caret)
```

Helper functions:
```{r}
# e^x/(e^x+1) for input x
expit <- function(x) exp(x)/(exp(x) + 1)

# root mean squared error for true y and prediction yhat
rmse <- function(y,yhat) sqrt(mean((y-yhat)^2))
# accuracy for predictions yhat and true binary outcome y
acc <- function(y, yhat) mean(y == yhat)
```

## Cross-Validation

Split the data frame into 4 parts. This sets up the indexing.
```{r}
faithfulGrid <- data.frame(waiting = seq(35,105))
n <- dim(faithful)[1]
k <- 4
randomInd <- sample(n)
splitSize <- n/k  # Only works if n/k is an integer. Use caret package!
```

Train on 3/4 of data 4 times, leaving out a different 1/4 for testing each time. Here we plot the 4 estimated models. The indexing is an absolute nightmare, as you can see, and this is the special case when n/k is an integer. Indexing with a minus sign first *removes* indices.
```{r}
linMods <- lapply(1:k, function(i) lm(eruptions ~ waiting,
data = faithful[randomInd[-((splitSize*(i-1)+1):(splitSize*i))],]))

deg <- 20
polMods <- lapply(1:k, function(i) lm(eruptions ~ poly(waiting,deg),
data = faithful[randomInd[-((splitSize*(i-1)+1):(splitSize*i))],]))

par(mfrow=c(2,2))
for (i in 1:k) {
  plot(faithful$waiting, faithful$eruptions)
  lines(faithfulGrid$waiting, predict(linMods[[i]], faithfulGrid), col = 'red')
  lines(faithfulGrid$waiting, predict(polMods[[i]], faithfulGrid), col = 'blue')
}
```

RMSEs from CV: test on the left-out 1/4 of the data. You will not have to do this indexing yourself.
```{r}
linRMSEs <- sapply(1:k,
                   function(i) rmse(faithful$eruptions[randomInd[(splitSize*(i-1)+1):(splitSize*i)]],
                                    predict(linMods[[i]],
                                            faithful[randomInd[(splitSize*(i-1)+1):(splitSize*i)],])))

polRMSEs <- sapply(1:k,
                   function(i) rmse(faithful$eruptions[randomInd[(splitSize*(i-1)+1):(splitSize*i)]],
                                    predict(polMods[[i]],
                                            faithful[randomInd[(splitSize*(i-1)+1):(splitSize*i)],])))
mean(linRMSEs)
mean(polRMSEs)
```

Caret package makes CV (cross-validation easier).
```{r}
folds4 <- createFolds(sample(n), k)
linModsCaret <- lapply(folds4, function(fold) lm(eruptions ~ waiting,
                                                 data = faithful[-fold,]))
linRMSEsCaret <- sapply(1:k, function(i) rmse(faithful$eruptions[folds4[[i]]],
                                              predict(linModsCaret[[i]],
                                                      faithful[folds4[[i]],])))

mean(linRMSEsCaret)
```
Still isn't pretty, though.


## Binary Outcomes

We know how to handle binary and categorical predictor variables. How about outcomes?
```{r}
votes115 <- read_csv('../data/S115_votes.csv')
members115 <- read_csv('../data/S115_members.csv')
joint115 <- inner_join(votes115, members115, by = 'icpsr')  # link records on senator ID
dim(joint115)
unique(joint115$cast_code)
joint115 <- joint115 %>% filter(cast_code %in% c(1,6),
                                party_code %in% c(100,200))
joint115$party <- ifelse(joint115$party_code == 100, 'D', 'R')
joint115$yay <- (joint115$cast_code == 1)
```

How likely are senators from each party to vote to pass bills? 115th congress (2017-2019).
```{r}
# Logistic regression
partyMod <- glm(yay ~ party, data = joint115, family = 'binomial')
summary(partyMod)
```
Which party controlled the Senate 2017-2019?

What do the coefficient estimates mean?
```{r}
cat('Prob D votes yay:', expit(coef(partyMod)[1]),
    '\nProb R votes yay:', expit(sum(coef(partyMod))), '\n\n')
```

We didn't need a fancy model for this, actually:
```{r}
joint115 %>% group_by(party) %>% summarize(yayProb = mean(yay))
```

It gets interesting when continuous covariates are included:
```{r}
nomMod <- glm(yay ~ nominate_dim1*nominate_dim2, 
              data = joint115, family = 'binomial')
summary(nomMod)
```

Predict for a different congress (103rd, 1993-1995). Extrapolation or no?
```{r}
votes103 <- read_csv('../data/S103_votes.csv')
members103 <- read_csv('../data/S103_members.csv')
joint103 <- inner_join(votes103, members103, by = 'icpsr')
dim(joint103)
unique(joint103$cast_code)
joint103 <- joint103 %>% filter(cast_code %in% c(1,6),
                                party_code %in% c(100,200))
joint103$party <- ifelse(joint103$party_code == 100, 'D', 'R')
joint103$yay <- (joint103$cast_code == 1)
```

```{r}
predProbsParty <- expit(predict(partyMod, joint103))
acc(joint103$yay, predProbsParty > 0.5)

predProbsNom <- expit(predict(nomMod, joint103))
acc(joint103$yay, predProbsNom > 0.5)

mean(joint115$yay)
acc(joint103$yay, 1)
```

## Plot Generating for Slides

In case you're curious how I make slides.

Linear regression:
```{r}
x <- rnorm(100)
y <- 1 + x + rnorm(100)
plot(x, y, main = '')
abline(1, 1, col = 'red')
```

Logistic regression:
```{r}
x <- rnorm(100)
probs <- expit(1 + x + rnorm(100))
y <- rbinom(100, 1, probs)
plot(x, y, main = '', ylab = 'y, p-hat')
lines(seq(-3,3,.1), expit(1 + seq(-3,3,.1)), col = 'red')
```

