---
title: "Week 6 Lecture"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
titanic <- read.csv("http://guerzhoy.princeton.edu/201s20/titanic.csv")
```


### Cross-validation


We sometimes want to try different models, and pick the best one. Generally, we would evaluate the models by measuring the performance on the *validation* set, pick the best one, and then evaluate the final model on the test set.

If we had used the validation set as the test set, we would be overestimating how well the model would perform on new data -- basically, we'd be overfitting on the test set.


### Trying out different models: variable selection

Let's try to predict survival on the Titanic again. We'll try different things -- `Pclass`, `as.factor(Pclass)`, `log(Age)`...

We'll use the validation set to select the model that works best on the validation set. To get an estimate of the performance on "new" data (i.e., data from the Titanic that we haven't looked at yet), we'll use the test set.

```{r}
set.seed(0)
idx <- sample(1:length(titanic$Survived))
titanic.train <- titanic[idx[1:50],]
titanic.valid <- titanic[idx[300:600],]
titanic.test <- titanic[idx[601:800],]
```


```{r}
GetCorrectClassificationRate <- function(fit, dat){
  pred <- predict(fit, newdata=dat, type = "response") > .5
  mean(pred == dat$Survived)
}

fit0 <- glm(Survived ~ 1, family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit0, titanic.valid)

fit1 <- glm(Survived ~ Sex + Age + Pclass, family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit1, titanic.valid)

fit2 <- glm(Survived ~ Sex + Age + as.factor(Pclass), family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit2, titanic.valid)

fit3 <- glm(Survived ~ Sex + log(Age), family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit3, titanic.valid)

fit3 <- glm(Survived ~ Sex + Age + Fare, family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit3, titanic.valid)

fit4 <- glm(Survived ~ Sex + Age + Fare + Pclass, family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit4, titanic.valid)

fit5 <- glm(Survived ~ Sex + Age + Fare + Pclass + Parents.Children.Aboard, family = binomial, data = titanic.train)
GetCorrectClassificationRate(fit5, titanic.valid)


GetCorrectClassificationRate(fit4, titanic.test)
```





### `replicate` 

`replicate` can be used to repeat the same call multiple times. Unless there is randomness involved, that's not very useful: you'll get the same output every time. But with randomness involved, we can obtain important insights.

First, let's see the syntax:

```{r}
square <- function(x){
  x**2
}

replicate(10, square(5))
```

Now, let's try this with randomness:

```{r}
CoinFlip <- function(){
  sample(c("Heads", "Tails"), size = 1)
}
set.seed(0)
replicate(10, CoinFlip())
```

We can now use `replicate` to see how the classification performance on the `titanic` dataset depends on the selection of the training set.

```{r}
GetCorrectClassificationRate <- function(fit, dat){
  pred <- predict(fit, newdata=dat, type = "response") > .5
  mean(pred == dat$Survived)
}


GetValidPerformance <- function(valid.size, titanic){
  idx <- sample(1:length(titanic$Survived))
  titanic.train <- titanic[idx[1:500],]
  titanic.valid <- titanic[idx[501:(501 + valid.size)],]
  fit <- glm(Survived ~ Pclass,family = binomial, data = titanic.train)
  GetCorrectClassificationRate(fit, titanic.valid)
}

GetValidPerformance.fixtrained <- function(train.idx, rest.idx, valid.size, titanic){
  
  titanic.train <- titanic[train.idx,]
  valid.idx <- sample(rest.idx)[1:valid.size]
  titanic.valid <- titanic[valid.idx,]
  fit <- glm(Survived ~ Pclass,family = binomial, data = titanic.train)
  GetCorrectClassificationRate(fit, titanic.valid)
}


set.seed(0)
idx <- sample(1:length(titanic$Survived))
train.idx <- idx[1:500]
rest.idx <- idx[501:length(idx)]

perf.data <- data.frame(valid.perf = replicate(1000, GetValidPerformance.fixtrained(train.idx, rest.idx, 50, titanic)))

ggplot(data = perf.data, mapping = aes(x=valid.perf)) + 
  geom_bar()
```