---
title: "Titanic -- Extrabionmial Variation"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown
```{r}
#setwd("/media/guerzhoy/Windows/STA303/lectures/lec5_pre")
setwd("c:/STA303/lectures/lec5_pre")

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

Let's cook up a dataset with overdispersion. Suppose that we have 200 women in first class, whose probability of survival is 90%, and 200 women in second class, whose probability of survival is 10%

Suppose men's probability of survival is 50%, and there are 50 men

Suppose we haven't recorded the class information.

```{r}
survived<-c(rbinom(200, size=1, prob=.9), rbinom(200, size=1, prob=.1), rbinom(50, size=1, prob=.5))
survived<-c(rbinom(200, size=1, prob=.4), rbinom(200, size=1, prob=.4), rbinom(50, size=1, prob=.5))
sex  <-   c(rep("Female", 200)          , rep("Female", 200)          , rep("Male", 25)   , rep("Male", 25))
class<-   c(rep(1, 200)                 , rep(2, 200)                 , rep(1, 25)        , rep(2, 25))

df<- data.frame(survived=survived, sex=sex, class=class)


fit_quasi <- glm(survived~sex, family=quasibinomial, data=df)
summary(fit_quasi)
```

Here's a more direct demonstration:
```{r}
x<-c(rbinom(100, size=1, prob=.9), rbinom(100, size=1, prob=.1)) ; c(var(x),mean(x)*(1-mean(x)))
```

As you can see, it's usually not a huge effect, with just the one predictor variable.

## Example of overdispersion in the Titanic dataset
Omitting the last name, using the first letter of the last name

```{r}
titanic$first_letter <- sapply(titanic$name, function(name) substr(name, 1,1))
titanic <- titanic[!(titanic$first_letter %in% c('v', 'U', 'Y')),]
summary(glm(survived~sex+log(age)+first_letter, family=quasibinomial, data=titanic))
```

As you can see, the overdispersion parameter here is 1.06. This makes the results slightly less significant.


Let's see which residual are responsible for this

```{r}
fit<-glm(survived~sex+log(age)+first_letter, family=quasibinomial, data=titanic)
plot(fit, 1)
titanic[which(fit$residuals>10),]
```