---
title: "Probability, Odds, and Simple Logistic Regression"
subtitle: "Notes and in-class exercises"
format: 
  html:
    embed-resources: true
    toc: true
---


You can download the .qmd file for this activity [here](../activity_templates/15+16-Logistic-Regression-I.qmd) and open in R-studio. The rendered version is posted in the [course website](https://mutasim221b.github.io/Mac-STAT-155-Sp-26/) (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins.



```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')
```





# Learning goals

By the end of this lesson, you should be able to:

- Distinguish between probabilities and odds, and convert one to the other
- Make appropriate visualizations for displaying relationships between multiple categorical variables (mosaic plots, stacked bar plots, etc.)
- Explain the differences between linear regression and logistic regression for modeling binary outcomes
- Construct simple logistic regression models in R
- Interpret coefficients in simple logistic regression models
- Use simple logistic regression models to make predictions
- Describe the form (shape) of relationships on the log odds, odds, and probability scales


## Readings and videos

Go through the following reading or videos **before** class:

- Reading: Sections 2.5, Section 6.2 introduction, and Sections 4.1-4.3 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)

- Videos:
  - [Prob vs. Odds vs. Log Odds](https://www.youtube.com/watch?v=FgGGYv3oO6s) ([script](https://drive.google.com/file/d/1QLl1-6T3jEO_HkeF4Qo4-qIdI5cR7yoy/view?usp=sharing))
  - [Calculating Probability and Odds from 2x2 Tables](https://youtu.be/LsLK8C-ORbk) ([script](https://drive.google.com/file/d/158QHXFZFQBkSgcqNS-FC0WMQ2kmvAq1k/view?usp=sharing))
  - [Logistic regression](https://youtu.be/yi7I2vmxbiQ) ([slides](https://drive.google.com/file/d/1o5WwSnShEazkV_HTWtOCeK4W0snrsep3/view?usp=sharing))


**File organization:** Save this file in the "Activities" subfolder of your "STAT155" folder.



# Live Note-taking!
- Try this [Coin-toss Applet](https://www.rossmanchance.com/applets/2021/oneprop/OneProp.htm).



# Exercises on Probability and Odds 

**Context:** To begin formally learning about probabilities and odds, we’ll be exploring a dataset containing information on 2,500 singleton (i.e. not twins) births in King County, Washington in 2001. Each row contains information from one birth parent, and there are no birth parents included in the dataset more than once. 

The main research question this study aimed to answer was whether the [First Steps](https://kingcounty.gov/en/dept/dph/health-safety/health-centers-programs-services/maternity-support-wic/maternity-support-services-infant-case-management) program in King County improved birth outcomes for women from socioeconomically disadvantaged backgrounds. We'll attempt to answer this research question using the information available to us as we go!

The variables in this dataset we'll look at more closely for each birth parent are:

- `age`: age of birth parent at time of birth (years)
- `parity`: number of children the birth parent has given birth to before
- `married`: indicator for whether the birth parent is currently married (1 = yes, 0 = no)
- `bwt`: birthweight of the child (in grams)
- `smokeN`: number of cigarettes smoked per day during pregnancy
- `drinkN`: number of alcoholic drinks consumed per day during pregnancy
- `firstep`: indicator for whether the birth parent participated in the “First Steps” pregnancy program
- `gestation`: number of weeks at which birth parent gave birth

Run the code below to read in the `firststeps` data, and create a few new variables that we'll explore as well.


```{r}
library(tidyverse)
firststeps <- read_csv("https://mac-stat.github.io/data/firststeps.csv") %>%
  mutate(firstchild = ifelse(parity == 0, "Yes", "No"), # Is this the first child this parent has had?
         low_bwt = ifelse(bwt < 2500, "low", "not low"),
         preterm = ifelse(gestation < 37, "Yes", "No")) # short gestational period
```




## Exercise 1: Exploring First Steps enrollment and Gestational Age

A baby born prior to 37 weeks is considered premature. In figuring out whether we have evidence that the First Steps program is associated with better birth outcomes than those not in the First Steps program, we can look at whether the individuals in the program are more likely to have preterm babies.

Below, we make a 2x2 table in R:

```{r}
# 2x2 Table: preterm vs. First Steps
firststeps %>% 
    count(preterm, firstep)
```

You may be wondering why this is called a 2x2 table, when it looks as though the table has four rows and three columns. The data can be re-arranged (and usually is, in a formal report) as follows...

```{r}
table(firststeps$preterm, firststeps$firstep)
```

... but it's much cleaner to code up the original way! 

a. How many birth parents were enrolled in the First Steps program? Which rows did you use to calculate this number?

> 343 + 60 = 403 parents were enrolled in the First Steps program; used both rows of the table where `firstep = 1`.

b. What *percentage* of people in the study were enrolled in the First Steps program? **Recall: there were 2500 participants! You can confirm this by adding up the entire third column of the table**

> 16.12% of people in the study were enrolled in First Steps!

```{r eval = TRUE}
403 / 2500
```

c. How many birth parents *who were enrolled in First Steps* had a premature baby?

> 60 birth parents

d. What *percentage* of birth parents in First Steps had a premature baby? Think carefully about the numerator and denominator you use to calculate this!

> 14.89% of birth parents in First Steps had a premature baby.

```{r eval = TRUE}
60 / 403
```

e. What *percentage* of birth parents who had a premature baby were enrolled in First Steps? Think carefully about the numerator and denominator you use to calculate this! 

> The total number of birth parents who had a premature baby was 218 + 60 = 278. Of those. 60 were enrolled in First Steps. Therefore, 21.58% of birth parents who had a premature baby were enrolled in First Steps.

```{r eval = TRUE}
60/278
```



Congratulations! If you've made it to this point, you already know from the class notes what *marginal* and *conditional* probabilities are. Formally, 

- a *marginal* probability, denoted $P(A)$ for an event $A$, is the probability that $A$ occurs overall. You calculated the marginal probability that people were enrolled in First Steps in part (b)! In this case, the denominator used to calculate the probability was the total number of people in the study.

- a *conditional* probability, denoted $P(A | B)$ for events $A$ and $B$, is the probability that $A$ occurs **given** that event $B$ occurs. You calculated the conditional probability that a premature baby was born **given** that a parent was in First Steps in part (d)! In this case, the denominator used to calculate the probability was the total number of birth parents in the First Steps program. You also calculated a conditional probability in part (e).

Using formal probability notation, write the probabilities you calculated in parts (b), (d), and (e) as

$$
b. P(\text{First Steps}) = .1612
$$

$$
d. P(\text{Preterm} | \text{First Steps}) = .1488
$$

$$
e. P(\text{First Steps} | \text{Preterm}) = .2158
$$


Note that the conditional probabilities calculated in parts (d) and (e) are not the same! This is because which event you condition on alters the denominator, and the event you're interested in alters the numerator.

f. To determine if gestational age differed by enrollment in First Steps, we'll want to calculate the conditional probability that a baby is born prematurely given First Steps enrollment (Done!.1488), *and* given that a parent is *not* enrolled in First Steps. Use the 2x2 table to calculated this conditional probability.

```{r}
table(firststeps$preterm, firststeps$firstep)
```

$$
P(\text{Preterm} | \text{Not in First Steps}) = 218 / (218 + 1879) = 0.103958
$$

g. A ratio of conditional probabilities, where the conditioning event is the same for both, tells us how many times *more likely* an event is to occur for one group compared to another. Calculate how many times more likely a birth parent enrolled in First Steps is to have a premature baby compared to birth parents not enrolled in First Steps.

$$
\frac{(\text{Preterm} | \text{First Steps})}{P(\text{Preterm} | \text{Not in First Steps})} = .1488 / 0.103958 = 1.43
$$

h. Write a two-sentence summary, appropriate for a general audience, summarizing your results in terms of a *ratio of probabilities*. Does gestational age appear to differ greatly by First Steps enrollment? What does this imply about the effectiveness of the First Steps program, if anything?

> Parents in this study in the First Steps program are 1.43 times more likely to have a premature birth than those not enrolled in the First Steps program, indicating that gestational age does differ by First Steps enrollment. This implies that enrollment in the First Steps program may not be associated with better birth outcomes, as measured by gestational age.



i. To go along with your summary, let's make a visualization! There are three basic options for visualization two categorical variables. All are perfectly valid, but some may be more useful to read than others, and display different information.

You'll see one other fancier option (called a `mosaic` plot) in the next activity. 

```{r}
# Side-by-side bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar(position = "dodge") +
  theme_classic()

# Stacked bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar() +
  theme_classic()

# Stacked relative frequency bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar(position = "fill") +
  theme_classic()
```

**Bonus Question:** Which of the above three plots allows you to directly see the conditional probabilities we calculated previously?

> The stacked relative frequency bar chart! 



## Exercise 2: Exploring First Steps enrollment and Low birthweights

Another birth outcome we can consider when comparing those enrolled in the First Steps program to those not enrolled is birth weight. A baby is considered to have *low* birth weight when birth weight is less than 2500 grams.


a. Fill in the code below to make a table comparing `low_bwt` to `firsteps`.

```{r}
# 2x2 Table: low_bwt vs. First Steps
firststeps %>%
  count(low_bwt, firstep)
table(firststeps$low_bwt, firststeps$firstep)
```
b. Using the table from part (a), calculate the following conditional probabilities:

$$
P(\text{Low birth weight} | \text{First Steps}) = 25 / (25 + 378) = 0.062
$$

$$
P(\text{Normal birth weight} | \text{First Steps}) = 378 / (25 + 378) = 0.938
$$


$$
P(\text{Low birth weight} | \text{Not in First Steps}) = 102 / (102 + 1995) = 0.049
$$


$$
P(\text{Normal birth weight} | \text{Not in First Steps}) = 1995 / (102 + 1995) = 0.951
$$


An additional numerical summary that is often useful when working with indicator variables is *odds*. Odds are defined as

$$
Odds = \frac{p}{1 - p}
$$

where $p$ is the probability that an event occurs. Therefore, if we know $p$, we can calculate the odds that an event happens! Similarly, if we know the odds, we can calculate $p$ using

$$
p = Odds / (1 + Odds)
$$


We can also calculate odds from our 2x2 (or 3x2, 4x2, ...) tables. In colloquial terms, probabilities are "yes"'s over "total"'s, and odds are "yes"'s over "no's". In pseudo-math:

$$
p = \frac{Yes}{Total}, \quad Odds = \frac{Yes}{No}
$$
We'll see why odds are especially useful when we have binary outcome variables in a regression model in the next activity. For now, note that they're also commonly used in lots of contexts: sports, gambling, case-control studies, etc.

c. Using your answer to part (b), calculate the following *odds* 

$$
Odds(\text{Low birth weight} | \text{First Steps}) = 0.062 / (1 - 0.062) = 0.06609808
$$

$$
Odds(\text{Normal birth weight} | \text{First Steps}) = 0.938 / (1 - 0.938) = 15.12903
$$

$$
Odds(\text{Low birth weight} | \text{Not in First Steps}) = 0.049 / (1 - 0.049) = 0.05152471
$$
$$
Odds(\text{Normal birth weight} | \text{Not in First Steps}) = 0.951 / (1 - 0.951) = 19.40816
$$

d. A ratio of odds (called an **odds ratio**, unsurprisingly) tells us how many times *higher* or *greater* the odds are that an event occurs, comparing one group to another. This might sound irritatingly circular. The key here is that while odds ratios *do* allow us to compare binary/indicator outcomes from one group to one another, they *do not* tell us how much *more likely* an event is to occur comparing those same groups. This is distinct from ratios of probabilities!

Calculate the ratio of the odds of having a low-birth-weight baby, comparing those in the First Steps program to those not in the First Steps program (i.e., how many times higher/lower is the odds of having a low-birth-weight baby among those in First Steps as compared to those not in First Steps?)


```{r eval = TRUE}
0.06609808 / 0.05152471
```



e. Write a two-sentence summary, appropriate for a general audience, summarizing your results in terms of an *odds ratio*. Does birth weight appear to differ greatly by First Steps enrollment? What does this imply about the effectiveness of the First Steps program, if anything?

> The odds of having a low birth weight baby are 1.28 times higher for those enrollment in First Steps compared to those not in First Steps. Just as in Exercise 1, this implies that the First Steps program may not be associated with improved birth outcomes (with the same caveats as given in the answer to 1 (h)).

f. To go along with your summary, add code below to make one of the three visualization options we tried out in Exercise 1.

```{r}
# Stacked relative frequency bar chart (with some fancy aesthetics)
firststeps %>%
  mutate(Birthweight = low_bwt %>% str_to_title()) %>%
  ggplot(aes(firstep, fill = Birthweight)) +
  geom_bar(position = "fill") +
  theme_classic() +
  scale_fill_viridis_d(option = "H") +
  labs(x = "First Steps", title = "Birthweight by First Steps Enrollment") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Enrolled", "Enrolled")) 

```




# Probability and Odds: Additional Practice 

## Exercise 3: Conditional vs. Marginal probabilities

Suppose we select a person at random from the entire global population. For each of the following probabilities, which do you think is bigger? Explain your reasoning.

a. P(lung cancer) or P(lung cancer | smoker)

> Your response here

b. P(likes McDonald's) or P(likes McDonald's | vegetarian)

> Your response here

c. P(smart | Mac grad) or P(Mac grad | smart)

> Your response here

## Exercise 4: Probability practice

Let's explore whether birthweight of a baby varies by whether or not it was the *first* child that a mother had, *and* whether this relationship differs by First Steps enrollment. We make a table below:

```{r}
firststeps %>%
  count(firstchild, low_bwt, firstep)
```

a. What is the probability that a mother enrolled in First steps who is having their first child, has a baby who is born at a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(___ | ___) = ?

b. What is the probability that a mother *not* enrolled in First steps who is having their first child, has a baby who is born at a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(___ | ___) = ?

c. What is the probability that a mother's first child has a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(___ | ___) = ?

d. How many times *more likely* is a child to be born at a low birthweight, comparing children who are the first born to those not first born?

> P(___ | ___) / P(___ | ___) = ?



\
\
\
\

# Exercises on Logistic Univariate

**Context:** The *Titanic* was a British passenger ship that famously sank in 1912 after hitting an iceberg in the North Atlantic Ocean. Approximately 2200 passengers were on board the Titanic, and it's estimated that 1500 of them did not survive the crash. Historians have worked diligently to collect data on the passengers that were aboard the Titanic. 

We have data for 1313 passengers, where the following information is available for each passenger:

- `Name`: name
- `PClass`: ticket class (1st, 2nd, 3rd)
- `Age`: age (years)
- `Sex`: binary sex (female, male)
- `Survived`: indicator that the passenger survived (1 = survived, 0 = died)

Our question of interest is: *how do different factors relate to survival?*

In the Console, run `install.packages("ggmosaic")` to install the `ggmosaic` package that we'll be using to make a specialized type of plot.

```{r}
library(tidyverse)
library(ggmosaic)

titanic <- read_csv("https://mac-stat.github.io/data/titanic.csv")
head(titanic)
```

## Exercise 1: Exploring age

Did younger passengers tend to have higher survival rates than older passengers?

Visualizing the relationship between a binary response and a quantitative predictor can be tricky. We will take a few approaches here.

1. Create a boxplot where one box corresponds to the age distribution of survivors and the second to that of non-survivors.
2. Create density plots with separate colors for the survivors and non-survivors.
3. The remainder of the code below creates a plot of the fraction who survived at each age. (Since we have a large data set and multiple (though sometimes not many) observations at most ages, we can manually calculate the survival fraction.

After inspecting the plots, summarize what you learn.

```{r}
# Create a boxplot
# Note that you'll need to force R to view Survived as a binary categorical variable by using x = factor(Survived) instead of just x = Survived in the aes() part of your plot

ggplot(titanic, aes(x = factor(Survived), y = Age)) +
    geom_boxplot()

# Can flip the boxplot on its side too
ggplot(titanic, aes(y = factor(Survived), x = Age)) +
    geom_boxplot()


# Create a density plot (you'll need to use factor(Survived) again)
ggplot(titanic, aes(x = Age, color = factor(Survived))) +
    geom_density()

# Use the code below to create a plot of the fraction who survived at each age
titanic_summ <- titanic %>% 
    group_by(Age) %>%
    summarize(frac_survived = mean(Survived))

ggplot(titanic_summ, aes(x = Age, y = frac_survived)) +
    geom_point() +
    geom_smooth(se = FALSE)
```
> The boxplot doesn’t clearly indicate a difference in the age distributions across survivors and non-survivors, but we do notice from the density plot that there is a greater density of younger passengers among the survivors. We also see from the last plot that younger passengers tend to have a higher survival chance.


## Exercise 2: Exploring sex and ticket class

Were males or females more likely to survive? Did 1st class passengers tend to survive more than 2nd and 3rd class passengers?

The code below creates plots that allow us to explore how `Sex` and `PClass` relate to survival. The first two plots are standard bar plots that use color to indicate what fraction of each group survived. The last two plots are mosaic plots that are much like the standard bar plots, but the width of the bars reflects the distribution of the x-axis variable. (The widest bar is the most prevalent category.)

Summarize what you learn about the relationship between sex, ticket class, and survival.

```{r}
# Standard bar plots
ggplot(titanic, aes(x = Sex, fill = factor(Survived))) +
    geom_bar(position = "fill")

ggplot(titanic, aes(x = PClass, fill = factor(Survived))) +
    geom_bar(position = "fill")

```
> Females were more likely to survive than males.
1st class was most likely to survive, followed by 2nd then 3rd class.


## Exercise 3: Linear regression model

For now we will focus on exploring the relationship between (ticket) class and survival. 

Let's tabulate survival across classes. We can tabulate across two variables by providing both variables to `count()`:

```{r}
titanic %>% 
    count(PClass, Survived)
```

a. Using the `count()` output, we can fill in the following contingency table:

Class     | Died     | Survived | Total 
----------|----------|----------|-------
1st Class | 129      |  193     | 322
2nd Class | 160      |  119     | 279
3rd Class | 573      |  138     | 711
Total     | 862      |  450     | 1312

b. Using the table, we can estimate the followings:
    - the probability of surviving among 1st class passengers: 193/322 = 0.599
    - the probability of surviving among 2nd class passengers: 119/279 = 0.427
    - the probability of surviving among 3rd class passengers: 138/711 = 0.194
    - the difference in the probability of surviving, comparing 2nd class passengers to 1st class passengers (i.e., how much lower is the probability of 2nd class passengers as compared to 1st class passengers?): 119/279 - 193/322 = -0.173
    - the difference in the probability of surviving, comparing 3rd class passengers to 1st class passengers  (i.e., how much lower is the probability of 3rd class passengers as compared to 1st class passengers?): 138/711 - 193/322 = -0.405


c. After fitting the linear regression model below, write out the model formula using correct notation. Explain carefully what it means to talk about the expected/average value of a binary variable.

```{r}
lin_mod <- lm(Survived ~ PClass, data = titanic)
summary(lin_mod)
```

>  This model can be written as: $E[Survived | PClass] = \beta_0 + \beta_1 PClass2nd + \beta_2 PClass3rd$.
    - In the context of a binary variable, the expected value/average is the same as the probability that the variable equals one. To see an example of this, calculate the average of this list of 0's and 1's: (0,0,1,1,0,1,0,1). Now calculate the proportion of 1's. What do you notice?
    - This means that we can also write this model as follows: $P[Survived = 1 | PClass] = \beta_0 + \beta_1 PClass2nd + \beta_2 PClass3rd$




d. Write an interpretation of each of the coefficients in your linear regression model. How do your coefficient estimates compare to your answers in part b?

> The coefficient estimates are the differences in probability from part b!
    - `(Intercept)`: the estimated probability of survival for passengers in 1st class is 0.599 (59.9%)
    - `PClass2nd`: the difference in the estimated probability of survival comparing passengers in 1st class to passengers in 2nd class is 0.173 (17.3%), where passengers in 1st class have the higher estimated survival probability
        - OR... comparing passengers in 1st class to passengers in 2nd class, the difference in the proportion of passengers that survived is 0.173 (17.3%), with 1st class having a higher proportion of passengers that survived
        - OR... the probability of survival is 17.3% lower among passengers in 2nd class than it is among passengers in 1st class
    - `PClass3rd`: the difference in the estimated probability of survival comparing passengers in 1st class to passengers in 3rd class is 0.405 (40.5%), where passengers in 1st class have the higher estimated survival probability


> Live Note Taking! (Logistic Model)


## Exercise 4: Logistic regression model (categorical predictor) 

a. Refer back to your contingency table from Exercise 3a. Using your table, estimate the following:
    - the odds of surviving among 1st class passengers: 193/129 = 1.496
    - the odds of surviving among 2nd class passengers: 119/160 = 0.744
    - the odds of surviving among 3rd class passengers
    - the ratio of the odds of surviving, comparing 2nd class passengers to 1st class passengers (i.e., how many times higher/lower is the odds of survival among 2nd class passengers as compared to 1st class passengers?): 0.744/1.496 = 0.497
    - the ratio of the odds of surviving, comparing 3rd class passengers to 1st class passengers
    
Hints: Remember that we can also calculate odds from our 2x2 (or 3x2, 4x2, ...) tables. In colloquial terms, probabilities are "yes"'s over "total"'s, and odds are "yes"'s over "no's". In pseudo-math:

$$
p = \frac{Yes}{Total}, \quad Odds = \frac{Yes}{No}
$$


Class     | Died     | Survived | Total 
----------|----------|----------|-------
1st Class | 129      |  193     | 322
2nd Class | 160      |  119     | 279
3rd Class | 573      |  138     | 711
Total     | 862      |  450     | 1312

> Place your answer here, create R code chunk as needed

b. After fitting the logistic regression model below, write out the model formula using correct notation.


```{r}
log_mod <- glm(Survived ~ PClass, data = titanic, family = "binomial")
coef(summary(log_mod))
```

```{r eval = TRUE}
# Calculations for exponentiating coefficients
exp(0.4028778)
exp(-0.6989281)
exp(-1.8265098)
```


d. Write an interpretation of each of the *exponentiated* coefficients in your logistic regression model. Think carefully about what we are modeling when we fit a logistic regression model. How do these exponentiated coefficient estimates compare to your answers in part a?

> Place your answer here



## Exercise 5: Logistic regression model (quantitative predictor)

Now we will explore how to interpret a *quantitative* predictor in a logistic regression model.

a. After fitting the logistic regression model below, write out the model formula using correct notation.

```{r}
log_mod <- glm(Survived ~ Age, data = titanic, family = "binomial")
coef(summary(log_mod))
```

$\log(Odds[Survived = 1 | Age]) = \beta_0 + \beta_1 Age$

b. Write an interpretation of each of the *exponentiated* coefficients in this logistic regression model.

- exp(Intercept): $exp(-0.0814)=0.92$ --> the estimated odds of survival among passengers who are 0 years old is 0.92 (i.e., passengers who are 0 years old are 0.92 times more likely to survive than they are to die--so very slightly *more* likely to die)
- exp(Age): $exp(-0.0088)=0.99$ --> For every 1-year increase in a passenger's age, the estimated odds of survival *decrease* by about 1%.





## Exercise 6: Linear vs. logistic modeling

To highlight a key difference between linear vs. logistic modeling, consider the following linear and logistic regression models of survival with sex and age as predictors in addition to ticket class.

```{r}
lin_mod2 <- lm(Survived ~ PClass + Sex + Age, data = titanic)
coef(summary(lin_mod2))

log_mod2 <- glm(Survived ~ PClass + Sex + Age, data = titanic, family = "binomial")
coef(summary(log_mod2))
```

a. Use the linear regression model to predict the probability of survival for Rose (a 17 year old female in 1st class) and Jack (a 20 year old male in 3rd class). Show your work.

```{r eval = TRUE}
## predict for Rose
## (by hand)
1.130523 + (-0.207434)*0 + (-0.393344)*0 + (-0.501326)*0 + (-0.006005)*17 

## (using predict)
predict(lin_mod2, newdata = data.frame(PClass = "1st", Sex = "female", Age = 17))

## predict for Jack
## (by hand)
1.130523 + (-0.207434)*0 + (-0.393344)*1 + (-0.501326)*1 + (-0.006005)*20

## (using predict)
predict(lin_mod2, newdata = data.frame(PClass = "3rd", Sex = "male", Age = 20))
```

b. Now use the logistic regression model to predict the survival probability for Rose and Jack. Show your work. (Hint: use the logistic regression model to obtain the predicted log odds, exponentiate to get the odds, and then convert to probability.)

```{r eval = TRUE}
## predict for Rose
## (by hand)
log_odds_rose <- 3.75966210 + (-1.29196240)*0 + (-2.52141915)*0 + (-2.63135683)*0 + (-0.03917681)*17
odds_rose <- exp(log_odds_rose)
odds_rose/(1+odds_rose)

## (using predict)
predict(log_mod2, newdata = data.frame(PClass = "1st", Sex = "female", Age = 17), type = "response")

## predict for Jack
## (by hand)
log_odds_jack <- 3.75966210 + (-1.29196240)*0 + (-2.52141915)*1 + (-2.63135683)*1 + (-0.03917681)*20
odds_jack <- exp(log_odds_jack)
odds_jack/(1+odds_jack)

## (using predict)
predict(log_mod2, newdata = data.frame(PClass = "3rd", Sex = "male", Age = 20), type = "response")
```

c. Comment on differences that you notice in the predictions from parts a and b.

> Our linear model predicted that Rose's probability of survival was over 100% (which doesn't make sense). The predictions for Jack are fairly similar: 10.2% based on our logistic model and 11.6% based on our linear model.


\
\
\

# Additional exercises {-}


## Exercise 7: Calculating probabilities and odds

In Exercise 6, we used the `predict()` function to calculate the odds and probabilities of success for climbers that do and don't use oxygen:

```{r}
# odds of success for a climber that doesn't use oxygen
exp(predict(climbing_model_2, newdata = data.frame(oxygen_used = FALSE)))

# probability of success for a climber that doesn't use oxygen
predict(climbing_model_2, newdata = data.frame(oxygen_used = FALSE), type = "response")

# odds of success for a climber that uses oxygen
exp(predict(climbing_model_2, newdata = data.frame(oxygen_used = TRUE)))

# probability of success for a climber that uses oxygen
predict(climbing_model_2, newdata = data.frame(oxygen_used = TRUE), type = "response")
```

Calculate these 4 values *by hand*, starting from the model formula:

log(odds of `success`) = -1.327993 + 2.903864 `oxygen_usedTRUE`




## Exercise 8: Interpreting logistic regression coefficients

In general (not just for climbers), let's summarize how to interpret logistic regression coefficients.
Review the details below and convince yourself that these general properties are consistent with your work above.
Let $Y$ be a binary categorical response variable (eg: Y = 1 if the event of interest happens and is 0 otherwise).
Then a logistic regression model of $Y$ by $Y$ is

$$\log(\text{odds}) = \beta_0 + \beta_1 X$$

- Coefficient interpretation for quantitative X   
    $$\begin{split}    
    \beta_0     & = \text{ LOG(ODDS) when } X=0 \\
    e^{\beta_0} & = \text{ ODDS when } X=0 \\
    \beta_1     & = \text{ change in LOG(ODDS) per 1 unit increase in } X \\
    e^{\beta_1} & = \text{ multiplicative change in ODDS per 1 unit increase in } X \text{ (ie. the odds ratio, comparing the odds for X + 1 vs X)} \\
    \end{split}$$
- Coefficient interpretation for categorical X    
    $$\begin{split}    
    \beta_0     & = \text{ LOG(ODDS) at the reference category } \\
    e^{\beta_0} & = \text{ ODDS at the reference category } \\
    \beta_1     & = \text{ unit change in LOG(ODDS) relative to the reference} \\
    e^{\beta_1} & = \text{ multiplicative change in ODDS relative to the reference } \text{ (ie. the odds ratio, comparing the odds for the given category vs the reference category)} \\
    \end{split}$$






## Exercise 9: R code

A lot of R code was used to produce the plots and other output above.
The code is not the point.
It's more important to understand the logistic regression concepts, and the *output* of the code -- you can always look up code.
With that in mind, review the code to get a sense for the structure!



## Exercise 10: Optional math challenge

Consider a simple logistic regression model:

$$\log(\text{odds}) = \beta_0 + \beta_1 X$$

Prove that $e^{\beta_1}$ indicates the MULTIPLICATIVE change in ODDS for a 1-unit increase in X, i.e. that it's equivalent to the odds ratio for "X+1" vs "X".     
 


\
\
\
\


# Probability and Odds Solutions



```{r eval = TRUE, echo = FALSE}
library(tidyverse)
firststeps <- read_csv("https://mac-stat.github.io/data/firststeps.csv") %>%
  mutate(firstchild = ifelse(parity == 0, "Yes", "No"), # Is this the first child this parent has had?
         low_bwt = ifelse(bwt < 2500, "low", "not low"),
         preterm = ifelse(gestation < 37, "Yes", "No")) # short gestational period
```


## Exercise 1: Exploring First Steps enrollment and Gestational Age


```{r eval = TRUE}
# 2x2 Table: preterm vs. First Steps
firststeps %>% 
    count(preterm, firstep)
```


a. 343 + 60 = 403 parents were enrolled in the First Steps program. I used both rows of the table where `firstep = 1`.

b. 16.12% of people in the study were enrolled in First Steps!

```{r eval = TRUE}
403 / 2500
```

c. 60 birth parents

d. 14.89% of birth parents in First Steps had a premature baby.

```{r eval = TRUE}
60 / 403
```


e. The total number of birth parents who had a premature baby was 218 + 60 = 278. Of those. 60 were enrolled in First Steps. Therefore, 21.58% of birth parents who had a premature baby were enrolled in First Steps.

```{r eval = TRUE}
60/278
```


Using formal probability notation, we can write

> b. $P(\text{First Steps})$ = .1612

> d. $P(\text{Preterm} | \text{First Steps})$ = .1488

> e. $P(\text{First Steps} | \text{Preterm})$ = .2158


f. 

> $P(\text{Preterm} | \text{Not in First Steps})$ = 218 / (218 + 1879) = 0.103958

g. 

$$
\frac{(\text{Preterm} | \text{First Steps})}{P(\text{Preterm} | \text{Not in First Steps})} = .1488 / 0.103958 = 1.43
$$

h. Parents in this study in the First Steps program are 1.43 times more likely to have a premature birth than those not enrolled in the First Steps program, indicating that gestational age does differ by First Steps enrollment. This implies that enrollment in the First Steps program may not be associated with better birth outcomes, as measured by gestational age.

*Note:* However, you may argue that this is not a fair comparison, or that this summary is not what researchers were actually interested in! Ideally, we would compare birth outcomes from mothers in the First Steps program to the birth outcomes from *those same mothers* not in the First Steps program, to determine if the program made a positive impact. This idea hints at a sub-field of statistics called **causal inference** and the idea of a **counterfactual** ("what would have happened if..."). Take more statistics classes to learn about other methods for approaching this question!


i. 

```{r eval = TRUE}
# Side-by-side bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar(position = "dodge") +
  theme_classic()

# Stacked bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar() +
  theme_classic()

# Stacked relative frequency bar chart
firststeps %>%
  ggplot(aes(firstep, fill = preterm)) +
  geom_bar(position = "fill") +
  theme_classic()
```

**Bonus Question:** Which of the above three plots allows you to directly see the conditional probabilities we calculated previously?

> The stacked relative frequency bar chart! 



## Exercise 2: Exploring First Steps enrollment and Low birthweights

a. 

```{r eval = TRUE}
# 2x2 Table: low_bwt vs. First Steps
firststeps %>%
  count(low_bwt, firstep)
```

b. 

> $P(\text{Low birth weight} | \text{First Steps})$ = 25 / (25 + 378) = 0.062

> $P(\text{Normal birth weight} | \text{First Steps})$ = 378 / (25 + 378) = 0.938



> $P(\text{Low birth weight} | \text{Not in First Steps})$ = 102 / (102 + 1995) = 0.049

> $P(\text{Normal birth weight} | \text{Not in First Steps})$ = 1995 / (102 + 1995) = 0.951


c. 

> $Odds(\text{Low birth weight} | \text{First Steps})$ = 0.062 / (1 - 0.062) = 0.06609808

> $Odds(\text{Normal birth weight} | \text{First Steps})$ = 0.938 / (1 - 0.938) = 15.12903



> $Odds(\text{Low birth weight} | \text{Not in First Steps})$ = 0.049 / (1 - 0.049) = 0.05152471

> $Odds(\text{Normal birth weight} | \text{Not in First Steps})$ = 0.951 / (1 - 0.951) = 19.40816


d. 

```{r eval = TRUE}
0.06609808 / 0.05152471
```


e. The odds of having a low birth weight baby are 1.28 times higher for those enrollment in First Steps compared to those not in First Steps. Just as in Exercise 1, this implies that the First Steps program may not be associated with improved birth outcomes (with the same caveats as given in the answer to 1 (h)).

f. To go along with your summary, add code below to make one of the three visualization options we tried out in Exercise 1.

```{r eval = TRUE}
# Stacked relative frequency bar chart (with some fancy aesthetics)
firststeps %>%
  mutate(Birthweight = low_bwt %>% str_to_title()) %>%
  ggplot(aes(firstep, fill = Birthweight)) +
  geom_bar(position = "fill") +
  theme_classic() +
  scale_fill_viridis_d(option = "H") +
  labs(x = "First Steps", title = "Birthweight by First Steps Enrollment") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Enrolled", "Enrolled")) 
```

## Exercise 3: Conditional vs. Marginal probabilities

Suppose we select a person at random from the entire global population. For each of the following probabilities, which do you think is bigger? Explain your reasoning.

a. P(lung cancer | smoker) is likely bigger, since lung cancer is more rare in the general population than it is among smokers.

b. P(likes McDonald's) is likely bigger, since vegetarians don't likely like McDonald's very much (few options that they can eat).

c. P(smart | Mac grad) is likely bigger, because there are *very* few Mac grads relative to the global population. Lots of people are smart, few are Mac grads.

## Exercise 4: Probability practice

Let's explore whether birth weight of a baby varies by whether or not it was the *first* child that a mother had, *and* whether this relationship differs by First Steps enrollment. We make a table below:

```{r eval = TRUE}
firststeps %>%
  count(firstchild, low_bwt, firstep)
```

a. What is the probability that a mother enrolled in First steps who is having their first child, has a baby who is born at a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(\text{Low Birthweight} | \text{First Steps}, \text{First Child}) = 12 / (12 + 180) = 0.0625

b. What is the probability that a mother *not* enrolled in First steps who is having their first child, has a baby who is born at a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(\text{Low Birthweight} | \text{Not in First Steps}, \text{First Child}) = 59 / (59 + 915) = 0.06057495

c. What is the probability that a mother's first child has a low birthweight? Calculate your answer, and write it using formal probability notation.

> P(\text{Low Birthweight} | \text{First Child}) = (59 + 12) / (59 + 12 + 915 + 180) = 0.06089194

d. How many times *more likely* is a child to be born at a low birthweight, comparing children who are the first born to those not first born?

> P(\text{Low Birthweight} | \text{First Child}) / P(\text{Low Birthweight} | \text{Not First Child}) = ((59 + 12) / (59 + 12 + 915 + 180)) / ((43 + 13) / (43 + 13 + 1080 + 198)) = 1.450533


\
\
\
\



# Logistic Univariate Solutions


```{r eval = TRUE, echo = FALSE}
library(tidyverse)
library(ggmosaic)
titanic <- read_csv("https://mac-stat.github.io/data/titanic.csv")
```

## Exercise 1: Exploring age

The boxplot doesn't clearly indicate a difference in the age distributions across survivors and non-survivors, but we do notice from the density plot that there is a greater density of younger passengers among the survivors. We also see from the last plot that younger passengers tend to have a higher survival chance.

```{r eval = TRUE}
# Create a boxplot
ggplot(titanic, aes(x = factor(Survived), y = Age)) +
    geom_boxplot()

# Can flip the boxplot on its side too
ggplot(titanic, aes(y = factor(Survived), x = Age)) +
    geom_boxplot()

# Create a density plot
ggplot(titanic, aes(x = Age, color = factor(Survived))) +
    geom_density()

# Use the code below to create a plot of the fraction who survived at each age
titanic_summ <- titanic %>% 
    group_by(Age) %>%
    summarize(frac_survived = mean(Survived))

ggplot(titanic_summ, aes(x = Age, y = frac_survived)) +
    geom_point() +
    geom_smooth(se = FALSE)
```



## Exercise 2: Exploring sex and ticket class

- Females were more likely to survive than males.
- 1st class was most likely to survive, followed by 2nd then 3rd class.

```{r eval = TRUE}
# Standard bar plots
ggplot(titanic, aes(x = Sex, fill = factor(Survived))) +
    geom_bar(position = "fill")

ggplot(titanic, aes(x = PClass, fill = factor(Survived))) +
    geom_bar(position = "fill")

# Mosaic plots
ggplot(data = titanic %>% mutate(Survived = as.factor(Survived))) +
    geom_mosaic(aes(x = product(Sex), fill = Survived))

ggplot(data = titanic %>% mutate(Survived = as.factor(Survived))) +
    geom_mosaic(aes(x = product(PClass), fill = Survived))
```



## Exercise 3: Linear regression model

```{r eval = TRUE}
titanic %>% 
    count(PClass, Survived)
```

a. 

Class     | Died     | Survived | Total 
----------|----------|----------|-------
1st Class | 129      |  193     | 322
2nd Class | 160      |  119     | 279
3rd Class | 573      |  138     | 711
Total     | 862      |  450     | 1312

b. 
    - the probability of surviving among 1st class passengers: 193/322 = `r round(193/322,3)`
    - the probability of surviving among 2nd class passengers: 119/279 = `r round(119/279,3)`
    - the probability of surviving among 3rd class passengers: 138/711 = `r round(138/711,3)`
    - the difference in the probability of surviving, comparing 2nd class passengers to 1st class passengers (i.e., how much lower is the probability of 2nd class passengers as compared to 1st class passengers?): 119/279 - 193/322 = `r round(119/279 - 193/322, 3)`
    - the difference in the probability of surviving, comparing 3rd class passengers to 1st class passengers  (i.e., how much lower is the probability of 3rd class passengers as compared to 1st class passengers?): 138/711 - 193/322 = `r round(138/711 - 193/322, 3)`

c. This model can be written as: $E[Survived | PClass] = \beta_0 + \beta_1 PClass2nd + \beta_2 PClass3rd$.
    - In the context of a binary variable, the expected value/average is the same as the probability that the variable equals one. To see an example of this, calculate the average of this list of 0's and 1's: (0,0,1,1,0,1,0,1). Now calculate the proportion of 1's. What do you notice?
    - This means that we can also write this model as follows: $P[Survived = 1 | PClass] = \beta_0 + \beta_1 PClass2nd + \beta_2 PClass3rd$

```{r eval = TRUE}
lin_mod <- lm(Survived ~ PClass, data = titanic)
summary(lin_mod)
```

d. The coefficient estimates are the differences in probability from part b!
    - `(Intercept)`: the estimated probability of survival for passengers in 1st class is 0.599 (59.9%)
    - `PClass2nd`: the difference in the estimated probability of survival comparing passengers in 1st class to passengers in 2nd class is 0.173 (17.3%), where passengers in 1st class have the higher estimated survival probability
        - OR... comparing passengers in 1st class to passengers in 2nd class, the difference in the proportion of passengers that survived is 0.173 (17.3%), with 1st class having a higher proportion of passengers that survived
        - OR... the probability of survival is 17.3% lower among passengers in 2nd class than it is among passengers in 1st class
    - `PClass3rd`: the difference in the estimated probability of survival comparing passengers in 1st class to passengers in 3rd class is 0.405 (40.5%), where passengers in 1st class have the higher estimated survival probability



## Exercise 4: Logistic regression model (categorical predictor)

a. 
    - the odds of surviving among 1st class passengers: 193/129 = `r round(193/129,3)`
    - the odds of surviving among 2nd class passengers: 119/160 = `r round(119/160,3)`
    - the odds of surviving among 3rd class passengers: 138/573 = `r round(138/573,3)`
    - the ratio of the odds of surviving, comparing 2nd class passengers to 1st class passengers (i.e., how many times higher/lower is the odds of survival among 2nd class passengers as compared to 1st class passengers?): (119/160)/(193/129) = `r round((119/160) / (193/129), 3)`
    - the ratio of the odds of surviving, comparing 3rd class passengers to 1st class passengers: (138/573)/(193/129) = `r round((138/573) / (193/129), 3)`

b. $\log(Odds[Survived = 1 | PClass]) = \beta_0 + \beta_1 PClass2nd + \beta_2 PClass3rd$

```{r eval = TRUE}
log_mod <- glm(Survived ~ PClass, data = titanic, family = "binomial")

# These logistic coefficient estimates are NOT exponentiated
coef(summary(log_mod))
```

```{r eval = TRUE}
# Calculations for exponentiating coefficients
exp(0.4028778)
exp(-0.6989281)
exp(-1.8265098)
```


d. These exponentiated coefficient estimates compare to your the odds and odds ratios in part a!
    - exp(Intercept): the estimated odds of survival among passengers in first class is 1.496 (i.e., passengers in first class are 1.496 times more likely to survive than they are to die)
    - exp(PClass2nd): we estimate that the odds of survival for passengers in 2nd class are only 0.50 times as high as the odds of survival among passengers in 1st class (i.e., the odds of survival are 2 times higher among passengers in 1st class than they are among passengers in 2nd class)
    - exp(PClass3rd): we estimate that the odds of survival for passengers in 3rd class are only 0.16 times as high as the odds of survival among passengers in 1st class (i.e., the odds of survival are 1/0.16 = `r round(1/0.1609744, 2)` times higher among passengers in 1st class than they are among passengers in 3rd class)

## Exercise 5: Logistic regression model (quantitative predictor)

a. After fitting the logistic regression model below, write out the model formula using correct notation.

```{r eval = TRUE}
log_mod <- glm(Survived ~ Age, data = titanic, family = "binomial")
coef(summary(log_mod))
```

$\log(Odds[Survived = 1 | Age]) = \beta_0 + \beta_1 Age$

b. Write an interpretation of each of the *exponentiated* coefficients in this logistic regression model.

- exp(Intercept): $exp(-0.0814)=0.92$ --> the estimated odds of survival among passengers who are 0 years old is 0.92 (i.e., passengers who are 0 years old are 0.92 times more likely to survive than they are to die--so very slightly *more* likely to die)
- exp(Age): $exp(-0.0088)=0.99$ --> For every 1-year increase in a passenger's age, the estimated odds of survival *decrease* by about 1%.

## Exercise 6: Linear vs. logistic modeling

To highlight a key difference between linear vs. logistic modeling, consider the following linear and logistic regression models of survival with sex and age as predictors in addition to ticket class.

```{r eval = TRUE}
lin_mod2 <- lm(Survived ~ PClass + Sex + Age, data = titanic)
coef(summary(lin_mod2))

log_mod2 <- glm(Survived ~ PClass + Sex + Age, data = titanic, family = "binomial")
coef(summary(log_mod2))
```

a.

```{r eval = TRUE}
## predict for Rose
## (by hand)
1.130523 + (-0.207434)*0 + (-0.393344)*0 + (-0.501326)*0 + (-0.006005)*17

## (using predict)
predict(lin_mod2, newdata = data.frame(PClass = "1st", Sex = "female", Age = 17))

## predict for Jack
## (by hand)
1.130523 + (-0.207434)*0 + (-0.393344)*1 + (-0.501326)*1 + (-0.006005)*20

## (using predict)
predict(lin_mod2, newdata = data.frame(PClass = "3rd", Sex = "male", Age = 20))
```


b. 

```{r eval = TRUE}
## predict for Rose
## (by hand)
log_odds_rose <- 3.75966210 + (-1.29196240)*0 + (-2.52141915)*0 + (-2.63135683)*0 + (-0.03917681)*17
odds_rose <- exp(log_odds_rose)
odds_rose/(1+odds_rose)

## (using predict)
predict(log_mod2, newdata = data.frame(PClass = "1st", Sex = "female", Age = 17), type = "response")

## predict for Jack
## (by hand)
log_odds_jack <- 3.75966210 + (-1.29196240)*0 + (-2.52141915)*1 + (-2.63135683)*1 + (-0.03917681)*20
odds_jack <- exp(log_odds_jack)
odds_jack/(1+odds_jack)

## (using predict)
predict(log_mod2, newdata = data.frame(PClass = "3rd", Sex = "male", Age = 20), type = "response")
```

c. Our linear model predicted that Rose's probability of survival was over 100% (which doesn't make sense). The predictions for Jack are fairly similar: 10.2% based on our logistic model and 11.6% based on our linear model.


## Exercise 7: Calculating probabilities and odds

For people that don't use oxygen:

```{r}
# log(odds)
-1.327993

# odds
exp(-1.327993)

# probability
exp(-1.327993) / (exp(-1.327993) + 1)
```


For people that do use oxygen:

```{r}
# log(odds)
-1.327993 + 2.903864

# odds
exp(-1.327993 + 2.903864)

# probability
exp(-1.327993 + 2.903864) / (exp(-1.327993 + 2.903864) + 1)
```
