---
title: "🗓️ Week 03 -Fundamentals of Classification ([DS202](https://lse-dsi.github.io/DS202))"
subtitle: "Theme: Supervised Learning"
author: "Dr. Ghita Berrada"
date: "18 October 2024"
format: 
  html:
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 03!

Instead of slides, we are going to use a Quarto markdown as the lecture material for this week. This is to encourage you to follow along on your computers and help you practice taking notes of concepts, equations, and code using markdown.

⚙️ **Set up**

```{r import, message=FALSE, warning=FALSE}
# Tidyverse packages, spelled out individually for clarity
library(dplyr)
library(ggplot2)
library(archive)
library(readr)
library(stringr)
library(forcats)
library(tidyr)

# Tidymodels packages, spelled out individually for clarity
library(parsnip)
library(yardstick)
library(recipes)
library(workflows)
library(rsample) # to split data in training and test sets

```

# 🥅 Learning Objectives

What is the point of today?

The hope is that, by the end of today, you will learn the following things:

**1st Hour**

-   Review regression and understand the term 'goodness-of-fit'
-   Obtain a deeper insight into training and testing splits
-   Apply different metrics to assess the of regression algorithms

**2nd Hour**

-   Distinguish regression versus classification

-   Run a logistic regression (a classification method!) model

# 1: Regression once again

Let's go back to last week for a bit.

Remember that Stack Overflow Developer's Survey?

Let's download the 2022 version of it again.

```{r}
df<-read_csv(archive_read("https://info.stackoverflowsolutions.com/rs/719-EMH-566/images/stack-overflow-developer-survey-2022.zip", file = 3), col_types = cols())
```

Let's pre-process it to get it ready for modeling once again (i.e only selecting the data for the UK, removing the categories of developers irrelevant to our analysis, etc.)

```{r}
filtered_gender <- c("Man", "Woman", "Non-binary")
survey_results <- df %>% 
  dplyr::filter(
    Country == "United Kingdom of Great Britain and Northern Ireland",
    Employment == "Employed, full-time",
    ConvertedCompYearly > 0,
    ConvertedCompYearly < 2e6
  )

## identify non-ICs, to remove
managers_ctos <- survey_results %>%
  dplyr::filter(str_detect(DevType, "Engineering manager|Product manager|Senior executive/VP"))

## identify academics, to remove
academics <- survey_results %>%
  dplyr::filter(str_detect(DevType, "Academic researcher|Scientist|Educator"))

filtered_devtype <- c(
  "Other", "Student",
  "Marketing or sales professional"
)

results <- survey_results %>%
  anti_join(managers_ctos) %>%
  anti_join(academics) %>%
  transmute(ResponseId,
            EdLevel = fct_collapse(EdLevel,
                                   `Less than bachelor's` = c(
                                     "Primary/elementary school",
                                     "Secondary school (e.g. American high school, German Realschule or Gymnasium, etc.)",
                                     "Some college/university study without earning a degree",
                                     "Associate degree (A.A., A.S., etc.)"
                                   ),
                                   `Other` = "Something else",
                                   `Bachelor's degree` = "Bachelor’s degree (B.A., B.S., B.Eng., etc.)",
                                   `Graduate degree` = c(
                                     "Other doctoral degree (Ph.D., Ed.D., etc.)",
                                     "Master’s degree (M.A., M.S., M.Eng., MBA, etc.)",
                                     "Professional degree (JD, MD, etc.)"
                                    )
            ),
            DevType,
            RemoteWork,
            YearsCodePro = parse_number(YearsCodePro),
            Gender = case_when(
              str_detect(Gender, "Non-binary") ~ "Non-binary",
              TRUE ~ Gender
            ),
            ConvertedCompYearly
  ) %>%
  dplyr::filter(Gender %in% filtered_gender)

