Lecture 13

class: center, middle, inverse, title-slide

# Lecture 13
## Logistic regression
### Jung-Jin Lee
### Apr 21, 2020

---

## Odds

As a measure of likelihood of an event, the **odds** of an event are the ratio of the probability that the event occurs to the probability that the event cannot occur. In other words, for an event `$A$`, the odds of `$A$`, denoted `$\text{odds}(A)$`, is defined to be 
`$$\text{odds}(A) = \frac{P(A)}{1-P(A)}.$$`

**Example**: when rolling a fair die, the odds of getting 1 are `$\frac{1/6}{5/6} = \frac{1}{5}$` (or 1:5). The odds of getting an even number are `$\frac{1/2}{1/2} = 1$` (or 1:1).

It is notable that `$P(A)$` can be recovered from `$\text{odds}(A)$`:
`$$P(A)=\frac{\text{odds}(A)}{1+\text{odds}(A)}.$$`
---

## Titanic dataset

.pull-left-code[
- On April 15, 1912, the largest passenger liner ever made collided with an iceberg during her maiden voyage.

- When the Titanic sank it killed 1502 out of 2224 passengers and crew. 
  - One of the reasons that the shipwreck resulted in such loss of life was that there were not enough lifeboats for the passengers and crew.
  
  - Although there was some element of luck involved in surviving the sinking, some groups of people were more likely to survive than others.
]
.pull-right-plot[
<p align = "center">
  <img src = Titanic.png width = "240">
</p>  
]

---

## Titanic dataset

Download a comma separated file [titanic.csv](titanic.csv). It consists of 887 observations and 8 variables. You can read in the dataset using `read_csv()`, a function that comes with `tidyverse`.

```r
*titanic_ori <- read_csv(file = "titanic.csv")
head(titanic_ori)
```

```
## # A tibble: 6 x 8
##   Survived Pclass Name  Sex     Age `Siblings/Spous… `Parents/Childr…
##      <dbl>  <dbl> <chr> <chr> <dbl>            <dbl>            <dbl>
## 1        0      3 Mr. … male     22                1                0
## 2        1      1 Mrs.… fema…    38                1                0
## 3        1      3 Miss… fema…    26                0                0
## 4        1      1 Mrs.… fema…    35                1                0
## 5        0      3 Mr. … male     35                0                0
## 6        0      3 Mr. … male     27                0                0
## # … with 1 more variable: Fare <dbl>
```

---

## Titanic dataset

```r
names(titanic_ori)
```

```
## [1] "Survived"                "Pclass"                 
## [3] "Name"                    "Sex"                    
## [5] "Age"                     "Siblings/Spouses Aboard"
## [7] "Parents/Children Aboard" "Fare"
```

- `Survived`: did the passenger survive? (1 = survived)

- `Pclass`: passenger-class

- `Siblings/Spouses Aboard`: number of siblings/spouses aboard

- `Parents/Children Aboard`: number of parents/children aboard

- `Fare`: fare in UK pounds

```r
titanic <- titanic_ori %>%
  mutate(Pclass = factor(Pclass, levels = c(1, 2, 3))) %>%
  mutate(Sex = factor(Sex, levels = c("female", "male")))
```

---

## Titanic dataset

`Age` distribution by `Sex` and `Pclass`:

.pull-left-40[

```r
titanic %>% 
  ggplot(aes(Sex, Age)) + 
  geom_boxplot() + 
  facet_wrap(~Pclass) 
```
]

.pull-right-60[
<img src="Drexel_2020_Lecture_13_files/figure-html/unnamed-chunk-6-1.png" style="display: block; margin: auto;" />
]

---

## Logistic regression

How does `Age` affect the `Survived`?

.pull-left-55[

```r
fit.lm <- lm(Survived ~ Age, 
             data = titanic)
beta0 <- fit.lm$coefficients[1]
beta1 <- fit.lm$coefficients[2]
ggplot(titanic, aes(Age, Survived)) +
  geom_point() +
  geom_abline(intercept = beta0, 
              slope = beta1, 
              color = "red")
```
]

.pull-right-45[
<img src="Drexel_2020_Lecture_13_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" />
]

A simple linear regression does not work properly -- for example, what can be said about the survival of a 40-year-old passenger?

---

## Logistic regression

