---
title: "Titanic"
output: html_document
---

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

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

Here is the variable descriptions:
```
VARIABLE DESCRIPTIONS:
survival        Survival
                (0 = No; 1 = Yes)
pclass          Passenger Class
                (1 = 1st; 2 = 2nd; 3 = 3rd)
name            Name
sex             Sex
age             Age
sibsp           Number of Siblings/Spouses Aboard
parch           Number of Parents/Children Aboard
ticket          Ticket Number
fare            Passenger Fare
cabin           Cabin
embarked        Port of Embarkation
                (C = Cherbourg; Q = Queenstown; S = Southampton)

SPECIAL NOTES:
Pclass is a proxy for socio-economic status (SES)
 1st ~ Upper; 2nd ~ Middle; 3rd ~ Lower

Age is in Years; Fractional if Age less than One (1)
 If the Age is Estimated, it is in the form xx.5

With respect to the family relation variables (i.e. sibsp and parch)
some relations were ignored.  The following are the definitions used
for sibsp and parch.

Sibling:  Brother, Sister, Stepbrother, or Stepsister of Passenger Aboard Titanic
Spouse:   Husband or Wife of Passenger Aboard Titanic (Mistresses and Fiances Ignored)
Parent:   Mother or Father of Passenger Aboard Titanic
Child:    Son, Daughter, Stepson, or Stepdaughter of Passenger Aboard Titanic

Other family relatives excluded from this study include cousins,
nephews/nieces, aunts/uncles, and in-laws.  Some children travelled
only with a nanny, therefore parch=0 for them.  As well, some
travelled with very close friends or neighbors in a village, however,
the definitions do not support such relations.
```

Let's build some tables to explore the data a little bit

```{r}
mean(titanic$survived)
with(titanic, aggregate(survived, by=list(sex), mean))
with(titanic, aggregate(survived, by=list(pclass), mean))
with(titanic, aggregate(survived, by=list(pclass, sex), mean))
```

What about age? Actually, arguably the best thing to do to visualize the information is to run a logistic regression model! Let's do that. (See also the excellent visualization from the Statistical Sleuth text, in the slides.)

```{r}
fit <- glm(survived~age, family=binomial(link="logit"), data=titanic)
ages <- data.frame(age=titanic$age)
predictions <- predict(fit, newdata=ages, type="response", se=TRUE)

ggdata <- data.frame(age = ages$age,
                     pred = predictions$fit,
                     LL = predictions$fit-1.96*predictions$se.fit,
                     UL = predictions$fit+1.96*predictions$se.fit
                    )

ggplot(ggdata, aes(x = age, y = pred)) + 
   geom_ribbon(aes(ymin = LL, ymax = UL), alpha = .2) +
   geom_line(colour="blue", size=1)
```



We kind of already did, but let's now run the regression

```{r}
fit <- glm(survived~age, family=binomial(link="logit"), data=titanic)
summary(fit)
```

Let's first estimate the proportion of survivors in a very roundabout way:

```{r}
fit <- glm(survived~1, family=binomial(link="logit"), data=titanic)
summary(fit)
plogis(fit$coefficients["(Intercept)"])
```

Let's confirm that:
```{r}
mean(titanic$survived)
logit(mean(titanic$survived))
```


Predicting with age:
```{r}
fit <- glm(survived~age, family=binomial(link="logit"), data=titanic)
summary(fit)
```


Interetation of $\beta_0$: the predicted probability of survival for people aged 0.

We already did this when graphing the probabilities, but here it is again, in more detail


```{r}
fit <- glm(survived~age, family=binomial(link="logit"), data=titanic)
ages <- data.frame(age=titanic$age)
logodds <- predict(fit, newdata=ages, type="link", se=TRUE)
logodds$fit[1:10]
logodds$se.fit[1:10]
```

Here, we are getting the log-odds, together with their SEs. We can transform the implied interval to obtain the uncertainty for the probabilities

```{r}
plogis(logodds$fit[1])
c(plogis(logodds$fit[1]-1.96*logodds$se.fit[1]), plogis(logodds$fit[1]+1.96*logodds$se.fit[1]))
```

Note that this interval is not symmetric. Why?

(This is not really a CI, since the probability is never observed.)


Let's now visualize the difference between men and women, controlling for age

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

newdata <- data.frame(age=0:100, sex=rep(factor("male"), 101))
predictions_male <- predict(fit, newdata=newdata, type="response", se=TRUE)

newdata <- data.frame(age=0:100, sex=rep(factor("female"), 101))
predictions_female <- predict(fit, newdata=newdata, type="response", se=TRUE)


ggdata <- data.frame(age = newdata$age,
                     pred_male = predictions_male$fit,
                     mLL = predictions_male$fit-1.96*predictions_male$se.fit,
                     mUL = predictions_male$fit+1.96*predictions_male$se.fit,
                     pred_female = predictions_female$fit,
                     fLL = predictions_female$fit-1.96*predictions_female$se.fit,
                     fUL = predictions_female$fit+1.96*predictions_female$se.fit
                    )

ggplot(ggdata, aes(x = age)) + 
   geom_ribbon(aes(ymin = mLL, ymax = mUL, fill="male"), alpha = .2) + guides(fill=FALSE) +
   geom_line(aes(y = pred_male, color="male"), size=1) +
   geom_ribbon(aes(ymin = fLL, ymax = fUL, fill="female"), alpha = .2) + guides(fill=FALSE) +
   geom_line(aes(y = pred_female, color="female"), size=1) +
   labs(y = "Survival probs")


```


Also see a beautiful visualization of the Titanic data [here](https://plot.ly/ggplot2/geom_ribbon/).



### Likelihood Ratio Tests


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