results_parsed <- results %>%
  mutate(DevType = str_split(DevType, pattern = ";")) %>%
  unnest(DevType) %>%
  mutate(
    DevType = case_when(
      str_detect(str_to_lower(DevType), "data scientist") ~ "Data scientist",
      str_detect(str_to_lower(DevType), "data or business") ~ "Data analyst",
      str_detect(str_to_lower(DevType), "desktop") ~ "Desktop",
      str_detect(str_to_lower(DevType), "embedded") ~ "Embedded",
      str_detect(str_to_lower(DevType), "devops") ~ "DevOps",
      str_detect(DevType, "Engineer, data") ~ "Data engineer",
      str_detect(str_to_lower(DevType), "site reliability") ~ "DevOps",
      TRUE ~ DevType
    ),
    DevType = str_remove_all(DevType, "Developer, "),
    DevType = str_to_sentence(DevType),
    DevType = str_replace_all(DevType, "Qa", "QA"),
    DevType = str_replace_all(DevType, "Sre", "SRE"),
    DevType = str_replace_all(DevType, "Devops", "DevOps")
  ) %>%
  dplyr::filter(
    !DevType %in% filtered_devtype,
    !is.na(DevType)
  )

cleaned_df <- results_parsed %>%
  drop_na(YearsCodePro)%>%
  dplyr::filter(
    ConvertedCompYearly < 1e7,
    YearsCodePro < 60
  ) %>%
  dplyr::filter(Gender %in% c("Man", "Woman")) %>%
  select(-ResponseId) %>%
  mutate(ConvertedCompYearly = log(ConvertedCompYearly))
 
 
```

Last week for the sake of illustration, we very conveniently built a multivariate linear model on the whole dataset:

```{r}
library(tidymodels)

# Create a linear model
lm_spec <- 
  linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")


lm_fit <- 
  lm_spec %>%
  fit(ConvertedCompYearly ~ 0 + DevType + ., data=cleaned_df)

print(lm_fit)

#tidy(lm_fit)

glance(lm_fit)
```

In practice though, we generally don't do this! We want our models to generalize to unseen data. What if my survey was not about this year's data? Would the model still be valid then? What if we included data from another region: would the observations generalize?

That's why we split the data into training and test data to increase the robustness of our modeling.

**Training vs testing split**...

To assess the capabilities of this model to generalise to unseen data, I will split the data into training and test set. I have to be careful here though: as mentioned before, men are overrepresented in the data set while women and binaries are underrepresented so I have to split the data in such a way as to keep the proportion of genders intact.

```{r}
# Split the data with 75% being used to train the model
survey_split <- initial_split(cleaned_df, prop = 0.75,strata=Gender) #the option strata=Gender allows me to keep the original proportions of genders from the dataset 
# Create tibbles of the training and test set

survey_train <- training(survey_split)
survey_test <- testing(survey_split)
```

Now let's build linear models predicting salary:

```{r}
model1 <- 
  linear_reg() %>%
  set_engine("lm") %>%
  fit(ConvertedCompYearly ~ YearsCodePro, data = survey_train)
```

```{r}
model2 <- 
  linear_reg() %>%
  set_engine("lm") %>%
  fit(ConvertedCompYearly ~ 0 + DevType + ., data = survey_train)
```

```{r}
model3 <- 
  linear_reg() %>%
  set_engine("lm") %>%
  fit(ConvertedCompYearly ~ Gender, data = survey_train)