In logistic regression, we model 
`$$\log(\text{Odds of survival}) = \beta_0 + \beta_1\texttt{Age},$$`
--
that is, 
`$$\text{Odds of survival} = e^{\beta_0 + \beta_1\texttt{Age}},$$`
or equivalently,
`$$\frac{P(\text{Survival})}{1-P(\text{Survival})} = e^{\beta_0 + \beta_1\texttt{Age}}.$$`
--
This can be rewritten as 
`$$P(\text{Survival}) = \frac{e^{\beta_0 + \beta_1\texttt{Age}}}{1+e^{\beta_0 + \beta_1\texttt{Age}}}.$$`
--
The coefficients `$\beta_0$` and `$\beta_1$` are estimated by the function `glm()`.

---

## Logistic regression

```r
*fit.glm <- glm(Survived ~ Age, family = binomial(link = "logit"),
*              data = titanic)
fit.glm$coefficients
```

```
##  (Intercept)          Age 
## -0.209188757 -0.008774355
```

From this output, we conclude that

`$$\log(\text{Odds of survival}) = -0.209188757-0.008774355 \times \texttt{Age}$$`

For a 40-year-old passenger, the odds of survival is
`$$e^{-0.209188757-0.008774355 \times 40} = 0.571116.$$`
In other words, the probability of survival is
`$$\frac{e^{-0.209188757-0.008774355 \times 40}}{1+e^{-0.209188757-0.008774355 \times 40}}=36.35\%.$$`

---

## Interpretation of coefficients

Interpretation of `$\beta_1$` (-0.008774355 in the previous example):
--

- Odds of survival at age `$x$`:
`$$\text{Odds}(x) = e^{\beta_0 +\beta_1x}.$$`
--

- Odds of survival at age `$x+1$` (a unit increase in age):
`$$\text{Odds}(x+1) = e^{\beta_0 +\beta_1(x+1)}.$$`
--

- In other words, 
`$$\text{Odds}(x+1) = e^{\beta_0 +\beta_1(x+1)} = e^{\beta_0 +\beta_1x+\beta_1}=e^{\beta_1}\times\text{Odds}(x),\quad\text{or}$$`
`$$\frac{\text{Odds}(x+1)}{\text{Odds}(x)}=e^{\beta_1}.$$`
--

In general, `$e^{\beta_1}$` represents the ratio of two odds when there is a unit change in the predictor variable.

---

## Interpretation of coefficients

`$$\frac{\text{Odds}(x+1)}{\text{Odds}(x)}=e^{\beta_1}.$$`

- `$\beta_1 > 0 \Leftrightarrow e^{\beta_1} > 1 \Leftrightarrow \text{Odds}(x+1) > \text{Odds}(x) \Leftrightarrow$` increase in the predictor leads to increase in the odds.

- `$\beta_1 < 0 \Leftrightarrow e^{\beta_1} < 1 \Leftrightarrow \text{Odds}(x+1) < \text{Odds}(x) \Leftrightarrow$` increase in the predictor leads to decrease in the odds.

---

## Logistic regression with a categorical predictor

How does `Sex` affect the `Survived`?

```r
fit.glm2 <- glm(Survived ~ Sex, family = binomial(link = "logit"),
                data = titanic)
fit.glm2$coefficients
```

```
## (Intercept)     Sexmale 
##    1.056589   -2.505126
```

Interpretation:

- The odds of survival of female = `$e^{\beta_0} = e^{1.056589} = 2.877$`. This implies that the probability of survival of a female passenger equals `$\frac{e^{1.056589}}{1+e^{1.056589}}=74.2\%$`.

- On the other hand, the odds of survival of male = `$e^{\beta_0 + \beta_1} = e^{1.056589-2.505126} = 0.235$`. This implies that the probability of survival of a female passenger equals `$\frac{e^{1.056589-2.505126}}{1+e^{1.056589-2.505126}}=19.0\%$`.

---

## Logistic regression with a categorical predictor

```r
titanic %>%
  mutate(Survival = ifelse(Survived == 1, "yes", "no")) %>%
  ggplot(aes(Sex, fill = Survival)) +
  geom_bar()
```

---

## Logistic regression with a categorical predictor

```r
table(titanic$Survived, titanic$Sex)
```

```
##    
##     female male
##   0     81  464
##   1    233  109
```

```r
female_survival_prob <- 233/(81+233)
male_survival_prob <- 109/(464+109)
female_survival_prob; male_survival_prob
```

```
## [1] 0.7420382
```

