---
title: "Titanic -- Part 2"
output: html_document
---

```{r setup, include=TRUE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(ggplot2)
require(tigerstats)
require(xtable)
require(lattice)

#CHANGE THIS
#setwd("c:/STA303/lectures/lec5_pre")
setwd("/media/guerzhoy/Windows/STA303/lectures/lec5_pre")

titanic <- read.csv("titanic.csv")
titanic <- titanic[!is.na(titanic$age),]

```
## Drop in Deviance Tests

```{r}
full <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
reduced <- glm(survived~sex, family=binomial(link="logit"), data=titanic)
anova(reduced, full, test="LRT")
```
Let's compute this qunatity manually:

```{r}
1-pchisq(0.6694, df=1)
```

So we cannot reject the hypothesis that the coefficient for age is 0. (Of course, this is also confirmed by the Wald Test.)


Can we reject the hypothesis that the class doesn't matter?

```{r}
full <- glm(survived~age+pclass, family=binomial(link="logit"), data=titanic)
reduced <- glm(survived~age, family=binomial(link="logit"), data=titanic)
anova(reduced, full, test="LRT")
```

Apparently yes.

# AIC

Let's try AIC for 

```{r}
full <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
full
reduced <- glm(survived~sex, family=binomial(link="logit"), data=titanic)
reduced
```

AIC says use the reduced model. Intuition: the reduction in deviance is not large enough to compensate for the additional parameter. The Residual Deviance (just 2*Deviance in the slides) goes down by 1, but the number of parameters goes up by 1 (and gets multiplied by 2) when we go to the full model.

## Classification

Let's load the iris dataset.

```{r}
library(datasets)
data("iris")
iris_small <- subset(iris, Species=="virginica" | Species=="versicolor")
```

Plotting the iris dataset:

```{r}
ggplot(iris_small, aes(Petal.Width, Petal.Length)) + geom_point(aes(colour=factor(Species)))
```

Prediction:

```{r}
fit <- glm(Species~Petal.Width+Petal.Length, data=iris_small, family="binomial")

new_data <- data.frame(Petal.Width=1.75, Petal.Length=5)
plogis(predict(fit,new_data))
```

Predictions
```{r}
p_cutoff <- .5
iris_small$pred <- ifelse(fit$fitted.values > p_cutoff, "virginica", "versicolor")
```

How well did we do? We are computing how well the predictions match to the correct classes.

```{r}
mean(iris_small$pred==iris_small$Species)
```

Not bad!


Let's try this with the Titanic data

```{r}
p_cutoff <- 0.5  #cutoff for the pi
fit <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
titanic$pred <- ifelse(fit$fitted.values > p_cutoff, 1, 0)
mean(titanic$pred==titanic$survived)
```

What about the cut-offs? Let's say that you want to sell insurance to someone. Maybe it's more constly (because of the cost of risk) to sell insurance to someone who will die than to refuse to sell insurance to someone who will live.

```{r}
p_cutoff <- 0.5
fit <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
titanic$pred <- ifelse(fit$fitted.values > p_cutoff, 1, 0)

cost_sell_die <- 10
cost_refuse_survive <- 1

total_cost_sell_die <- cost_sell_die*sum((titanic$pred==1) & (titanic$survived==0))
total_cost_refuse_survive <- cost_refuse_survive*sum((titanic$pred==0) & (titanic$survived==1))
total_cost_sell_die+total_cost_refuse_survive
```

Let's try this with two difference cutoffs:

```{r}
p_cutoff <- 0.75
fit <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
titanic$pred <- factor(ifelse(fit$fitted.values > p_cutoff, 1, 0))

cost_sell_die <- 10
cost_refuse_survive <- 1

total_cost_sell_die <- cost_sell_die*sum((titanic$pred==1) & (titanic$survived==0))
total_cost_refuse_survive <- cost_refuse_survive*sum((titanic$pred==0) & (titanic$survived==1))
total_cost_sell_die+total_cost_refuse_survive
```


```{r}
p_cutoff <- 0.5
fit <- glm(survived~sex+age, family=binomial(link="logit"), data=titanic)
titanic$pred <- factor(ifelse(fit$fitted.values > p_cutoff, 1, 0))

cost_sell_die <- 10
cost_refuse_survive <- 1

total_cost_sell_die <- cost_sell_die*sum((titanic$pred==1) & (titanic$survived==0))
total_cost_refuse_survive <- cost_refuse_survive*sum((titanic$pred==0) & (titanic$survived==1))
total_cost_sell_die+total_cost_refuse_survive
```

Can also use the confusion matrix:
```{r}
with(titanic, table(pred, survived))
```

## Overfitting: Introduction (to be continued Wednesday)
```{r}
titanic_small <- titanic[1:20,]
fit<-glm(survived~pclass+sex+age+sibsp+parch+fare, data=titanic_small)
p_cutoff <- 0.5
titanic$pred <- ifelse(predict(fit, titanic) > p_cutoff, 1, 0)
titanic_small$pred <- ifelse(predict(fit, titanic_small) > p_cutoff, 1, 0)
mean(titanic_small$pred==titanic_small$survived, na.rm=TRUE)
mean(titanic$pred==titanic$survived, na.rm=TRUE)
```
Cross-validation
Split your dataset into a "training set" and a "validation set"
fit your model on the training set, and you check the performance on the "validation set"