```

How do we know if one linear model is better than another?

-   By comparing the **goodness-of-fit** of the models using standardised metrics
-   By comparing predictions of the models on a test set, **using the same unit of measurement**

## 1.1 Comparing models

We could use RMSE here if we wanted (RMSE has same unit of measurement as the variable you are trying to predict!) since we're always predicting the same variable so our units would be consistent throughout . But there's one metric that appears when you run `summary()` on a linear model that is agnostic to the unit of measurement: the $R^2$ metric. We'll use $R^2$ here.

```{r}
summary(model1$fit)
```

```{r}
summary(model2$fit)
```

```{r}
summary(model3$fit)
```

## 1.2 $R^2$ as a measures of goodness-of-fit

When you were asked to calculate, for example, the RMSE of a model when applied to the training set, you were calculating the goodness-of-fit.

> 📕 Goodness-of-fit is the concept of how well a model fits the data it was trained on.

The RMSE has the same unit of measurement as the variable you are trying to predict, but there are some measures that are agnostic to the unit of measurement and it might be useful to know them.

The R-squared metric uses the concept of residuals, which you have explored in the previous weeks. The $R^2$ metric is defined as:

$$
R^2 = 1 - \frac{RSS}{TSS}
$$

where $RSS = \sum_{i=1}^n (y_i - \hat{y}_i)^2$ (the sum of squares of the residuals) and $TSS = \sum_{i=1}^n (y_i - \bar{y})^2$ is the total sum of squares of all the data points used to fit the model [^1]. The closer $R^2$ is to 1, the better the model *explains the variance* in the data. The closer it is to 0, the worse the model is.

[^1]: See [An Introduction to Statistical Learning](https://www.statlearning.com/) by Hastie et al. Page 70 of the 2nd edition (freely available online).

The $R^2$ metric is closely related to the Pearson correlation coefficient ($\rho$):

$$
\rho = \sqrt{R^2}
$$

where $\rho$ is the Pearson correlation coefficient.

The interpretation is a little different: $\rho$ ranges from -1 to 1. Values of $|\rho|$ closer to 1 indicate that the two variables are more strongly **linearly correlated**. Values of $|\rho|$ closer to 0 indicate that the two variables are less strongly **linearly correlated**.

```{r}
# Function cor returns the Pearson correlation coefficient
cor(survey_train$ConvertedCompYearly, survey_train$YearsCodePro)
```

Note though that the Pearson coefficient only applies to numerical variables.

## ⚠️ Note on Metrics and Non-Linearity:

These two metrics are only useful in the context of linear models. It is definitely possible for two variables to be strongly associated with each other but in a *non-linear* way. Think of studying and exam scores: more studying leads to better scores, but only to a point (diminishing returns).

For now, we use $$R^2$$ to compare `model1`, `model2` and `model3` in a **linear context**. While not absolute proof ("all models are wrong, but some are useful"), **model2** seems superior based on these measures.

But *'How do I know **for sure**'*, I hear you ask? To more robustly assess performance, let's test all three models on unseen data using a **test set**. That's our next step!

# 2. Don't stick just to goodness-of-fit, though, go deeper

It is very easy to conclude by looking at measures of goodness-of-fit that `model1` is superior to `model2`. It is true that the first model is a better linear fit. However, what I want you to remember when comparing different models is the following:

> *All models are wrong, but just how useful are these models?*

## 2.1 Why use a test set?

The whole point of having a test set is to check if our model is as good as it sounds on paper. We want our models to be useful in the real world, not just when used on the training data.

Thinking this way, let's compare the predictions made by `model1` and `model2` and although I could still continue to compare them with unit-agnostic measures like $R^2$, I want to bring you back to the RMSE - it will be easier to interpret and make it more 'real-world'.

Predictions made by `model1` are already made in dollars:

```{r}
model1 %>% predict(survey_test) %>% head()
```

The predictions made by `model2` are also made in dollars:

```{r}
model2 %>% predict(survey_test) %>% head()
```

The same can be said for `model3` :

```{r}
model3 %>% predict(survey_test) %>% head()
```

## 2.2 Residuals plot (on test set)

### 2.2.1 Easy to do for `model1`

```{r message=FALSE}
plot_df <- model1 %>% augment(survey_test)

g <- (
  ggplot(plot_df, aes(x = ConvertedCompYearly, y=.resid)) +

  ## Geoms
  geom_hline(yintercept=0, linetype="dashed", color="#212121") +
  geom_point(size=2, alpha=0.3, stroke=1) +

  ## Customisation of scales and labels
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     limits=c(-1e5, 1e5)) +
  labs(x = "Salary (data)", 
       y = "Residual yearly salary",
       title="Residuals plot for model1") + 
    
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw() +
  theme(plot.title=element_text(size=rel(2.5)),
        plot.subtitle=element_text(size=rel(1.5)),
        axis.text=element_text(size=rel(1.5)),
        axis.title=element_text(size=rel(2)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        legend.title=element_text(size=rel(2)),
        legend.text=element_text(size=rel(1.5)),
        legend.key.size=unit(25, "pt")
       )
) 

g
```

### 2.2.2 Same work for `model2` and for `model3`

Let's create the residuals plot for `model2` now.

```{r message=FALSE}
plot_df <- model2 %>% augment(survey_test)

g <- (
  ggplot(plot_df, aes(x = ConvertedCompYearly, y=.resid)) +

  ## Geoms
  geom_hline(yintercept=0, linetype="dashed", color="#212121") +
  geom_point(size=2, alpha=0.3, stroke=1) +

  ## Customisation of scales and labels
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     limits=c(-1e5, 1e5)) +
  labs(x = "Salary (data)", 
       y = "Residual yearly salary",
       title="Residuals plot for model2") + 
    
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw() +
  theme(plot.title=element_text(size=rel(2.5)),
        plot.subtitle=element_text(size=rel(1.5)),
        axis.text=element_text(size=rel(1.5)),
        axis.title=element_text(size=rel(2)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        legend.title=element_text(size=rel(2)),
        legend.text=element_text(size=rel(1.5)),
        legend.key.size=unit(25, "pt")
       )
) 