```
## [1] 0.1902269
```

---

## Logistic regression with multiple predictors

How do `Age` and `Pclass` affect the `Survived`?

```r
fit.glm3 <- glm(Survived ~ Age + Pclass, 
                family = binomial(link = "logit"),
                data = titanic)
fit.glm3$coefficients
```

```
## (Intercept)         Age     Pclass2     Pclass3 
##  2.09929535 -0.03946933 -1.03811692 -2.28986274
```

Interpretation:

- When `Pclass` remains unchanged, the odds of survival of a patient with age `$x+1$` is `$e^{-0.03946933}=0.9612994$` times the odds of survival of a patient with age `$x$` (The older, the less likely to survive).

- When `Age` remains unchanged, the odds of survival of a patient in passenger-class `$2$` is `$e^{-1.03811692}=0.3541209$` times the odds of survival of a patient in passenger-class `$1$`.

- When `Age` remains unchanged, the odds of survival of a patient in passenger-class `$3$` is `$e^{-2.28986274}=0.1012804$` times the odds of survival of a patient in passenger-class `$1$`.

---

## Logistic regression with multiple predictors

```r
titanic %>%
  mutate(Survival = ifelse(Survived == 1, "yes", "no")) %>%
  ggplot(aes(Pclass, fill = Survival)) +
  geom_bar()
```

---

## Exercise

Recall the horseshoe crabs dataset. We want to model the existence of satellites of a female crab using the weight and spine condition of the  female crab.

```r
crabs <- read_tsv("crabs.tsv") %>%
  mutate(has.satell = ifelse(satell > 0, 1, 0)) %>%
  mutate(spine = factor(spine, levels = c("1", "2", "3"))) %>%
  mutate(weight = weight/1000) # weight in grams to kilograms
```
--

- What is the probability that a female crab with weight 4.1kg and spine condition 3 has at least one satellite?
- What is the probability that a female crab with weight 2.8kg and spine condition 1 has at least one satellite?

---

## Exercise

```r
fit.crab <- glm(has.satell ~ weight + spine,
                family = binomial(link = "logit"),
                data = crabs)
fit.crab$coefficients
```

```
## (Intercept)      weight      spine2      spine3 
## -3.65268080  1.78814346 -0.29967472  0.07082751
```
--

For female crabs with spine condition 3, the odds of having a satellite equals 
`$$e^{-3.65268080 + 1.78814346\times\texttt{weight}+0.07082751}.$$`
--

In particular if the weight is `$4.1$`, then the odds become
`$$e^{-3.65268080 + 1.78814346\times4.1+0.07082751}=42.50.$$`
--

This implies that the probability of having a satellite is 
`$$\frac{e^{-3.65268080 + 1.78814346\times4.1+0.07082751}}{1+e^{-3.65268080 + 1.78814346\times4.1+0.07082751}}=97.7\%.$$`

---

## Exercise

For female crabs with spine condition 1, the odds of having a satellite equals 
`$$e^{-3.65268080 + 1.78814346\times\texttt{weight}}.$$`
--

In particular if the weight is `$2.8$`, then the odds become
`$$e^{-3.65268080 + 1.78814346\times2.8}=3.87.$$`
--

This implies that the probability of having a satellite is 
`$$\frac{e^{-3.65268080 + 1.78814346\times2.8}}{1+e^{-3.65268080 + 1.78814346\times2.8}}=79.5\%.$$`

---

## Exercise

```r
crabs %>%
  mutate(has.satell = ifelse(has.satell == 1, "yes", "no")) %>%
  ggplot(aes(spine, fill = has.satell)) +
  geom_bar()
```

---

## Exercise

```r
crabs %>%
  mutate(has.satell = ifelse(has.satell == 1, "yes", "no")) %>%
  ggplot(aes(has.satell, weight)) +
  geom_boxplot()
```

---

## Final exam

- Tuesday, April 28. 5:30pm ~ 10:00pm

- Online (PDF/Word)

- Cumulative (more than 70% from material after Exam 2)

```r
num1 <- c(20, 8, 1, 14, 11)
num2 <- c(25, 15, 21)
char1 <- LETTERS[num1]
char2 <- LETTERS[num2]
IWannaSay <- paste(paste(char1, collapse = ""),
                   paste(char2, collapse = ""))
print(IWannaSay)
```

```
## [1] "THANK YOU"
```

```