g
```

Then, let's create the residuals plot but `model3` now.

```{r message=FALSE}
plot_df <- model3 %>% augment(survey_test)

g <- (
  ggplot(plot_df, aes(x = ConvertedCompYearly, y=.resid)) +

  ## Geoms
  geom_hline(yintercept=0, linetype="dashed", color="#212121") +
  geom_point(size=2, alpha=0.3, stroke=1) +

  ## Customisation of scales and labels
  scale_x_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3)) +
  scale_y_continuous(labels=scales::unit_format(unit = "K", scale = 1e-3),
                     limits=c(-1e5, 1e5)) +
  labs(x = "Salary (data)", 
       y = "Residual yearly salary",
       title="Residuals plot for model3") + 
    
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw() +
  theme(plot.title=element_text(size=rel(2.5)),
        plot.subtitle=element_text(size=rel(1.5)),
        axis.text=element_text(size=rel(1.5)),
        axis.title=element_text(size=rel(2)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        legend.title=element_text(size=rel(2)),
        legend.text=element_text(size=rel(1.5)),
        legend.key.size=unit(25, "pt")
       )
) 

g
  
```

## 2.3 Other metrics (test set)

### 2.3.2 RMSE

You know this one already!

**Model 1**

```{r}
model1 %>%
  augment(survey_test) %>%
  rmse(ConvertedCompYearly, .pred) 
```

**Model 2**

```{r}
model2 %>%
  augment(survey_test) %>%
  rmse(ConvertedCompYearly, .pred) 
```

```{r}
model3 %>%
  augment(survey_test) %>%
  rmse(ConvertedCompYearly, .pred) 
```

When it comes to this particular metric, `model2` is still the winner here!

### 2.3.2 MAPE

RMSE is already normalized and provides a way to compare the usual fluctuations of your predicted variable (e.g in this case the usual fluctations i.e variance of the annual salary of a developer). Another one is MAPE or Mean Absolute Percentage Error.

There is a handy [`yardstick::mape()` function](https://yardstick.tidymodels.org/reference/mape.html) to calculate this for us, but the equation is as follows:

$$
\text{MAPE} = \frac{1}{n}\sum_{i=1}^n{\left|\frac{y_i - \hat{y_i}}{y_i}\right|}
$$

That is: on average, how much does the residual differs from the actual baseline value?

```{r}
model1_test_results <-
    model1 %>%
    augment(survey_test) %>%
    mape(ConvertedCompYearly, .pred)
model2_test_results <-
    model2 %>%
    augment(survey_test) %>%
    mape(ConvertedCompYearly, .pred)
model3_test_results <-
    model2 %>%
    augment(survey_test) %>%
    mape(ConvertedCompYearly, .pred)
model1_test_results
model2_test_results
model3_test_results
```

Here, `model2` and `model3` are tied. So overall, `model2` is our clear winner.

⚠️ **WARNING:** R-squared [will always improve](https://lse-dsi.github.io/lse-ds202-course-notes/blog/posts/r-squared.html) if you add more predictors, even if the predictors are random. So, if you are comparing models with different predictors, use the adjusted R-squared instead.

## 

# 4. Classification: a different type of supervised learning model

## 4.1 Making binary predictions

Now say I have a brand new ambition: to predict whether developers have more than 10 years of professional experience i.e they're not rookies (`Yes` or `No`).

Allow me to start over with my pre-processing and to create a target variable `not_rookie`:

```{r message=FALSE, error=FALSE}

dataset <- 
    cleaned_df %>% 
    mutate(
      not_rookie = if_else(YearsCodePro > 10, "Yes", "No")
    ) %>%
    tidyr::drop_na(not_rookie)


# run View(dataset) to view the new dataset
```

As before, let's split our data in training/test set:

```{r}
# Split the data with 75% being used to train the model
survey_split <- initial_split(dataset, prop = 0.75,strata=Gender) #the option strata=Gender allows me to keep the original proportions of genders from the dataset 
# Create tibbles of the training and test set

survey_train <- training(survey_split)
survey_test <- testing(survey_split)
```

```{r}
 survey_train %>% group_by(not_rookie) %>% tally()
```

Most are not rookies!

OK, but how do I write a model for the `not_rookies` variable?

## 4.2 Could I do linear regression?

In theory, yes, but only if you convert the `Yes` and `No` **labels** to numbers. Say: 1 and 0, respectively.

```{r}
model_yes_no <-
  linear_reg() %>%
  set_engine("lm") %>%
  fit(not_rookie ~ ConvertedCompYearly, 
      data=survey_train %>% mutate(not_rookie=if_else(not_rookie == "Yes", 1, 0)))
model_yes_no
```

What does this mean? Well, since linear regression is a **regression** algorithm, i.e., it operates on numbers, I can do as before and plot the predictions versus the actual values:

```{r}
model_yes_no %>%
augment(new_data =survey_test%>% mutate(not_rookie=if_else(not_rookie == "Yes", 1, 0))) %>%
ggplot(aes(x = ConvertedCompYearly, y=.resid)) +
geom_point()
```

NO. Linear Regression is not suitable for predicting labels. It would be more useful if we were to predict **probabilities** instead of single numbers.

## 4.3. Logistic Regression

Consider a **binary** response:

$$
Y = \begin{cases}
0 \\
1
\end{cases}
$$

We model the probability that $Y = 1$ using the **logistic function** (aka. sigmoid curve):

$$
Pr(Y = 1|X) = p(X) = \frac{e^{\beta_0 + \beta_1X}}{1 + e^{\beta_0 + \beta_1 X}}
$$

This function is mathematically guaranteed to be bounded by the interval $[0, 1]$.

### 4.3.1 Insight into the logistic function

Let's see how this mathematical function looks like:

```{r}
plot_logistic_function <- function(X_range, beta0, beta1, color){
  
  df <- expand.grid(X=X_range, beta0=beta0, beta1=beta1)
  
  df <- 
    df %>% group_by(beta0, beta1) %>% 
    transform(X=X, Y=exp(beta0+beta1*X)/(1+exp(beta0+beta1*X)))
  
  df <- df %>% 
    transform(beta0=factor(beta0, ordered = TRUE),
              beta1=factor(beta1, ordered = TRUE))
  
  g <- 
    if(missing(color)){
      ggplot(df, aes(x=X, y=Y))
    }else if(color=="beta0"){
      ggplot(df, aes(x=X, y=Y, color=beta0, linetype=beta0))
    }else{
      ggplot(df, aes(x=X, y=Y, color=beta1, linetype=beta1))
    }
  
  g <- 
    g + 
    geom_line(linewidth=1.4) +
  
    xlab("X") + ylab("Y") +
  
    # THEMING
    theme_bw() +
  	theme(
          axis.title.x = element_text(size=24,margin = margin(t = 10, r = 0, b = 0, l = 0)),
          axis.title.y = element_text(size=24, margin = margin(t = 0,  r = 10, b = 0, l = 0)),
          legend.text = element_text(size=18),
          legend.position="right",
          plot.title=element_text(size=26),
          axis.text.x = element_text(size=18),
          axis.text.y = element_text(size=18),
          panel.grid.minor = element_blank(),
          aspect.ratio = 1/2)
  return(g)
}
```

```{r}
plot_logistic_function(X_range=(-400:400)/100, beta0=c(-2, 0, 1), beta1=c(1), color="beta0") + 
labs(title=expression(paste("If I keep ", ~beta[1]==1," fixed and vary ", ~beta[0] )))
```

```{r}
plot_logistic_function(X_range=seq(-4, 4, 0.1), beta0=c(0), beta1=c(-1, 0, 1, 2), color="beta1") + 
labs(title=expression(paste("If I keep ", ~beta[0]==0," fixed and vary ", ~beta[1] )))
```

### 4.3.2 Building a model

Looking up the [list of models available on `tidymodels`](https://www.tidymodels.org/find/parsnip/), I find that there are many implementations of logistic regression. Here I am going to use the one from the `glm` package.

That is, something like:

```{r}
#| error: true
logistic_model1 <-
  logistic_reg() %>%
  set_engine("glm") %>%
  fit(not_rookie ~ ConvertedCompYearly, data=survey_train)
```

However, if you try to run the code above, you will get an error saying:

> Error in `check_outcome()`: ! For a classification model, the outcome should be a `factor`, not a `character`.

So we need to convert our `not_rookie` from a character ("Yes" vs "No") variable to a factor, a categorical variable.

Let's go back to creating our dataset:

```{r message=FALSE, error=FALSE}

dataset <- 
    dataset <- 
    cleaned_df %>% 
    mutate(
      not_rookie = if_else(YearsCodePro > 10, "Yes", "No"),
      not_rookie = factor(not_rookie, levels=c("No", "Yes"))
    ) %>%
    tidyr::drop_na(not_rookie)
  

# Split the data with 75% being used to train the model
survey_split <- initial_split(dataset, prop = 0.75,strata=Gender) #the option strata=Gender allows me to keep the original proportions of genders from the dataset 
# Create tibbles of the training and test set

survey_train <- training(survey_split)
survey_test <- testing(survey_split)

# Run View(dataset) to browse the new data
```

```{r}
set.seed(123)
logistic_model1 <-
  logistic_reg() %>%
  set_engine("glm") %>% 
  fit(not_rookie ~ ConvertedCompYearly, data=survey_train)

# Look at the coefficients
logistic_model1$fit
```

It worked! And the output is quite similar to that of a linear regression.

Let's see what happens when I augment the dataset with predictions:

```{r}
logistic_model1 %>% augment(survey_train)%>% arrange() %>% head()
```

I have several columns that relate to predictions!

### 4.3.3 Plot the model

Let's plot the output of the model so we can better understand what it represents:

```{r}
plot_df <- logistic_model1 %>% augment(survey_train)

min_x <- plot_df %>% filter(ConvertedCompYearly == min(ConvertedCompYearly)) %>% slice(1)
max_x <- plot_df %>% filter(ConvertedCompYearly == max(ConvertedCompYearly)) %>% slice(1)

g <- (
  ggplot(plot_df, aes(x = ConvertedCompYearly, y=.pred_Yes)) +

  ## Geoms
  geom_point(size=2, alpha=0.3, stroke=1) +

  ## Customisation of scales and labels
  scale_x_continuous(labels=scales::label_dollar(prefix = "$", suffix="K", scale=1e-3)) +
  labs(x = "Developers salary (in$)", 
       y = "Probability of being a rookie",
       title="Logistic regression is a good model for this problem",
       subtitle="The predictions are bounded between 0 and 1",
       caption="However, I am not truly sure what to make of this yet!") + 
  ylim(0, 1) +
    
  ## Theming. 
  ## Increase/decrease the font sizes according to your screen
  theme_bw() +
  theme(plot.title=element_text(size=rel(2.5)),
        plot.subtitle=element_text(size=rel(1.5)),
        axis.text=element_text(size=rel(1.5)),
        axis.title=element_text(size=rel(2)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        legend.title=element_text(size=rel(2)),
        legend.text=element_text(size=rel(1.5)),
        legend.key.size=unit(25, "pt"),
        plot.caption=element_text(size=rel(1.2), margin=margin(t=15))
       )
)
g
```

### 4.3.4 How do we interpret the coefficients?

We know that:

$$
\begin{equation}
\beta_0 \approx -4.4224\\
\beta_1 \approx 0.3394
\end{equation}
$$

This means that:

1.  In the case where the annual salary was zero, the probability of being a rookie is:

$$
\frac{e^{-4.4224}}{1 + e^{-4.4224}} = 0.721
$$

2.  For every 1 USD increase in the annual salary, the probability of being a rookie increases by a factor of

    $$
    \frac{e^{-4.4224 + 0.3394}}{1 + e^{-4.4224 + 0.3394}}= 0.01658
    $$

🤔 "*increased by a factor of... a number that is much lower than 1".* Yes, I know, this can be very confusing.

## 4.4 Making predictions

At what point do you say: Given these probabilities, I can confidently say who is a rookie and who isn't. By default, the **threshold** is set at probability=0.5.

```{r}
logistic_model1 %>% 
  augment(new_data=survey_train) %>%
  summarise(min_prob=min(.pred_Yes),
            q1=quantile(.pred_Yes, 0.25),
            median=quantile(.pred_Yes, 0.5),
            q3=quantile(.pred_Yes, 0.75),
            max=max(.pred_Yes))
```

But the threshold could change...

We'll see next time how to evaluate this model more systematically and classify in different ways.

# What's Next

In the lab next week, you will be tasked with creating a logistic regression and more...