---
title: "Intermediate Regression in R"
author: "Joschka Schwarz"
---

```{r}
#| include: false
source(here::here("R/setup-ggplot2-tie.R"))
options(dplyr.summarise.inform = FALSE)
```

Linear regression and logistic regression are the two most widely used statistical models and act like master keys, unlocking the secrets hidden in datasets. This course builds on the skills you gained in "Introduction to Regression in R", covering linear and logistic regression with multiple explanatory variables. Through hands-on exercises, you’ll explore the relationships between variables in real-world datasets, Taiwan house prices and customer churn modeling, and more. By the end of this course, you’ll know how to include multiple explanatory variables in a model, understand how interactions between variables affect predictions, and understand how linear and logistic regression work.

# Parallel Slopes

Extend your linear regression skills to "parallel slopes" regression, with one numeric and one categorical explanatory variable. This is the first step towards conquering multiple linear regression.

## Parallel slopes linear regression

Theory. Coming soon ...


**1. Parallel slopes linear regression**

Hi, I'm Richie. Welcome!

**2. The previous course**

This course builds on the skills from the previous course.

**3. From simple regression to multiple regression**

In that last course, you performed linear and logistic regression with a single explanatory variable. This time, you'll learn to fit models that include multiple explanatory variables. This is sometimes called "multiple regression".Including more explanatory variables in the model often gives you more insight into the relationship between the explanatory variables and the response, and can provide more accurate predictions. It's an important step towards mastering regression.

**4. The course contents**

Here's the plan. In Chapter 1, you'll explore parallel slopes linear regression. This is a special case of multiple linear regression, with one numeric explanatory variable and one categorical explanatory variable.Chapter 2 introduces interactions between variables and covers Simpson's Paradox, a counter-intuitive result affecting models containing categorical explanatory variables.Chapter 3 extends linear regression to even more explanatory variables, and gives some deeper insight into how linear regression works.Finally, Chapter 4 introduces multiple logistic regression, the logistic distribution, and digs into how logistic regression works.

**5. The fish dataset**

Here's the fish dataset from the previous course. Each row represents a fish, the mass is the response variable, and there is one numeric and one categorical explanatory variable.

**6. One explanatory variable at a time**

Recall that you run a linear regression by calling lm, passing a formula and a data frame. The formula has the response variable on the left and the explanatory variable on the right, with the variables separated by a tilde.Here you can see mass modeled against length. Printing the model shows the model coefficients. With a single numeric explanatory variable, you get one intercept coefficient and one slope coefficient.Let's change the explanatory variable to species. Recall that when you have a categorical explanatory variable, the coefficients are a little easier to understand if you use "plus zero" to tell R not to include an intercept in the model.Now you get one intercept coefficient for each category. That is, one coefficient for each species of fish.

**7. Both variables at same time**

To include both explanatory variables in the model, you combine them in the right-hand side of the formula, separated with a plus, just like you did with the zero.This time there is one slope coefficient, and one intercept coefficient for each category in the categorical variable.

**8. Comparing coefficients**

Examining the coefficients of each model, it's clear that the numbers are different.Notice that the slope coefficient for length, labeled length_cm, changes from thirty five to forty three once you include species in the model as well. The intercept coefficients for each species show an even bigger change. For example, once you add length into the model, bream, labeled speciesBream, changes from six hundred and eighteen to minus six hundred and seventy two.

**9. Visualization: 1 numeric explanatory var**

Here's the standard visualization for a linear regression with a numeric explanatory variable. You draw a scatter plot, then use geom_smooth with method equals "lm" to add a linear trend line. Setting se to FALSE prevents a standard error ribbon being drawn.

**10. Visualization: 1 categorical explanatory var**

For a categorical response, there are a few possible plots. The simplest one is to draw a box plot for each category. The model coefficients are the means of each category, which I've added using stat_summary with the fun-dot-y argument set to mean. shape equals 15 makes the mean point square.

**11. Visualization: both explanatory vars**

With a numeric and a categorical explanatory variable, you can draw a scatter plot as before. ggplot2 doesn't have an easy way to plot the model results, but fortunately one is provided by the moderndive package.In the plot the lines are parallel to each other. Consequently, this type of regression is nicknamed "parallel slopes regression", and the function to draw the lines is geom_parallel_slopes.

**12. Let's practice!**

Let's begin.

## Fitting a parallel slopes linear regression

<!-- 
LO: Fit a linear regression model with two explanatory variables, without interactions.
-->
In <a href="https://learn.datacamp.com/courses/introduction-to-regression-in-r">Introduction to Regression in R</a>, you learned to fit linear regression models with a single explanatory variable. In many cases, using only one explanatory variable limits the accuracy of predictions. That means that to truly master linear regression, you need to be able to include multiple explanatory variables.

The case when there is one numeric explanatory variable and one categorical explanatory variable is sometimes called a "parallel slopes" linear regression due to the shape of the predictions &ndash; more on that in the next exercise.

**Data**

Here, you'll revisit the Taiwan real estate dataset. Recall the meaning of each variable.

|Variable                |Meaning                                                             |
|:-----------------------|:-------------------------------------------------------------------|
|`dist_to_mrt_station_m` |Distance to nearest MRT metro station, in meters.                   |
|`n_convenience`         |No. of convenience stores in walking distance.                      |
|`house_age_years`       |The age of the house, in years, in 3 groups.                        |
|`price_twd_msq`         |House price per unit area, in New Taiwan dollars per meter squared. |

**Steps**

1. Using the `taiwan_real_estate` dataset, model the house price (in TWD per square meter) versus the number of nearby convenience stores.

```{r, message=FALSE}
# Load packages
library(fst)
library(tibble)

# Load data
taiwan_real_estate <- read_fst("data/taiwan_real_estate2.fst") |> as_tibble()

# Fit a linear regr'n of price_twd_msq vs. n_convenience
mdl_price_vs_conv <- lm(price_twd_msq ~ n_convenience, data = taiwan_real_estate)

# See the result
mdl_price_vs_conv
```

2. Model the house price (in TWD per square meter) versus the house age (in years). Don't include an intercept term.

```{r}
# Fit a linear regr'n of price_twd_msq vs. house_age_years, no intercept
mdl_price_vs_age <- lm(price_twd_msq ~ house_age_years + 0, data = taiwan_real_estate)

# See the result
mdl_price_vs_age
```

3. Model the house price (in TWD per square meter) versus the number of nearby convenience stores plus the house age (in years). Don't include an intercept term.

```{r}
# Fit a linear regr'n of price_twd_msq vs. n_convenience 
# plus house_age_years, no intercept
mdl_price_vs_both <- lm(price_twd_msq ~ n_convenience + house_age_years + 0, data = taiwan_real_estate)

# See the result
mdl_price_vs_both
```

Perfect parallel slopes fitting! To combine multiple explanatory variables in the regression formula, separate them with a `+`.

## Interpreting parallel slopes coefficients

<!-- 
LO: Interpret the coefficients (one slope and several intercepts) for a parallel slopes regression.
-->
For linear regression with a single numeric explanatory variable, there is an intercept coefficient and a slope coefficient. For linear regression with a single categorical explanatory variable, there is an intercept coefficient for each category.

In the "parallel slopes" case, where you have a numeric and a categorical explanatory variable, what do the coefficients mean?

> *Question*
> ---
> Look at the coefficients of `mdl_price_vs_both`. What is the meaning of the `n_convenience` coefficient?<br>
> <br>
> ⬜ The mean number of nearby convenience stores is `0.79`.<br>
> ⬜ For a house with zero nearby convenience stores, the expected house price is `0.79` TWD per square meter.<br>
> ✅ For each additional nearby convenience store, the expected house price, in TWD per square meter, increases by `0.79`.<br>
> ⬜ For each additional `0.79` nearby conveniences stores, the expected house price increases by 1 TWD per square meter.<br>

> *Question*
> ---
> What is the meaning of the `"0 to 15 years"` coefficient?<br>
> <br>
> ⬜ For a house aged 0 to 15 years, the mean number of nearby convenience stores is `9.41`.<br>
> ✅ For a house aged 0 to 15 years with zero nearby convenience stores, the expected house price is `9.41` TWD per square meter.<br>
> ⬜ For each additional year of house age, the expected house price, in TWD per square meter, increases by `9.41`.<br>
> ⬜ For each additional `15` years of house age, the expected house price increases by `9.41` TWD per square meter.<br>

Insightful interpretation! The model has one slope coefficient, and three intercept coefficients (one for each possible value of the categorical explanatory variable).

## Visualizing each explanatory variable

<!-- 
LO: Draw a scatter plot with a linear regession trend line, and a box plot.
-->
Being able to see the predictions made by a model makes it easier to understand. In the case where there is only one explanatory variable, ggplot lets you do this without any manual calculation or messing about.

To visualize the relationship between a numeric explanatory variable and the numeric response, you can draw a scatter plot with a linear trend line.

To visualize the relationship between a categorical explanatory variable and the numeric response, you can draw a box plot.

**Steps**

1. Using the `taiwan_real_estate` dataset, plot the house price versus the number of nearby convenience stores.
2. Make it a scatter plot.
3. Add a smooth linear regression trend line without a standard error ribbon.

```{r}
# Load packages
library(ggplot2)

# Using taiwan_real_estate, plot price_twd_msq vs. n_convenience
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq)) +
  # Add a point layer
  geom_point() +
  # Add a smooth trend line using linear regr'n, no ribbon
  geom_smooth(method = "lm", se = FALSE)
```

4. Using the `taiwan_real_estate` dataset, plot the house price versus the house age.
5. Make it a box plot.

```{r}
# Using taiwan_real_estate, plot price_twd_msq vs. house_age_years
ggplot(taiwan_real_estate, aes(house_age_years, price_twd_msq)) +
  # Add a box plot layer
  geom_boxplot()
```

Veracious visualizing! With a single numeric explanatory variable, the predictions form a single straight line. With a single categorical explanatory variable, the predictions are the means for each category.

## Visualizing parallel slopes

<!-- 
LO: Draw a scatter plot with parallel slope predictions.
-->
The two plots in the previous exercise gave very different predictions: one gave a predicted response that increased linearly with a numeric variable; the other gave a fixed response for each category. 
The only sensible way to reconcile these two conflicting predictions is to incorporate both explanatory variables in the model at once.

When it comes to a linear regression model with a numeric and a categorical explanatory variable, `ggplot2` doesn't have an easy, "out of the box" way to show the predictions. Fortunately, the `moderndive` package includes an extra geom, `geom_parallel_slopes()` to make it simple.

**Steps**

1. Using the `taiwan_real_estate` dataset, plot house prices versus the number of nearby convenience stores, colored by house age.
2. Make it a scatter plot.
3. Add parallel slopes, without a standard error ribbon.

```{r}
# Load packages
library(moderndive)

# Using taiwan_real_estate, plot price_twd_msq vs. n_convenience
# colored by house_age_years
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq, color = house_age_years)) +
  # Add a point layer
  geom_point() +
  # Add parallel slopes, no ribbon
  geom_parallel_slopes(se = FALSE)
```

Parallel slope paradise! The "parallel slope" model name comes from the fact that the prediction for each category is a slope, and all those slopes are parallel.

## Predicting parallel slopes

Theory. Coming soon ...


**1. Predicting parallel slopes**

Predicting responses is perhaps the most useful feature of regression models. With two explanatory variables, the code for prediction has one subtle difference from the case with a single explanatory variable.

**2. The prediction workflow 1**

The prediction workflow starts with choosing values for explanatory variables. You pick any values you want, and store them in a data frame or tibble.For a single explanatory variable, the data frame has one column. Here, it's a sequence of lengths from 5cm to 60cm, in steps of 5cm.For multiple explanatory variables, it's the same process, but there's a useful trick. expand_grid from the tidyr package returns a data frame of all combinations of its inputs.Here, you have 5cm and each fish species, 10cm and each fish species, through to 60cm and each fish species.

**3. The prediction workflow 2**

Next you add a column of predictions to the data frame. To calculate the predictions, call predict, passing the model and the explanatory data. Here's the code for one explanatory variable.With two or more explanatory variables, other than the model variable name, the code is exactly the same!

**4. Visualizing the predictions**

Just as in the single explanatory variable case, we can visualize the predictions from the model by adding another geom_point layer and setting the data argument to prediction_data.I set the size and shape arguments to make the predictions big square points. A good sign that this worked is that the prediction points lie along the lines calculated by ggplot.

**5. Manually calculating predictions**

In the previous course, you saw how to manually calculate the predictions. The coefficients function extracts the coefficients from the model.The intercept is the first coefficient, and the slope is the second coefficient.Then the response value is the intercept plus the slope times the explanatory variable.

**6. Coefficients for parallel slopes**

For the parallel slopes model, there is an added complication. Each category of the categorical variable has a different intercept. Due to the way we specified the model, the slope coefficient is the first one.

**7. Choosing an intercept with ifelse()**

You can choose this intercept using if-else statements, but this becomes clunky when you have lots of categories.With just four categories, these nested calls to ifelse are hard to write and hard to read. It's a recipe for buggy code.

**8. case_when()**

dplyr has a function called case_when that simplifies the code. Each argument to case_when is a formula, just like you use when specifying a model. On the left-hand side, you have a logical condition. On the right-hand side, you have the value to give to those rows where the condition is met.This is very abstract, so let's look at how we use it for predictions.

**9. Choosing an intercept with case_when()**

The first argument to case_when has a logical condition to check for rows where the species is Bream. On the right-hand side of the formula, we give those rows the value of the bream intercept.Then we repeat this for the other species.  This code does the same thing as the ifelse code, but I find it easier to write and to read.

**10. The final prediction step**

The final step is to calculate the response. As before, the response is the intercept plus the slope times the numeric explanatory variable. This time, the intercept is different for different rows.

**11. Compare to predict()**

The model predicts some negative masses, which isn't a good sign. Let's check that we got the right answer by calling predict.You can see that the predictions are the same numbers as the mass column that we calculated, so our calculations are correct. It's just that this model performs poorly for small fish lengths.

**12. Let's practice!**

Time for you to make predictions.

## Predicting with a parallel slopes model

<!-- 
LO: Make predictions with multiple explanatory variables.
-->
While ggplot can automatically show you model predictions, in order to get those values to program with, you'll need to do the calculations yourself. 

Just as with the case of a single explanatory variable, the workflow has two steps: create a data frame of explanatory variables, then add a column of predictions. To make sure you've got the right answer, you can add your predictions to the ggplot with the `geom_parallel_slopes()` lines.

**Steps**

1. Make a grid of explanatory data, formed from combinations of the following variables.

    * `n_convenience` should take the numbers zero to ten.
    * `house_age_years` should take the unique values of the `house_age_years` column of `taiwan_real_estate`.

```{r}
# Load packages
library(tidyr)

# Make a grid of explanatory data
explanatory_data <- expand_grid(
  # Set n_convenience to zero to ten
  n_convenience = 0:10,
  # Set house_age_years to the unique values of that variable
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# See the result
explanatory_data
```

2. Add a column to the `explanatory_data` named for the response variable, assigning to `prediction_data`.
3. The response column contain predictions made using `mdl_price_vs_both` and `explanatory_data`.

```{r}
# Load packages
library(dplyr)

# From previous step
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# Add predictions to the data frame
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_both, explanatory_data)
  )

# See the result
prediction_data
```

4. Update the plot to add a point layer of predictions. Use the `prediction_data`, set the point size to 5, and the point shape to 15.

```{r}
# From previous steps
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_both, explanatory_data)
  )

taiwan_real_estate %>% 
  ggplot(aes(n_convenience, price_twd_msq, color = house_age_years)) +
  geom_point() +
  geom_parallel_slopes(se = FALSE) +
  # Add points using prediction_data, with size 5 and shape 15
  geom_point(data = prediction_data, size = 5, shape = 15)
```

Excellent grid expansion of explanatory variables! The workflow for making predictions with multiple explanatory variables is the same as for a single variable, except that you can use `expand_grid()` to get a more complete set of predictions.

## Manually calculating predictions

<!-- 
LO: Can calculate parallel slopes linear regression predictions manually.
-->
As with simple linear regression, you can manually calculate the predictions from the model coefficients. The only change for the parallel slopes case is that the intercept is different for each category of the categorical explanatory variable. That means you need to consider the case when each each category occurs separately.

**Steps**

1. Get the coefficients from `mdl_price_vs_both`, assigning to `coeffs`.
2. Assign each of the elements of `coeffs` to the appropriate variable.

```{r}
# Get the coefficients from mdl_price_vs_both
coeffs <- coefficients(mdl_price_vs_both)

# Extract the slope coefficient
slope <- coeffs[1]

# Extract the intercept coefficient for 0 to 15
intercept_0_15 <- coeffs[2]

# Extract the intercept coefficient for 15 to 30
intercept_15_30 <- coeffs[3]

# Extract the intercept coefficient for 30 to 45
intercept_30_45 <- coeffs[4]
```

3. Add columns to `explanatory_data`.

    * To choose the `intercept`, in the case when `house_age_years` is `"0 to 15"`, choose `intercept_0_15`. In the case when `house_age_years` is `"15 to 30"`, choose `intercept_15_30`. Do likewise for `"30 to 45"`.
    * Manually calculate the predictions as the intercept plus the slope times `n_convenience`.

```{r}
# From previous step
coeffs <- coefficients(mdl_price_vs_both)
slope <- coeffs[1]
intercept_0_15 <- coeffs[2]
intercept_15_30 <- coeffs[3]
intercept_30_45 <- coeffs[4]

prediction_data <- explanatory_data %>% 
  mutate(
    # Consider the 3 cases to choose the intercept
    intercept = case_when(
      house_age_years == "0 to 15" ~ intercept_0_15,
      house_age_years == "15 to 30" ~ intercept_15_30,
      house_age_years == "30 to 45" ~ intercept_30_45 
    ),
    # Manually calculate the predictions
    price_twd_msq = intercept + slope * n_convenience
  )

# See the results
prediction_data
```

This is the case when you've succeeded! Even for more complicated linear regression models, the prediction just involve adding and multiplying coefficients for different cases of explanatory variables.

## Assessing model performance

Theory. Coming soon ...


**1. Assessing model performance**

The big benefit of using more than one explanatory variable in a model is that you can sometimes get a better fit than when you use a single explanatory variable.

**2. Model performance metrics**

In the last course, you saw two metrics for measuring model performance: the coefficient of determination, and the residual standard error. The coefficient of determination, sometimes known as the R-squared value, measures how good the regression's prediction line fits the observed values, and a larger number is better.The residual standard error, sometimes abbreviated to RSE, is - loosely speaking - the typical size of the residuals.Let's see if these metrics improve when both explanatory variables are included in the fish model.

**3. Getting the coefficient of determination**

To easily get the coefficient of determination, load the dplyr and broom packages.Recall that broom's glance function retrieves model-level metrics as a data frame. Then dplyr's pull function can be used to extract the metric we want. The coefficient of determination is called r-dot-squared. For the mass versus species model, the coefficient of determination is zero-point-seven-two, where zero is the worst possible fit and one is a perfect fit.For the mass versus length model, the coefficient of determination is better, at zero-point-eight-two.For the mass versus both model, the coefficient of determination is even higher, at zero-point-nine-seven.Using this metric, the model with both explanatory variables is the best one, since it has the highest coefficient of determination.

**4. Adjusted coefficient of determination**

Adding more explanatory variables often increases the coefficient of determination for a model, but there is a problem.Including too many explanatory variables in your model can lead to a phenomenon called overfitting. That's when your model is optimized to provide the best fit for that particular dataset, but no longer reflects the general population. In this case, the model would be overfit if it performed well on this fish dataset, but badly on a different fish dataset.A variant metric called adjusted coefficient of determination includes a small penalty term for each additional explanatory variable to compensate for this effect. Its a better metric than the plain coefficient of determination.Its equation is based on the plain coefficient of determination, the number of observations, and the number of explanatory variables, including interactions.The penalty is big enough to worry about if the plain coefficient of determination is small, or if the number of explanatory variables is a sizable fraction of the number of observations.To get this metric, we retrieve the adj-dot-r-dot-squared element from the glanced model.

**5. Getting the adjusted coefficient of determination**

To see the effect of penalization, let's look at the unadjusted and adjusted coefficients side-by-side.Since each model only contains one or two explanatory variables, the effect is tiny.

**6. Getting the residual standard error**

The code to get the residual standard error is the same as before, but this time we pull out sigma. The mass versus species model has an RSE of just over three hundred.The mass versus length model has an RSE of about one hundred and fifty. Finally, the mass versus both model has an RSE of just over one hundred, meaning that it typically gets the mass wrong by about one hundred grams. Since that number is the lowest of the three, by this metric, the mass versus both model is best.That means that all metrics indicate that the model with two explanatory variables is better than the models with just one explanatory variable.

**7. Let's practice!**

Your turn to look at some model metrics.

## Comparing coefficients of determination

<!-- 
LO: Understands that using multiple explanatory variables can increase the coefficient of determination, providing a better fit.
-->
Recall that the coefficient of determination is a measure of how well the linear regression line fits the observed values. An important motivation for including several explanatory variables in a linear regression is that you can improve the fit compared to considering only a single explanatory variable.

Here you'll compare the coefficient of determination for the three Taiwan house price models, to see which gives the best result.

**Steps**

1. Get the unadjusted and adjusted coefficients of determination for `mdl_price_vs_conv` by glancing at the model, then selecting the `r.squared` and `adj.r.squared` values. 
2. Do the same for `mdl_price_vs_age` and `mdl_price_vs_both`.

```{r}
# Load packages
library(broom)

mdl_price_vs_conv %>% 
  # Get the model-level coefficients
  glance() %>% 
  # Select the coeffs of determination
  select(r.squared, adj.r.squared)

# Get the coeffs of determination for mdl_price_vs_age
mdl_price_vs_age %>% 
  glance() %>% 
  select(r.squared, adj.r.squared)

# Get the coeffs of determination for mdl_price_vs_both
mdl_price_vs_both %>% 
  glance() %>% 
  select(r.squared, adj.r.squared)
```

> *Question*
> ---
> Which model does the adjusted coefficient of determination suggest gives a better fit?<br>
> <br>
> ⬜ `mdl_price_vs_conv`.<br>
> ⬜ `mdl_price_vs_age`.<br>
> ✅ `mdl_price_vs_both`.<br>
> ⬜ All models are equally good.<br>
> ⬜ Adjusted coefficient of determination doesn't measure the goodness of fit of a regression model.<br>

Magnificent model metric comparing! When both explanatory variables are included in the model, the adjusted coefficient of determination is higher, resulting in a better fit.

## Comparing residual standard error

<!--
LO: Understand that including multiple explanatory variables can lower the RSE, improving the fit.
-->
The other common metric for assessing model fit is the residual standard error (RSE), which measures the typical size of the residuals. 

In the last exercise you saw how including both explanatory variables into the model increased the coefficient of determination. How do you think using both explanatory variables will change the RSE?

**Steps**

1. Get the residual standard error for `mdl_price_vs_conv` by glancing at the model, then pulling the `sigma` value. 
2. Do the same for `mdl_price_vs_age`. 
3. Do the same for `mdl_price_vs_both`.

```{r}
mdl_price_vs_conv %>% 
  # Get the model-level coefficients
  glance() %>% 
  # Pull out the RSE
  pull(sigma)

# Get the RSE for mdl_price_vs_age
mdl_price_vs_age %>% 
  glance() %>% 
  pull(sigma)

# Get the RSE for mdl_price_vs_both
mdl_price_vs_both %>% 
  glance() %>% 
  pull(sigma)
```

> *Question*
> ---
> Which model does the RSE suggest gives more accurate predictions?<br>
> <br>
> ⬜ `mdl_price_vs_conv`.<br>
> ⬜ `mdl_price_vs_age`.<br>
> ✅ `mdl_price_vs_both`.<br>
> ⬜ Both models are equally good.<br>
> ⬜ RSE doesn't measure the accuracy of a regression model.<br>

Resplendent use of RSE! By including both explanatory variables in the model, a lower RSE was achieved, indicating a smaller difference between the predicted responses and the actual responses.

# 2. Interactions

Explore the effect of interactions between explanatory variables. Considering interactions allows for more realistic models that can have better predictive power. You&#39;ll also deal with Simpson&#39;s Paradox: a non-intuitive result that arises when you have multiple explanatory variables.

## Models for each category

Theory. Coming soon ...


**1. Models for each category**

The parallel slopes model enforced a common slope for each category. That's not always the best option.

**2. 4 categories**

Recall that the fish dataset had four different species of fish.One way to give each species a different slope is to run a separate model for each of these.

**3. Splitting the dataset**

There are many smart ways of splitting a dataset into parts and computing on each part. In base-R, you can use split and lapply. In dplyr, you can use nest_by and mutate.We aren't going to do that. Instead, let's filter for each species one at a time and assign the result to individual variables.I'm choosing this approach partly because I don't want fancy code to get in the way of reasoning about models, and partly because running regression models is such a fundamental task that you need to be able to write the code without thinking, and that takes practice. With this approach, you'll be writing the modeling code for every category in the dataset.

**4. 4 models**

Now we have four datasets, we can run four models. Again, there's no fancy looping, just the same model four times.Observe that each model gives a different intercept and a different slope.

**5. Explanatory data**

To make predictions with these models, we first have to create a data frame of explanatory variables. The good news is that since each model has the same explanatory variable, you only have to write this code once. Give your copy and paste fingers a rest.

**6. Making predictions**

Predicting follows the now familiar flow. Add a column with mutate, name it after the response variable, call predict with the model and the explanatory data. The only difference in each case is the model variable.It isn't necessary for calculating the predictions, but to make the plotting code you are about to see easier, I've also included the species in each prediction dataset.

**7. Visualizing predictions**

Here's the standard ggplot for showing linear regression predictions. geom_point makes it a scatter plot, and geom_smooth with method equals "lm" provides prediction lines.Unlike the parallel slopes case, each line has its own slope. This is achieved by setting the color aesthetic.

**8. Adding in your predictions**

To sanity check our predictions, we add them to the plot to see if they align with what ggplot2 calculated. The size and shape are changed to help them stand out. Thankfully, each line of squares follows ggplot's trend lines.

**9. Coefficient of determination**

An important question here is are these models better?The coefficient of determination for a model on the whole fish dataset is point-nine-two.Now here's the coefficient of determination for each of the individual models. The pike number is higher, indicating a better fit, though the numbers for the other models are lower.

**10. Residual standard error**

Here's the residual standard error for the whole dataset model, one hundred and three.For the individual models, this time the pike residual standard error is higher, indicating larger differences between actual and predicted values, but the other models show an improvement over the whole dataset model.This mixed performance result is quite common: the whole dataset model benefits from the increased power of more rows of data, whereas individual models benefit from not having to satisfy differing components of data.

**11. Let's practice!**

Let's try this on the housing data.

## One model per category

<!-- 
LO: {{This is mostly setup for the next two exercises, so not much learning.}} Can create a model for each category of a categorical variable.
-->
The model you ran on the whole dataset fits some parts of the data better than others. It's worth taking a look at what happens when you run a linear model on different parts of the dataset separately, to see if each model agrees or disagrees with the others.

**Steps**

1. Filter `taiwan_real_estate` for rows where `house_age_years` is `"0 to 15"`, assigning to `taiwan_0_to_15`. 
2. Repeat this for the `"15 to 30"` and `"30 to 45"` house age categories.

```{r}
# Filter for rows where house age is 0 to 15 years
taiwan_0_to_15 <- taiwan_real_estate %>%
  filter(house_age_years == "0 to 15")

# Filter for rows where house age is 15 to 30 years
taiwan_15_to_30 <- taiwan_real_estate %>%
  filter(house_age_years == "15 to 30")

# Filter for rows where house age is 30 to 45 years
taiwan_30_to_45 <- taiwan_real_estate %>%
  filter(house_age_years == "30 to 45")
```

3. Run a linear regression of `price_twd_msq` versus `n_convenience` using the `taiwan_0_to_15` dataset.
4. Repeat this for `taiwan_15_to_30` and `taiwan_30_to_45`.

```{r}
# From previous step
taiwan_0_to_15 <- taiwan_real_estate %>%
  filter(house_age_years == "0 to 15")
taiwan_15_to_30 <- taiwan_real_estate %>%
  filter(house_age_years == "15 to 30")
taiwan_30_to_45 <- taiwan_real_estate %>%
  filter(house_age_years == "30 to 45")

# Model price vs. no. convenience stores using 0 to 15 data
mdl_0_to_15 <- lm(price_twd_msq ~ n_convenience, data = taiwan_0_to_15)

# Model price vs. no. convenience stores using 15 to 30 data
mdl_15_to_30 <- lm(price_twd_msq ~ n_convenience, data = taiwan_15_to_30)

# Model price vs. no. convenience stores using 30 to 45 data
mdl_30_to_45 <- lm(price_twd_msq ~ n_convenience, data = taiwan_30_to_45)

# See the results
mdl_0_to_15
mdl_15_to_30
mdl_30_to_45
```

Multiple model magic! You now have three models giving three different answers. Let's visualize the differences, then figure out how to reconcile the differences.

## Predicting multiple models

<!-- 
LO: Make predictions from a simple linear regression model.
-->
In order to see what each of the models for individual categories are doing, it's helpful to make predictions from them. The flow is exactly the same as the flow for making predictions on the whole model, though remember that you only have a single explanatory variable in these models (so `expand_grid()` isn't needed.)

**Steps**

1. Create a tibble of explanatory data, setting `n_convenience` to a vector from zero to ten, assigning to `explanatory_data_0_to_15`.

```{r}
# Create a tibble of explanatory data, setting
# no. of conv stores to 0 to 10
explanatory_data <- tibble(
  n_convenience = 0:10
)
```

2. Add a column of predictions named `price_twd_msq` to `explanatory_data`, using `mdl_0_to_15` and `explanatory_data`. Assign to `prediction_data_0_to_15`.
3. Repeat this for the 15 to 30 year and 30 to 45 year house age categories.

```{r}
# From previous step
explanatory_data <- tibble(
  n_convenience = 0:10
)

# Add column of predictions using "0 to 15" model and explanatory data 
prediction_data_0_to_15 <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_0_to_15, explanatory_data))

# Same again, with "15 to 30"
prediction_data_15_to_30 <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_15_to_30, explanatory_data))

# Same again, with "30 to 45"
prediction_data_30_to_45 <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_30_to_45, explanatory_data))
```

Predictions, predictions, predictions! Now that you have predictions for each model, let's see how they compare.

## Visualizing multiple models

<!-- 
LO: Visualize multiple linear models on the same plot.
-->
In the last two exercises, you ran models for each category of house ages separately, then calculated predictions for each model. Now it's time to visualize those predictions to see how they compare.

When you use `geom_smooth()` in a ggplot with an aesthetic that splits the dataset into groups and draws a line for each group (like the `color` aesthetic), you get multiple trend lines. This is the same as running a model on each group separately, so we get a chance to test our predictions against ggplot's.

**Steps**

1. Using `taiwan_real_estate`, plot `price_twd_msq` versus `n_convenience` colored by `house_age_years`.
2. Add a point layer.
3. Add smooth trend lines for each color using the linear regression method and turning off the standard error ribbon.

```{r}
# Using taiwan_real_estate, plot price vs. no. of conv. stores
# colored by house age
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq, color = house_age_years)) +
  # Make it a scatter plot
  geom_point() +
  # Add smooth linear regression trend lines, no ribbon
  geom_smooth(method = "lm", se = FALSE)
```

4. Extend the plot by adding the prediction points from `prediction_data_0_to_15`. Color them red, with size 3 and shape 15.
5. Add prediction points from `prediction_data_15_to_30`, colored green, size 4, and shape 15.
6. Add prediction points from `prediction_data_30_to_45`, colored blue, size 4, and shape 15.

```{r}
# Extend the plot to include prediction points
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq, color = house_age_years)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  # Add points using prediction_data_0_to_15, colored red, size 3, shape 15
  geom_point(data = prediction_data_0_to_15, color = "red", size = 3, shape = 15) +
  # Add points using prediction_data_15_to_30, colored green, size 3, shape 15
  geom_point(data = prediction_data_15_to_30, color = "green", size = 3, shape = 15) +
  # Add points using prediction_data_30_to_45, colored blue, size 3, shape 15
  geom_point(data = prediction_data_30_to_45, color = "blue", size = 3, shape = 15)
```

Three cheers for three trend lines! It's a good sign that our predictions match those of ggplot's. Notice that the 30 to 45 year house age group has a much shallower slope compared to the other lines.

## Assessing model performance

<!-- 
LO: Can extract model performance metrics from linear regression model objects.
-->
To test which approach is best &ndash; the whole dataset model or the models for each house age category &ndash; you need to calculate some metrics. Here's, you'll compare the coefficient of determination and the residual standard error for each model.

**Steps**

1. Get the coefficient of determination for `mdl_all_ages`, `mdl_0_to_15`, `mdl_15_to_30`, and `mdl_30_to_45`.

```{r}
# Get the coeff. of determination for mdl_all_ages
mdl_all_ages <- lm(price_twd_msq ~ n_convenience, data = taiwan_real_estate)
mdl_all_ages %>% 
  glance() %>% 
  pull(r.squared)

# Get the coeff. of determination for mdl_0_to_15
mdl_0_to_15 %>% 
  glance() %>% 
  pull(r.squared)

# Get the coeff. of determination for mdl_15_to_30
mdl_15_to_30 %>% 
  glance() %>% 
  pull(r.squared)

# Get the coeff. of determination for mdl_30_to_45
mdl_30_to_45 %>% 
  glance() %>% 
  pull(r.squared)
```

2. Get the residual standard error for `mdl_all_ages`, `mdl_0_to_15`, `mdl_15_to_30`, and `mdl_30_to_45`.

```{r}
# Get the RSE for mdl_all
mdl_all_ages %>% 
  glance() %>% 
  pull(sigma)

# Get the RSE for mdl_0_to_15
mdl_0_to_15 %>% 
  glance() %>% 
  pull(sigma)

# Get the RSE for mdl_15_to_30
mdl_15_to_30 %>% 
  glance() %>% 
  pull(sigma)

# Get the RSE for mdl_30_to_45
mdl_30_to_45 %>% 
  glance() %>% 
  pull(sigma)
```

Meritable model metric measuring! It seems that both metrics for the 15 to 30 age group model are much better than those for the whole dataset model, but the models for the other two age groups are similar to the whole dataset model. Thus using individual models will improve predictions for 15 to 30 age group.

## One model with an interaction

Theory. Coming soon ...

**1. One model with an interaction**

Messing about with different models for different bits of your dataset is a pain. A better solution is to specify a single model that contains intercepts and slopes for each category. This is achieved through specifying interactions between explanatory variables.

**2. What is an interaction?**

To understand the idea of interactions between explanatory variables, consider what we know about the fish dataset. Different fish species have different mass to length ratios. In statistical terms, we can say that the effect that length has on the expected mass of the fish varies between species. That means that length and species interact. More generally, if the effect of one explanatory variable on the expected response has different values dependent on the values of another explanatory variable, then those two explanatory variables interact.

**3. Specifying interactions**

You've seen how to include multiple explanatory variables in a formula using plus, for example, length plus species. To include an interaction between those variables, you simply swap the plus for a times. I've called this syntax implicit because you didn't write down what interactions are needed - R figures that out itself. Usually this concise syntax is best, but occasionally you may wish to explicitly document which interactions are included in the model. The explicit syntax is to add each explanatory variable separated by plus then add a third term with both explanatory variables separated by a colon. The result is exactly the same, so choosing a syntax depends on personal preference: do you like brevity or detail?

**4. Running the model**

Here's the formula in a model, and there are eight coefficients. As you saw in the models with a categorical explanatory variable, the coefficients are tricky to understand. The (Intercept) coefficient is the intercept for the first species, namely bream. The length coefficient is the slope for the bream. Then the (Intercept) coefficient plus the speciesPerch coefficient is the intercept for perch. And the length coefficient plus the length-speciesPerch coefficient is the slope for perch. It's a mess.

**5. Easier to understand coefficients**

Ironically, to get easier to understand coefficients, we need to make the formula harder to read. This is the same model specified differently. On the right-hand side of the formula, you can see the categorical explanatory variable, species, then an interaction between the two explanatory variables, then zero to remove the global intercept. Now we get an intercept coefficient for each species shown on the top row, and a slope coefficient for each species shown on the bottom row.

**6. Familiar numbers**

You've seen all these coefficient values before. If we examine the coefficients from the model on the bream data from the previous video, you can see that the intercept and slope are the same as in the model we just made. The same is true for the other three species. In fact, the model with the interaction is effectively the same as fitting separate models for each category, only you get the convenience of not having to manage four sets of code.

**7. Let's practice!**

Let's interact with some exercises! 

## Specifying an interaction

<!--
LO: Specify an interaction term in a model formula using times or colon syntax.
-->
So far you used a single parallel slopes model, which gave an OK fit for the whole dataset, then three separate models for each house age category, which gave a better fit for each individual category, but was clunky because you had three separate models to work with and explain. Ideally, you'd have a single model that had all the predictive power of the individual models.

Defining this single model is achieved through adding interactions between explanatory variables. R's formula syntax is flexible, and gives you a couple of options, depending on whether you prefer concise code that is quick to  type and to read, or explicit code that describes what you are doing in detail.

**Steps**

1. Fit a linear regression of `price_twd_msq` versus `n_convenience` and `house_age_years` and their interaction, using the "times" syntax to implicitly generate the interaction.

```{r}
# Model price vs both with an interaction using "times" syntax
lm(price_twd_msq ~ n_convenience * house_age_years, data = taiwan_real_estate)
```

2. Fit a linear regression of `price_twd_msq` versus `n_convenience` and `house_age_years` and their interaction, using the "colon" syntax to explicitly generate the interaction.

```{r}
# Model price vs both with an interaction using "colon" syntax
lm(
  price_twd_msq ~ n_convenience + house_age_years + n_convenience:house_age_years, 
  data = taiwan_real_estate
)
```

Incredible interaction inclusion! Notice that the model coefficients are the same in each case. The formula syntax is very flexible, giving you a choice between 'easy to type' and 'very explicit', depending on your preference.

## Interactions with understandable coeffs

<!--
LO: Can interpret coefficients of linear regression with interaction term.
-->
The previous model with the interaction term returned coefficients that were a little tricky to interpret. In order clarify what the model is predicting, you can reformulate the model in a way that returns understandable coefficients. For further clarity, you can compare the results to the models on the separate house age categories (`mdl_0_to_15`, `mdl_15_to_30`, and `mdl_30_to_45`).

`taiwan_real_estate`, `mdl_0_to_15`, `mdl_15_to_30`, and `mdl_30_to_45` are available.

**Steps**

1. Fit a linear regression of `price_twd_msq` versus `house_age_years` plus an interaction between `n_convenience` and `house_age_years`, and no global intercept, using the `taiwan_real_estate` dataset.
2. For comparison, get the coefficients for the three models for each category: `mdl_0_to_15`, `mdl_15_to_30`, and `mdl_30_to_45`.

```{r}
# Model price vs. house age plus an interaction, no intercept
mdl_readable_inter <- lm(
  price_twd_msq ~ house_age_years + n_convenience:house_age_years + 0, 
  data = taiwan_real_estate
)

# See the result
mdl_readable_inter

# Get coefficients for mdl_0_to_15
coefficients(mdl_0_to_15)

# Get coefficients for mdl_15_to_30
coefficients(mdl_15_to_30)

# Get coefficients for mdl_30_to_45
coefficients(mdl_30_to_45)
```

> *Question*
> ---
> Which statement about the coefficients of `mdl_readable_inter` is true?<br>
> <br>
> ⬜ For house ages of 0 to 15 years, when there are zero nearby convenience stores, the expected house price is 6.87 TWD per square meter.<br>
> ✅ The expected increase in house price for each nearby convenience store is lowest for the 30 to 45 year age group.<br>
> ⬜ The expected increase in house price for each nearby convenience store is lowest for the 15 to 30 year age group.<br>
> ⬜ For house ages of 0 to 15 years, when there are zero nearby convenience stores, the expected house price is 0.83 TWD per square meter.<br>

Insightful intercept introspection! Sometimes fiddling about with how the model formula is specified makes it easier to interpret the coefficients. In this version, you can see how each category has its own intercept and slope (just like the 3 separate models had).

## Making predictions with interactions

Theory. Coming soon ...


**1. Making predictions with interactions**

Let's run through the prediction flow again, this time with the model containing an interaction. I'm hoping you get a sense of deja vu, because you've seen all this code before.

**2. The model with the interaction**

Here's the model of mass versus length and species with an interaction. I've used the version of the formula that gives the clearest coefficients, but I could also have used the simpler length times species syntax; it doesn't affect the predictions.

**3. The prediction flow, again**

Here's the code for the prediction flow. It's exactly the same as the code in the parallel slopes model. R will automatically take care of the interaction, so you don't need to change anything.I love it when things just work!The only thing to remember here is the use of expand_grid to get all the combinations of lengths and species.

**4. Visualizing the predictions**

Here's the plot of predictions - both ggplot's automatic prediction lines, and the points predicted from the model.The plot is identical to the one you saw in the first video of this chapter. This time however, the code is simpler because you don't have four separate models to worry about.

**5. Manually calculating the predictions**

To see how the predictions work, let's manually calculate them. First we get the coefficients from the model using the coefficients function.Then we use square bracket indexing to extract the four intercepts and the four slopes.

**6. Manually calculating the predictions**

As usual, we add the predictions as a column named after the response variable. We'll use case_when to distinguish the calculations for each species.

**7. Manually calculating the predictions**

The first species is bream, so the left-hand side of the first formula filters for rows containing bream.

**8. Manually calculating the predictions**

The calculation for the bream rows is the bream-specific intercept plus the bream-specific slope times the length of the fish.

**9. Manually calculating the predictions**

Then we repeat this for the other three species. The calculated values are the same as those returned by the predict function.

**10. Let's practice!**

Time for you to make some predictions.

## Predicting with interactions

<!-- 
LO: Make predictions of a linear regression containing interaction terms.
-->
As with every other regression model you've created, the fun part is making predictions. Fortunately, the code flow for this case is the same as the one without interactions &ndash; R can handle calculating the interactions without any extra prompting from you. The only thing you need to remember is the trick for getting combinations of explanatory variables.

**Steps**

1. Make a grid of explanatory data, formed from combinations of the following variables.

    * `n_convenience` should take the numbers zero to ten.
    * `house_age_years` should take the unique values of the `house_age_years` column of `taiwan_real_estate`.

```{r}
# Make a grid of explanatory data
explanatory_data <- expand_grid(
  # Set n_convenience to zero to ten
  n_convenience = 0:10,
  # Set house_age_years to the unique values of that variable
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# See the result
explanatory_data
```

2. Add a column to the `explanatory_data`, assigning to `prediction_data`.
3. The column should be named after the response variable, and contain predictions made using `mdl_price_vs_both_inter` and `explanatory_data`.

```{r}
# From previous step
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# Add predictions to the data frame
mdl_price_vs_both_inter <- mdl_readable_inter

prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_both_inter, explanatory_data)
  )

# See the result
prediction_data
```

4. Using `taiwan_real_estate`, plot `price_twd_msq` versus `n_convenience`, colored by `house_age_years`.
5. Add a point layer.
6. Add smooth trend lines using linear regression, no standard error ribbon.
7. Add another point layer using `prediction_data`, with `size` `5` and `shape` `15`.

```{r}
# From previous step
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_both_inter, explanatory_data)
  )

# Using taiwan_real_estate, plot price vs. no. of convenience 
# stores, colored by house age
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq, color = house_age_years)) +
  # Make it a scatter plot
  geom_point() +
  # Add linear regression trend lines, no ribbon
  geom_smooth(method = "lm", se = FALSE) +
  # Add points from prediction_data, size 5, shape 15
  geom_point(data = prediction_data, size = 5, shape = 15)
```

Pretty good predicting! The code flow here is identical to the parallel slopes case in Chapter 1. Take a look again at that exercise to check.

## Manually calculating predictions with interactions

<!--
LO: Can calculate linear regression predictions manually for case with 1 numeric, 1 categorical explanatory variable and interaction.
-->
In order to understand how `predict()` works, it's time to calculate the predictions manually again. For this model, there are three separate lines to calculate for, and in each one, the prediction is an intercept plus a slope times the numeric explanatory value. The tricky part is getting the right intercept and the right slope for each case.

**Steps**

1. Get the coefficients from `mdl_price_vs_both_inter`, assigning to `coeffs`.
2. Get the three intercept coefficients from `coeffs`, assigning to `intercept_0_15`, `intercept_15_30`, and `intercept_30_45`.
3. Get the three slope coefficients from `coeffs`, assigning to `slope_0_15`, `slope_15_30`, and `slope_30_45`.

```{r}
# Get the coefficients from mdl_price_vs_both_inter
coeffs <- coefficients(mdl_price_vs_both_inter)

# Get the intercept for 0 to 15 year age group
intercept_0_15 <- coeffs[1]

# Get the intercept for 15 to 30 year age group
intercept_15_30 <- coeffs[2]

# Get the intercept for 30 to 45 year age group
intercept_30_45 <- coeffs[3]

# Get the slope for 0 to 15 year age group
slope_0_15 <- coeffs[4]

# Get the slope for 15 to 30 year age group
slope_15_30 <- coeffs[5]

# Get the slope for 30 to 45 year age group
slope_30_45 <- coeffs[6]
```

4. Add a `price_twd_msq` column to `explanatory_data` containing the predictions.

    * In the case when `house_age_years` is `"0 to 15"`, choose the appropriate intercept plus the appropriate slope times the number of nearby convenience stores.
    * Do likewise for the cases where the house age is `"15 to 30"` and `"30 to 45"`.

```{r}
# From previous step
coeffs <- coefficients(mdl_price_vs_both_inter)
intercept_0_15 <- coeffs[1]
intercept_15_30 <- coeffs[2]
intercept_30_45 <- coeffs[3]
slope_0_15 <- coeffs[4]
slope_15_30 <- coeffs[5]
slope_30_45 <- coeffs[6]

prediction_data <- explanatory_data %>% 
  mutate(
    # Consider the 3 cases to choose the price
    price_twd_msq = case_when(
      house_age_years == "0 to 15" ~ intercept_0_15 + slope_0_15 * n_convenience,
      house_age_years == "15 to 30" ~ intercept_15_30 + slope_15_30 * n_convenience,
      house_age_years == "30 to 45" ~ intercept_30_45 + slope_30_45 * n_convenience 
    ) 
  )

# See the result
prediction_data
```

Magnificent manual prediction! As before, the prediction just involved adding and multiplying coefficients and explanatory values.

## Simpson's Paradox

Theory. Coming soon ...


**1. Simpson's Paradox**

This chapter looked at the difference between models of a whole dataset and individual models for each category. For some datasets, this can lead to a nonintuitive result known as Simpson's Paradox.

**2. A most ingenious paradox!**

Simpson's paradox is that the trend given by a model on the whole dataset is very different from the trends in subsets of the data.That's pretty abstract, so let's try an example.

**3. Synthetic Simpson data**

Here's a synthetic - that is, made up - dataset designed to demonstrate the paradox.Each row has an x and y coordinate, and the dataset is split into five groups, labeled A to E.

**4. Linear regressions**

Fitting a linear regression of y versus x to the whole dataset shows a positive slope of one point seven five.However, fitting a model that includes the group and an interaction shows something completely different.The bottom row of coefficients contains the slope for each group. Every group has a negative slope, apparently contradicting the fact that the whole dataset has a positive slope.Let's visualize the dataset to try and reconcile these opposite coefficients.

**5. Plotting the whole dataset**

This is the now-standard scatter plot with a linear regression trend line.As x increases, so does y, resulting in a positive slope over the whole dataset.

**6. Plotting by group**

Amending the plot to color the lines by group shows that within each group, y decreases as x increases.

**7. Reconciling the difference**

One moral of this story is that it's helpful to try and visualize your dataset. This is especially true if different models give conflicting results.Some common advice for how to choose which model is best is correct but annoying.It depends on the dataset and what question you are trying to answer.A useful corollary is that you should decide on a question to answer before you start fitting models.

**8. Test score example**

Thinking up examples where the grouped model is best is fairly easy. Here's the same synthetic dataset as before, with different axis labels.If x is the number of hours spent playing video games each month, and y is the score on a test, modeling the whole dataset suggests that playing more video games is related to a higher test score.If we reveal that each group represents the age of the child taking the test, it changes the interpretation. Now older children score more highly in the test, and playing lots of video games is related to a lower score.

**9. Infectious disease example**

Coming up with examples where the model of the whole dataset is more useful than the model split by group is much harder. So much so that I had to ask the internet for suggestions.One proposed example was that for an infectious disease, the infection rate tends to be higher when the population density is higher. In this plot, each point represents a neighborhood in a city.Splitting by city reveals that the highest density areas of each city have lower infection rates. However, this may be due to other things that you haven't included in the model, like the wealth and demographics of the residents.That's an interesting insight, but "increasing population density is related to increased infection rate" is arguably more important.

**10. Reconciling the difference, again**

Unfortunately, resolving these model disagreements is messy. Often, the models including the groups will contain insight that you'd miss otherwise.The disagreements between the models may reveal that you need even more explanatory variables to understand why they are different.Finally, I'm going to repeat the correct but annoying advice: to choose the best model you need contextual information about what your dataset means and what question you are trying to answer.

**11. Simpson's paradox in real datasets**

Such a clear case of Simpson's paradox is very rare. Subtle differences between models are more common. A slope may go to zero instead of changing its direction. You may only see the effect in some groups, but not all of them.

**12. Let's practice!**

Time to play with the paradox.

## Modeling eBay auctions

<!-- 
LO: Can fit and visualize a linear regression on the whole dataset.
-->
Sometimes modeling a whole dataset suggests trends that disagree with models on separate parts of that dataset. This is known as Simpson's paradox. In the most extreme case, you may see a positive slope on the whole dataset, and negative slopes on every subset of that dataset (or the other way around).

Over the next few exercises, you'll look at <a href="http://www.modelingonlineauctions.com/datasets">eBay auctions</a> of Palm Pilot M515 PDA models.

|variable       |meaning                        |
|:--------------|:------------------------------|
|`price`        |Final sale price, USD          |
|`openbid`      |The opening bid, USD           |
|`auction_type` |How long did the auction last? |

**Steps**

1. Look at the structure of the `auctions` dataset and familiarize yourself with its columns.
2. Fit a linear regression model of `price` versus `openbid`, using the `auctions` dataset. *Look at the coefficients.*

```{r}
# Load data
auctions <- readRDS("data/auctions_mod.rds")

# Take a glimpse at the dataset
glimpse(auctions)

# Model price vs. opening bid using auctions
mdl_price_vs_openbid <- lm(price ~ openbid, data = auctions)

# See the result
mdl_price_vs_openbid
```

3. Using auctions, plot `price` versus `openbid` as a scatter plot with linear regression trend lines (no ribbon). *Look at the trend line.*

```{r}
# Using auctions, plot price vs. opening bid as a 
# scatter plot with linear regression trend lines
ggplot(auctions, aes(openbid, price)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```

Amazing auction modeling! The slope coefficient is small enough that it might as well be zero. That is, opening bid appears to have no effect on the final sale price for Palm Pilots.

## Modeling each auction type

<!-- 
LO: Can reason about Simpson's Paradox.
-->
You just saw that the opening bid price appeared not to affect the final sale price of Palm Pilots in the eBay auctions. Now let's look at what happens when you model the three auction types (3 day, 5 day, and 7 day) separately.

**Steps**

1. Fit a linear regression model of `price` versus `openbid` and `auction_type`, with an interaction, using the `auctions` dataset. *Look at the coefficients.*

```{r}
# Fit linear regression of price vs. opening bid and auction 
# type, with an interaction.
mdl_price_vs_both <- lm(
  price ~ auction_type + openbid:auction_type + 0, # or price ~ auction_type * openbid
  data = auctions
)

# See the result
mdl_price_vs_both
```

2. Using `auctions`, plot `price` versus `openbid`, colored by `auction_type`, as a scatter plot with linear regression trend lines (no ribbon). *Look at the trend lines.*

```{r}
# Using auctions, plot price vs. opening bid colored by
# auction type as a scatter plot with linear regr'n trend lines
ggplot(auctions, aes(openbid, price, color = auction_type)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```

> *Question*
> ---
> Which statement about the model resolves Simpson's Paradox?<br>
> <br>
> ⬜ The model of the whole dataset showed no expected change in final sale price due to opening bid. Since this model includes all the data, we should believe it and conclude that there is no effect from opening bid price.<br>
> ⬜ The model including auction type showed that final sale price increases in the 5 day auction category. Since this model is more specific, we should believe it and conclude that there is an effect from opening bid price.<br>
> ⬜ The two models disagree, so we can't conclude anything from the models.<br>
> ✅ The two models disagree, and the best model to take advice from depends upon the question we are trying to answer.<br>

Super Simpson solving! Interpreting models is a subtle art, and your conclusions need to be based on the question you are trying to answer. Here, the answer to 'Does opening bid affect final sale price?' is *no* overall, but the answer to 'Does opening bid price affect final sale price for any type of auction?' is *yes, for 5 day auctions*.

# 3. Multiple Linear Regression

See how modeling, and linear regression in particular, makes it easy to work with more than two explanatory variables. Once you&#39;ve mastered fitting linear regression models, you&#39;ll get to implement your own linear regression algorithm.

## Two numeric explanatory variables

Theory. Coming soon ...


**1. Two numeric explanatory variables**

In the previous chapters, the models had one numeric and one categorical explanatory variable. Let's see what changes if you have two numeric explanatory variables instead.

**2. Visualizing 3 numeric variables**

Two numeric explanatory variables plus a numeric response variable gives three numeric variables to plot. Since scatter plots are designed to show relationships between two numeric variables, it takes more thinking about. There are two common choices.Either we draw a 3D scatter plot, or a 2D scatter plot, using color for the response variable.

**3. Another column for the fish dataset**

Let's revisit the fish dataset, which now has an extra numeric column, the height of the fish in centimeters.

**4. 3D scatter plot**

To make a 3D scatter plot, you can use scatter3D from the plot3D package. Unlike ggplot, this simply requires three numeric vectors for the x, y, and z coordinates.Writing the name of the dataset and using the dollar operator to access columns three times is tedious, so for code like this I prefer magrittr's dollar pipe.This code does the same thing, but is easier to write and read. The dollar pipe treats each of the arguments in the following function call as though they are columns in the data frame given to the left of the pipe.

**5. 3D scatter plot**

The code is nicer  now, but there is a bigger problem, namely that the plot is impossible to interpret.Some cleanup work like labeling axes can make things easier, but the fundamental problem is that screens are two dimensional, so 3D plots always suffer perspective issues. The only way to circumvent this is to create an interactive plot that the audience can rotate to explore the data from different angles.Maybe virtual reality will solve this one day, but for now let's move on to the next type of plot.

**6. 2D scatter plot, color for response**

The next plot type to explore uses color for the response variable. This is a standard 2D scatter plot so we can use ggplot2.It's an improvement, but ggplot's default color scale doesn't make it easy to pick out differences in blues.

**7. Viridis color scales**

ggplot has a set of color scales called "viridis" that provide easier to distinguish colors.scale_color_viridis_c is used for continuous scales, where you have numeric data.This plot uses the "inferno" palette option, which moves from black through blue and red to yellow. As you move up and to the right in the plot, the colors get brighter, representing heavier fish.

**8. Modeling with 2 numeric explanatory variables**

Although plotting was harder with this extra explanatory variable, modeling isn't much different. The explanatory variables on the right of the formula are separated with a plus, as before.You get a global intercept coefficient, and one slope coefficient for each explanatory variable.

**9. The prediction flow**

The prediction flow is no different. Create a grid of explanatory values with expand_grid, then add a column of predictions with mutate and predict.

**10. Plotting the predictions**

The plotting code also remains the same, though the results look different.The color grid gives a nice overview of how the response variable changes over the plane of the explanatory variables. The heaviest fish are in the top-right, where they are long and tall.

**11. Including an interaction**

To include an interaction in the model, the only change is to replace the plus in the formula with a times. This gives you one extra slope term for the effect of the interaction between the two explanatory variables.

**12. The prediction flow again**

The prediction flow is exactly the same as before. Boring, but pleasingly easy.

**13. Plotting the predictions**

The plotting code is identical, but the colors on the plot are slightly different.In this case, the colors of the square prediction points closely match the colors of the nearby circular data points, which is a nice visual indicator that the model is a good fit.

**14. Let's practice!**

Back to the housing dataset.

## 3D visualizations

<!-- 
LO: Can visualize three continuous variables.
-->
Since computer screens and paper are both two-dimensional objects, most plots are best suited to visualizing two variables at once. For the case of three continuous variables, you can draw a 3D scatter plot, but perspective problems usually make it difficult to interpret. There are some "flat" alternatives that provide easier interpretation, though they require a little thinking about to make.

**Steps**

1. With the `taiwan_real_estate` dataset, draw a 3D scatter plot of the number of nearby convenience stores on the x-axis, the **square-root** of the distance to the nearest MRT stop on the y-axis, and the house price on the z-axis.

```{r}
# Load packages
library(plot3D)
library(magrittr)

# With taiwan_real_estate, draw a 3D scatter plot of
# no. of conv. stores, sqrt dist to MRT, and price
taiwan_real_estate %$% 
  scatter3D(n_convenience, sqrt(dist_to_mrt_m), price_twd_msq)
```

2. With the `taiwan_real_estate` dataset, draw a scatter plot of the square-root of the distance to the nearest MRT stop versus the number of nearby convenience stores, colored by house price.
3. Use the continuous viridis color scale, using the `"plasma"` option.

```{r}
# Using taiwan_real_estate, plot sqrt dist to MRT vs. 
# no. of conv. stores, colored by price
ggplot(
  taiwan_real_estate, 
  aes(n_convenience, sqrt(dist_to_mrt_m), color = price_twd_msq)
) + 
  # Make it a scatter plot
  geom_point() +
  # Use the continuous viridis plasma color scale
  scale_color_viridis_c(option = "plasma")
```

Three cheers for plotting three variables! 3D scatter plots are usually a pain to easily interpret due to problems with perspective. The best alternative for displaying a third variable involves using colors.

## Modeling 2 numeric explanatory variables

<!--
LO: Can model and predict with two numeric explanatory variables.
-->
You already saw how to make a model and predictions with a numeric and a categorical explanatory variable. The code for modeling and predicting with two numeric explanatory variables in the same, other than a slight difference in how to specify the explanatory variables to make predictions against.

Here you'll model and predict the house prices against the number of nearby convenience stores and the square-root of the distance to the nearest MRT station.

**Steps**

1. Fit a linear regression of house price versus the number of convenience stores and the square-root of the distance to the nearest MRT stations, without an interaction, using the `taiwan_real_estate` dataset.

```{r}
# Fit a linear regression of price vs. no. of conv. stores
# and sqrt dist. to nearest MRT, no interaction
mdl_price_vs_conv_dist <- lm(
  price_twd_msq ~ n_convenience + sqrt(dist_to_mrt_m), data = taiwan_real_estate
)

# See the result
mdl_price_vs_conv_dist
```

2. Create expanded grid of explanatory variables with number of convenience stores from 0 to 10 and the distance to the nearest MRT station as a sequence from 0 to 80 in steps of 10, all squared (0, 100, 400, …, 6400). Assign to `explanatory_data`.
3. Add a column of predictions to `explanatory_data` using `mdl_price_vs_conv_dist` and `explanatory_data`. Assign to `prediction_data`.

```{r}
# From previous step 
mdl_price_vs_conv_dist <- lm(price_twd_msq ~ n_convenience + sqrt(dist_to_mrt_m), data = taiwan_real_estate)

# Create expanded grid of explanatory variables with
# no. of conv. stores and  dist. to nearest MRT
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  dist_to_mrt_m = seq(0, 80, 10) ^ 2
)

# Add predictions using mdl_price_vs_conv_dist and explanatory_data
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_conv_dist, explanatory_data)
  )

# See the result
prediction_data
```

4. Extend the plot to add a layer of points using the prediction data, colored yellow, with size 3.

```{r}
# From previous steps
mdl_price_vs_conv_dist <- lm(price_twd_msq ~ n_convenience + sqrt(dist_to_mrt_m), data = taiwan_real_estate)
explanatory_data <- expand_grid(n_convenience = 0:10, dist_to_mrt_m = seq(0, 80, 10) ^ 2)
prediction_data <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_price_vs_conv_dist, explanatory_data))

# Add predictions to plot
ggplot(
  taiwan_real_estate, 
  aes(n_convenience, sqrt(dist_to_mrt_m), color = price_twd_msq)
) + 
  geom_point() +
  scale_color_viridis_c(option = "plasma") +
  # Add prediction points colored yellow, size 3
  geom_point(data = prediction_data, color = "yellow", size = 3)
```

Nice numeric modeling! The modeling and prediction flow for two numeric variables is just as it was for the previous case that included a categorical variable. R automatically handles this different scenario.

## Including an interaction

<!--
LO: Can model and predict with two numeric explanatory variables with an interaction.
-->
Just as in the case with one numeric and one categorical explanatory variable, it is possible that numeric explanatory variables can interact. With this model structure, you'll get a third slope coefficient: one for each explanatory variable and one for the interaction. 

Here you'll run and predict the same model as in the previous exercise, but this time including an interaction between the explanatory variables.

**Steps**

1. Fit a linear regression of house price versus the number of convenience stores and the square-root of the distance to the nearest MRT stations, *with* an interaction, using the `taiwan_real_estate` dataset.

```{r}
# Fit a linear regression of price vs. no. of conv. stores
# and sqrt dist. to nearest MRT, with interaction
mdl_price_vs_conv_dist <- lm(
  price_twd_msq ~ n_convenience * sqrt(dist_to_mrt_m), data = taiwan_real_estate
)

# See the result
mdl_price_vs_conv_dist
```

2. Create expanded grid of explanatory variables with number of convenience stores from 0 to 10 and the distance to the nearest MRT station as a sequence from 0 to 80 in steps of 10, all squared (0, 100, 400, …, 6400). Assign to `explanatory_data`.
3. Add a column of predictions to `explanatory_data` using `mdl_price_vs_conv_dist` and `explanatory_data`. Assign to `prediction_data`.

```{r}
# From previous step 
mdl_price_vs_conv_dist <- lm(price_twd_msq ~ n_convenience * sqrt(dist_to_mrt_m), data = taiwan_real_estate)

# Create expanded grid of explanatory variables with
# no. of conv. stores and  dist. to nearest MRT
explanatory_data <- expand_grid(
  n_convenience = 0:10,
  dist_to_mrt_m = seq(0, 80, 10) ^ 2
)

# Add predictions using mdl_price_vs_conv_dist and explanatory_data
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_conv_dist, explanatory_data)
  )

# See the result
prediction_data
```

4. Extend the plot to add a layer of points using the prediction data, colored yellow, with size 3.

```{r}
# From previous steps
mdl_price_vs_conv_dist <- lm(price_twd_msq ~ n_convenience * sqrt(dist_to_mrt_m), data = taiwan_real_estate)
explanatory_data <- expand_grid(n_convenience = 0:10, dist_to_mrt_m = seq(0, 80, 10) ^ 2)
prediction_data <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_price_vs_conv_dist, explanatory_data))

# Add predictions to plot
ggplot(
  taiwan_real_estate, 
  aes(n_convenience, sqrt(dist_to_mrt_m), color = price_twd_msq)
) + 
  geom_point() +
  scale_color_viridis_c(option = "plasma")+
  # Add prediction points colored yellow, size 3
  geom_point(data = prediction_data, color = "yellow", size = 3)
```

Insightful interaction inclusion! Although the predictions from the model are different now that the interaction is included, the code to create them is the same.

## More than 2 explanatory variables

Theory. Coming soon ...


**1. More than 2 explanatory variables**

Regression models aren't limited to two explanatory variables. Here, we'll consider three of them, and think about what happens when you increase that number even further.

**2. From last time**

In the last video you saw this scatter plot with the response variable, mass, indicated by color, and the explanatory variables shown on the x and y axes. You can see several distinct clusters of points. Perhaps these correspond to the different species of fish.We can check this by faceting on species.

**3. Faceting by species**

Giving each species its own panel with facet_wrap makes the groups of data more easily apparent.There is a noticeable strong positive correlation between length and height for each species of fish. The relationship between the explanatory variables and the response is harder to quantify because you can't determine colors as accurately as x and y coordinates, but for each species you can see that as fish get longer and taller they also get heavier.In general, while it is tricky to include more than three numeric variables in a scatter plot, you can include as many categorical variables as you like using faceting. However, more facets can make it harder to see an overall picture. Plotting rapidly becomes difficult as you increase the number of variables to display.

**4. Different levels of interaction**

By contrast, modeling doesn't get much harder as you increase the number of explanatory variables. Here, there are three explanatory variables and the only change is that you include an extra plus.The main tricky thing about including more explanatory variables in models is that there are more options regarding interactions. This model specifies no interactions between variables.You could also include 2-way, or "pairwise" interactions between the explanatory variables.A third option is to include a 3-way interaction between all the explanatory variables.This syntax quickly becomes cumbersome to type, so fortunately there are shortcuts.

**5. All the interactions**

You've already seen the concise syntax for including all the interactions. Simply swap the plus operators for times symbols. Both these formulas mean the same thing.You still need a plus before the zero to denote not including a global intercept term.

**6. Only 2-way interactions**

To get only 2-way interactions in the model, but not the 3-way interaction, you can use a new syntax, namely wrapping the explanatory variables in parentheses and raising them to the power of two.Here, the power of operator has a special meaning; it doesn't square the explanatory variable values. To do that, you need to wrap the terms in the I function, like you saw in the previous course.

**7. The prediction flow**

The prediction flow with an extra variable contains no surprises. This is exactly what you've seen before. Modeling code scales nicely with more variables.

**8. Visualizing predictions**

Likewise, the plotting code is the same as before. The colors of the square prediction points are similar to the colors of the nearby circular data points, indicating that the model provides a good fit.

**9. Let's practice!**

Let's make some models!

## Visualizing many variables

<!-- 
LO: Can visualize more than 3 variables at once.
-->
As you begin to consider more variables, plotting them all at the same time becomes increasingly difficult. In addition to using x and y scales for two numeric variables, you can use color for a third numeric variable, and you can use faceting for categorical variables. And that's about your limit before the plots become to difficult to interpret. There are some specialist plot types like correlation heatmaps and parallel coordinates plots that will handle more variables, but they give you much less information about each variable, and they aren't great for visualizing model predictions.

Here you'll push the limits of the scatter plot by showing the house price, the distance to the MRT station, the number of nearby convenience stores, and the house age, all together in one plot.

**Steps**

1. Using the `taiwan_real_estate` dataset, draw a scatter plot of `n_convenience` versus the square root of `dist_to_mrt_m`, colored by `price_twd_msq`.
2. Use the continuous viridis plasma color scale. 
3. Facet the plot, wrapping by `house_age_years`.

```{r}
# Using taiwan_real_estate, no. of conv. stores vs. sqrt of
# dist. to MRT, colored by plot house price
ggplot(
  taiwan_real_estate, 
  aes(sqrt(dist_to_mrt_m), n_convenience, color = price_twd_msq)
) +
  # Make it a scatter plot
  geom_point() +
  # Use the continuous viridis plasma color scale
  scale_color_viridis_c(option = "plasma") +
  # Facet, wrapped by house age
  facet_wrap(vars(house_age_years))
```

Sublime scatter plotting! As you add more variables, it becomes increasingly challenging to create an easily interpretable plot.

## Different levels of interaction

<!--
LO: Can specify models with three explanatory variables, including interactions.
-->
Once you have three explanatory variables, the number of options for specifying interactions increases. You can specify no interactions. You can specify 2-way interactions, which gives you model coefficients for each pair of variables. The third option is to specify all the interactions, which means the three 2-way interactions and and interaction between all three explanatory variables.

As the number of explanatory variables increases further, the number of interaction possibilities rapidly increases.

**Steps**

1. Fit a linear regression of house price versus `n_convenience`, the square-root of `dist_to_mrt_m`, and `house_age_years`. Don't include a global intercept, and don't include any interactions.

```{r}
# Model price vs. no. of conv. stores, sqrt dist. to MRT 
# station & house age, no global intercept, no interactions
mdl_price_vs_all_no_inter <- lm(
  price_twd_msq ~ n_convenience + sqrt(dist_to_mrt_m) + house_age_years + 0, 
  data = taiwan_real_estate
)

# See the result
mdl_price_vs_all_no_inter
```

2. Fit a linear regression of house price versus the square-root of `dist_to_mrt_m`, `n_convenience`, and `house_age_years`. Don't include a global intercept, but do include 2-way and 3-way interactions between the explanatory variables.

```{r}
# Model price vs. sqrt dist. to MRT station, no. of conv.
# stores & house age, no global intercept, 3-way interactions
mdl_price_vs_all_3_way_inter <- lm(
  price_twd_msq ~ sqrt(dist_to_mrt_m) * n_convenience * house_age_years + 0, 
  data = taiwan_real_estate
)

# See the result
mdl_price_vs_all_3_way_inter
```

3. Fit a linear regression of house price versus the square-root of `dist_to_mrt_m`, `n_convenience`, and `house_age_years`. Don't include a global intercept, but do include 2-way (not 3-way) interactions between the explanatory variables.

```{r}
# Model price vs. sqrt dist. to MRT station, no. of conv.
# stores & house age, no global intercept, 2-way interactions
mdl_price_vs_all_2_way_inter <- lm(
  price_twd_msq ~ (sqrt(dist_to_mrt_m) + n_convenience + house_age_years) ^ 2 + 0, 
  data = taiwan_real_estate
)

# See the result
mdl_price_vs_all_2_way_inter
```

Magnificent model fitting! The formula syntax is flexible enough to provide precise control over which interactions are specified.

## Predicting again

<!--
LO: Run a prediction workflow using multiple explanatory variables.
-->
You've followed the prediction workflow several times now with different combinations of explanatory variables. Time to try it once more on the model with three explanatory variables. Here, you'll use the model with 3-way interactions, though the code is the same when using any of the three models from the previous exercise.

`taiwan_real_estate` and `mdl_price_vs_all_3_way_inter` are available; `dplyr`, `tidyr` and `ggplot2` are loaded.

**Steps**

1. Make a grid of explanatory data, formed from combinations of the following variables.

    * `dist_to_mrt_m` should take a sequence from zero to eighty in steps of ten, all squared (0, 100, 400, …, 6400).
    * `n_convenience` should take the numbers zero to ten.
    * `house_age_years` should take the unique values of the `house_age_years` column of `taiwan_real_estate`.

```{r}
# Make a grid of explanatory data
explanatory_data <- expand_grid(
  # Set dist_to_mrt_m a seq from 0 to 80 by 10s, squared
  dist_to_mrt_m = seq(0, 80, 10) ^ 2,
  # Set n_convenience to 0 to 10
  n_convenience = 0:10,
  # Set house_age_years to the unique values of that variable
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# See the result
explanatory_data
```

2. Add a column to the `explanatory_data`, assigning to `prediction_data`.
3. The column should be named after the response variable, and contain predictions made using `mdl_price_vs_all_3_way_inter` and `explanatory_data`.

```{r}
# From previous step
explanatory_data <- expand_grid(
  dist_to_mrt_m = seq(0, 80, 10) ^ 2,
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)

# Add predictions to the data frame
prediction_data <- explanatory_data %>% 
  mutate(
    price_twd_msq = predict(mdl_price_vs_all_3_way_inter, explanatory_data)
  )

# See the result
prediction_data
```

4. Extend the plot to include predictions as points from `prediction_data`, with size 3 and shape 15.
5. *Look at the plot. What do the prediction points tell you?*

```{r}
# From previous step
explanatory_data <- expand_grid(
  dist_to_mrt_m = seq(0, 80, 10) ^ 2,
  n_convenience = 0:10,
  house_age_years = unique(taiwan_real_estate$house_age_years)
)
prediction_data <- explanatory_data %>% 
  mutate(price_twd_msq = predict(mdl_price_vs_all_3_way_inter, explanatory_data))

# Extend the plot
ggplot(
  taiwan_real_estate, 
  aes(sqrt(dist_to_mrt_m), n_convenience, color = price_twd_msq)
) +
  geom_point() +
  scale_color_viridis_c(option = "plasma") +
  facet_wrap(vars(house_age_years)) +
  # Add points from prediction data, size 3, shape 15
  geom_point(data = prediction_data, size = 3, shape = 15)
```

Premier predicting! The plot nicely shows that the house price decreases as the square-root of the distance to the nearest MRT station increases, and increases as the number of nearby convenience stores increases, and is higher for houses under 15 years old.

## How linear regression works

Theory. Coming soon ...


**1. How linear regression works**

Let's see how linear regression works. To keep things understandable, we'll stick to simple linear regression, with a single numeric explanatory variable.

**2. The standard simple linear regression plot**

Here is the standard scatter plot of a simple linear regression.

**3. Visualizing residuals**

Here's the same plot, showing the residuals. That is, the actual response minus the predicted response.For the best fit, we want those lines to be as short as possible. That is, we want a metric that measures the size of all the residuals, and we want to make that as small as possible.

**4. A metric for the best fit**

The simplest idea for that metric would be to add up all the residuals.This doesn't work because some residuals are negative, so they would make the total smaller instead of larger.Instead, we do the next easiest thing, which is to square each residual so that they are non-negative, and then add them up.This metric is called the sum of squares.The tricky part is determining which intercept and slope coefficients will result in the smallest sum of squares.

**5. A detour into numerical optimization**

To solve this problem, we need to take a detour into numerical optimization, which means finding the minimum point of a function.Consider this quadratic equation.  The y-value is x-squared minus x plus ten. The plot shows that the minimum point of the function occurs when x is a little above zero and y is a little below ten, but how can we find it exactly?

**6. Using calculus to solve the equation**

I can solve this with calculus by taking the derivative, setting that derivative to zero, rearranging for x, then substituting back into the original equation to find y.It gets the right answer. x is zero-point-five and y is nine-point-seven-five.However, not all equations can be solved in this analytic fashion, and ever since Newton and Liebniz invented calculus, mathematicians have been trying to find ways of avoiding it. In fact, one of the perks of being a data scientist is that you can just let R figure out how to find the minimum.

**7. optim()**

To perform numerical optimization in R, you need a function to minimize. In this case, it takes x as an input, and returns x-squared minus x plus ten.Then you call optim. The first argument is an initial guess at the answer. For more complicated functions, this is sometimes important, but here you could pick anything. The second argument is the function to call, without parentheses.In the output, par is the estimate for the x-value. It's close to the correct answer of zero-point-five, and if you need better accuracy, there are many options you can play with to improve the answer.value is the estimated y-value, which is spot on.The other pieces of output are diagnostic values, which we don't need here.

**8. Slight refinements**

There are two changes to the code which aren't needed here but will assist us in the linear regression case.The function passed to optim is only allowed to have one argument, so to optimize for multiple variables, you need to pass them as a numeric vector. Here, calc_quadratic takes a numeric vector called coeffs, then extracts the first element and calls it x.Secondly, passing a named vector as the initial guess to the par argument of optim makes the output easier to read. Now the par element in the output is named "x".

**9. A linear regression algorithm**

While lm() is hundreds of lines of code, you can implement simple linear regression for a specific dataset in just a few lines.You define a function that accepts the intercept and slope, and returns the sum of the squares of residuals. You'll have to use the trick of giving the function a single coeffs argument, then extracting the individual numbers.You'll perform the rest of the calculation yourself in the exercises.Then you call optim, passing an initial guess for the coefficients and your sum of squares function.That's it!

**10. Let's practice!**

Time to delve into linear regression's internals!

## The sum of squares

<!--
LO: Understands that linear regression involves minimizing the sum of squares metric.
-->
In order to choose the "best" line to fit the data, regression models need to optimize some metric. For linear regression, this metric is called the *sum of squares*.

In the dashboard, try setting different values of the intercept and slope coefficients. In the plot, the solid black line has the intercept and slope you specified. The dotted blue line has the intercept and slope calculated by a linear regression on the dataset.

> *Question*
> ---
> How does linear regression try to optimize the sum of squares metric?<br>
> <br>
> ⬜ Linear regression makes the sum of the squares of the differences between the actual responses and the predicted responses zero.<br>
> ⬜ Linear regression makes the sum of the squares of the differences between the actual responses and the predicted responses infinite.<br>
> ⬜ Linear regression maximizes the sum of the squares of the differences between the actual responses and the predicted responses.<br>
> ✅ Linear regression minimizes the sum of the squares of the differences between the actual responses and the predicted responses.<br>

Magic minimization! Sum of squares is a measure of how far the predited responses are from the actual responses, so a smaller number is better. 

## Linear regression algorithm

<!--
LO: Can write a simple linear regression algorithm for a specific dataset. 
-->
To truly understand linear regression, it is helpful to know how the algorithm works. The code for `lm()` is hundreds of lines because it has to work with any formula and any dataset. However, in the case of simple linear regression for a single dataset, you can implement a linear regression algorithm in just a few lines of code.

The workflow is

1. Write a script to calculate the sum of squares.
2. Turn this into a function.
3. Use R's general purpose optimization function find the coefficients that minimize this. 

The explanatory values (the `n_convenience` column of `taiwan_real_estate`) are available as `x_actual`.
The response values (the `price_twd_msq` column of `taiwan_real_estate`) are available as `y_actual`.

**Steps**

1. Set the intercept to ten.
2. Set the slope to one.
3. Calculate the predicted y-values as the intercept plus the slope times the actual x-values.
4. Calculate the differences between actual and predicted y-values.
5. Calculate the sum of squares. Get the sum of the differences in y-values, squaring each value.

```{r}
# set data
x_actual <- taiwan_real_estate$n_convenience
y_actual <- taiwan_real_estate$price_twd_msq

# Set the intercept to 10
intercept <- 10

# Set the slope to 1
slope <- 1

# Calculate the predicted y values
y_pred <- intercept + slope * x_actual

# Calculate the differences between actual and predicted
y_diff <- y_actual - y_pred

# Calculate the sum of squares
sum(y_diff ^ 2)
```

6. Complete the function body.

    * Get the intercept from the first element of `coeffs`.
    * Get the slope from the second element of `coeffs`.
    * Calculate the predicted y-values as the intercept plus the slope times the actual x-values.
    * Calculate the differences between actual and predicted y-values.
    * Calculate the sum of squares. Get the sum of the differences in y-values, squaring each value.

```{r}
calc_sum_of_squares <- function(coeffs) {
  # Get the intercept coeff
  intercept <- coeffs[1]

  # Get the slope coeff
  slope <- coeffs[2]

  # Calculate the predicted y values
  y_pred <- intercept + slope * x_actual

  # Calculate the differences between actual and predicted
  y_diff <- y_actual - y_pred

  # Calculate the sum of squares
  sum(y_diff ^ 2)
}
```

7. Optimize the sum of squares metric.

    * Call an optimization function.
    * Initially guess that the intercept is zero and the slope is zero by passing a named vector of parameters.
    * Use `calc_sum_of_squares` as the optimization function.

```{r}
# From previous step
calc_sum_of_squares <- function(coeffs) {
  intercept <- coeffs[1]
  slope <- coeffs[2]
  y_pred <- intercept + slope * x_actual
  y_diff <- y_actual - y_pred
  sum(y_diff ^ 2)
}

# Optimize the metric
optim(
  # Initially guess 0 intercept and 0 slope
  par = c(intercept = 0, slope = 0), 
  # Use calc_sum_of_squares as the optimization fn 
  fn = calc_sum_of_squares
)

# Compare the coefficients to those calculated by lm()
lm(price_twd_msq ~ n_convenience, data = taiwan_real_estate)
```

Outstanding optimization! The results you got here with just a few lines of code are very close to the finely-tuned results from `lm()`. All you needed was a function to calculate the sum of squares metric, and the `optim()` function worked its magic to find where this function had its minimum value.

# 4. Multiple Logistic Regression

Extend your logistic regression skills to multiple explanatory variables. Understand the logistic distribution, which underpins this form of regression. Finally, implement your own logistic regression algorithm.

## Multiple logistic regression

Theory. Coming soon ...

**1. Multiple logistic regression**

Let's switch from linear regression to logistic regression.

**2. Bank churn dataset**

We'll revisit the three-column bank churn dataset from the previous course. has_churned is the response, denoting whether or not the customer churned, time_since_first_purchase is a measure of the length of relationship with the customer, and time_since_last_purchase is a measure of the recency of activity of the customer. The explanatory variables have been transformed to protect commercially sensitive information.

    1 https://www.rdocumentation.org/packages/bayesQR/topics/Churn

**3. glm()**

Recall that to perform a logistic regression, there are two changes compared to a linear regression. Firstly, you call glm, for generalized linear models, rather than lm. Secondly, you include a family argument to specify the error distribution, set to binomial. You'll explore the binomial function later in the chapter. To extend logistic regression to multiple explanatory variables, you change the formula in the same way as linear regression, with a plus to ignore interactions, or a times to include interactions. There's no new syntax here.

**4. Prediction flow**

The prediction flow should also feel familiar, since you've seen all the techniques already. Use expand_grid to make a tibble of explanatory variables, then mutate to add a column of predictions. The only change from the linear regression case is that you need to specify type equals "response" in the call to predict.

**5. The four outcomes**

Recall that when the response variable has two possible values, there are four outcomes for the model. Either it correctly predicts positive and negative responses, or it gets it wrong with a false positive or false negative. We can quantify and visualize these four outcomes using a confusion matrix.

    1 https://campus.datacamp.com/courses/introduction-to-regression-in-r/simple-logistic-regression?ex=10

**6. Confusion matrix**

The code flow is the same as before. Get the actual responses from the dataset and the predicted responses from the model, rounded to give zeroes and ones. Then use table to get counts of each of the four outcomes. yardstick's conf_mat function converts it to a confusion matrix object. That lets you plot the result as a mosaic plot using autoplot, and display metrics like model accuracy, sensitivity and specificity, using summary.

**7. Visualization**

Visualizing the plot when you have multiple explanatory variables is trickier. As with linear regression visualizations, you can use faceting to provide different panels for categorical variables. For the case of two numeric explanatory variables, you can use the technique of mapping the response variable to color. A nice trick here is to give predicted probabilities less than zero-point-five one color, and predicted probabilities above zero-point-five another color. This is achieved with ggplot2's scale_color_gradient2 function, setting midpoint to zero-point-five. You'll see how it looks in the exercises.

**8. Let's practice!**

Let's get logistic! 


## Visualizing multiple explanatory variables

<!--
LO: Can visualize logistic data with two numeric explanatory variables.
-->
Logistic regression also supports multiple explanatory variables. Plotting has similar issues as the linear regression case: it quickly becomes difficult to include more numeric variables in the plot. Here we'll look at the case of two numeric explanatory variables, and the solution is basically the same as before: use color to denote the response. 

Here there are only two possible values of response (zero and one), and later when we add predicted responses, the values all lie between zero and one. Once you include predicted responses, the most important thing to determine from the plot is whether the predictions are close to zero, or close to one. That means that a 2-color gradient split at 0.5 is really useful: responses above 0.5 are one color, and responses below 0.5 are another color.

**Steps**

1. Using the `churn` dataset, plot the recency of purchase, `time_since_last_purchase`, versus the length of customer relationship, `time_since_first_purchase`, colored by whether or not the customer churned, `has_churned`.
2. Add a point layer, with transparency set to `0.5`.
3. Use a 2-color gradient, with midpoint `0.5`.
4. Use the black and white theme.

```{r}
# Load data
churn <- read_fst("data/churn.fst")

# Using churn, plot recency vs. length of relationship,
# colored by churn status
ggplot(
  churn, 
  aes(time_since_first_purchase, time_since_last_purchase, color = has_churned)
) +
  # Make it a scatter plot, with transparency 0.5
  geom_point(alpha = 0.5) +
  # Use a 2-color gradient split at 0.5
  scale_color_gradient2(midpoint = 0.5) +
  # Use the black and white theme
  theme_bw()
```

Cool colored scatter plot! The 2-color gradient is excellent for distinguishing the two cases of a positive and negative response.

## Logistic regression with 2 explanatory variables

<!--
LO: Can fit a logistic regression with multiple explanatory variables.
-->
To include multiple explanatory variables in logistic regression models, the syntax is the same as for linear regressions. The only change is the same as in the simple case: you run a *generalized* linear model with a binomial error family.

Here you'll fit a model of churn status with both of the explanatory variables from the dataset: the length of customer relationship and the recency of purchase.

**Steps**

1. Fit a logistic regression of churn status, `has_churned` versus length of customer relationship, `time_since_first_purchase` and recency of purchase, `time_since_last_purchase`, and an interaction between the explanatory variables.

```{r}
# Fit a logistic regression of churn status vs. length of
# relationship, recency, and an interaction
mdl_churn_vs_both_inter <- glm(
  has_churned ~ time_since_first_purchase * time_since_last_purchase, 
  data   = churn, 
  family = binomial
)

# See the result
mdl_churn_vs_both_inter
```

Laudable logistic regression! The formula syntax for including multiple explanatory variables is the same for logistic regression as it is for linear regression.

## Logistic regression prediction

<!--
LO: Can make predictions with a logistic regression model that includes multiple explanatory variables.
-->
As with linear regression, the joy of logistic regression is that you can make predictions. Let's step through the prediction flow one more time!

**Steps**

1. Create a grid of explanatory variables.

    * Set `time_since_first_purchase` to a sequence from minus two to four in steps of `0.1`.
    * Set `time_since_last_purchase` to a sequence from minus one to six in steps of `0.1`.

```{r}
# Make a grid of explanatory data
explanatory_data <- expand_grid(
  # Set len. relationship to seq from -2 to 4 in steps of 0.1
  time_since_first_purchase = seq(-2, 4, 0.1),
  # Set recency to seq from -1 to 6 in steps of 0.1
  time_since_last_purchase = seq(-1, 6, 0.1)
)

# See the result
explanatory_data
```

2. Add a column to `explanatory_data` named `has_churned` containing predictions using `mdl_churn_vs_both_inter` and `explanatory_data` with type `"response"`.

```{r}
# From previous steps
explanatory_data <- expand_grid(
  time_since_first_purchase = seq(-2, 4, 0.1),
  time_since_last_purchase = seq(-1, 6, 0.1)
)

# Add a column of predictions using mdl_churn_vs_both_inter
# and explanatory_data with type response
prediction_data <- explanatory_data %>% 
  mutate(
    has_churned = predict(
      mdl_churn_vs_both_inter, explanatory_data, type = "response"
    )
  )

# See the result
prediction_data
```

3. Extend the plot by adding points from `prediction_data` with size 3 and shape 15.

```{r}
# From previous steps
explanatory_data <- expand_grid(
  time_since_first_purchase = seq(-2, 4, 0.1),
  time_since_last_purchase = seq(-1, 6, 0.1)
)
prediction_data <- explanatory_data %>% 
  mutate(
    has_churned = predict(mdl_churn_vs_both_inter, explanatory_data, type = "response")
  )

# Extend the plot
ggplot(
  churn, 
  aes(time_since_first_purchase, time_since_last_purchase, color = has_churned)
) +
  geom_point(alpha = 0.5) +
  scale_color_gradient2(midpoint = 0.5) +
  theme_bw() +
  # Add points from prediction_data with size 3 and shape 15
  geom_point(data = prediction_data, size = 3, shape = 15)
```

Painless probability predictions! The prediction workflow should be familiar now: you've seen it many times before. The only thing to remember here is to set the prediction type to `"response"`.

## Confusion matrix

<!--
LO: Can generate a confusion matrix from a dataset and logistic regression model, and plot in a mosaic plot.
-->
When the response variable has just two outcomes, like the case of churn, the measures of success for the model are "how many cases where the customer churned did the model correctly predict?" and "how many cases where the customer didn't churn did the model correctly predict?". These can be found by generating a confusion matrix and calculating summary metrics on it. A mosaic plot is the natural way to visualize the results.

**Steps**

1. Get the actual responses from the `churn` dataset.
2. Get the predicted responses from the rounded, fitted values of `mdl_churn_vs_both_inter`.
3. Create a table of the actual and predicted response values.
4. Convert the table to a `conf_mat` confusion matrix object.

```{r}
# Load packages
library(yardstick)

# Get the actual responses from churn
actual_response <- churn$has_churned

# Get the predicted responses from the model
predicted_response <- round(fitted(mdl_churn_vs_both_inter))

# Get a table of these values
outcomes <- table(predicted_response, actual_response)

# Convert the table to a conf_mat object
confusion <- conf_mat(outcomes)

# See the result
confusion
```

5. "Automatically" plot the confusion matrix, `confusion`.
6. Get summary metrics from the confusion matrix. Remember that the churn event is in the second row/column of the matrix.

```{r}
# From previous step
actual_response <- churn$has_churned
predicted_response <- round(fitted(mdl_churn_vs_both_inter))
outcomes <- table(predicted_response, actual_response)
confusion <- conf_mat(outcomes)

# "Automatically" plot the confusion matrix
autoplot(confusion)

# Get summary metrics
summary(confusion, event_level = "second")
```

Metric magic! Generating a confusion matrix and calculating metrics like accuracy, sensitivity, and specificity is the standard way to measure how well a logistic model fits.

## The logistic distribution

Theory. Coming soon ...


**1. The logistic distribution**

In order to understand logistic regression, you need to know about the logistic distribution.

**2. Gaussian probability density function (PDF)**

Before we get to the logistic distribution, let's look at the Gaussian, or normal distribution. Hopefully, you are familiar with the famous "bell curve" of its probability density function, made with the dnorm function.For the purposes of regression, we care more about the area under this curve. By integrating the dnorm function - calculating the area underneath it - we get another curve, known as the cumulative distribution function.

**3. Gaussian cumulative distribution function (CDF)**

To get the cumulative distribution function, or CDF, you call pnorm instead of dnorm.The y-axis is near zero on the far left of the plot, and near one on the far right of the plot. This is a feature of the CDF curve for all distributions. When x has its minimum possible value, in this case minus infinity, y will be zero. When x has its maximum possible value, in this case infinity, y will be one.You can think of the CDF as a transformation from the values of x to probabilities. When x is one, the CDF is curve is at zero-point-eight-four. That means that for a normally distributed variable x, the probability that x is less than one is eighty-four percent.

**4. Gaussian inverse CDF**

Since the CDF transforms from x-values to probabilities, you also need a way to get back from probabilities to x-values. This is the inverse CDF.Here we have a new dataset with probabilities from nearly zero to nearly one. The inverse CDF is calculated with qnorm.The line plot you see is the same as the CDF plot from the previous slide, but with the x and y axes flipped.

**5. Distribution function names**

The function names for distribution curves all follow a pattern. The PDF function starts with "d", the CDF function starts with "p", and the inverse CDF starts with "q". Then the names end with an abbreviation of the distribution name.dlogis, plogis, and qlogis follow the same naming convention as the normal distribution functions.

**6. glm()'s family argument**

Performing linear regression with lm is the same as performing it with glm and setting the error distribution family to gaussian.Switching from linear regression to logistic regression is done by changing the family argument to binomial.So, what are these family arguments?

**7. gaussian()**

gaussian is a function. Calling it and wrapping the result in str shows the structure - it returns an object that contains several other functions.Between them, these functions contain all the details for turning a generalized regression into a specific type of regression like linear or logistic regression.

**8. linkfun and linkinv**

Two elements of the family object are especially important. The linkfun element provides a transformation of the response variables, and the linkinv element undoes that transformation.For the gaussian case, it is boring because each transformation is just the identity function: it returns the same value that you put into it, since no special transformation is needed. For logistic regression however, it becomes more exciting.

**9. Logistic PDF**

Here's the logistic distribution PDF. It looks a little bit like the Gaussian PDF, but the tails at the extreme left and right of the plot are fatter.

**10. Logistic distribution**

The CDF for the logistic distribution is also known as the logistic function. The two terms are interchangeable.It has a fairly simple equation: one divided by one plus e to the minus x.The inverse CDF is sometimes called the logit function; again, the terms are interchangeable.Its equation is the logarithm of one divided by one plus p.In order to see what these curves look like, you'll have to try the exercises.

**11. Let's practice!**

Let's look at those curves!

## Cumulative distribution function

<!--
LO: Knows what a logistic distribution CDF looks like.
-->
Understanding the logistic distribution is key to understanding logistic regression. Like the normal (Gaussian) distribution, it is a probability distribution of a single continuous variable. Here you'll visualize the *cumulative distribution function* (CDF) for the logistic distribution. That is, if you have a logistically distributed variable, `x`, and a possible value, `xval`, that `x` could take, then the CDF gives the probability that `x` is less than `xval`.

The logistic distribution's CDF is calculated with the logistic function (hence the name). The plot of this has an S-shape, known as a *sigmoid curve*. An important property of this function is that it takes an input that can be any number from minus infinity to infinity, and returns a value between zero and one.

**Steps**

1. Create a tibble containing three columns.

    * `x` values as a sequence from minus ten to ten in steps of `0.1`.
    * `logistic_x` made from `x` transformed with the logistic distribution CDF.
    * logistic_x_man` made from `x` transformed with a logistic function calculated from the equation <img src="https://render.githubusercontent.com/render/math?math={\color{red}cdf(x)=\frac{1}{(1%2Bexp(-x))}}">.
    * Check that both logistic transformations (`logistic_x` and `logistic_x_man`) have the same values with `all.equal()`.

```{r}
logistic_distn_cdf <- tibble(
  # Make a seq from -10 to 10 in steps of 0.1
  x = seq(-10, 10, 0.1),
  # Transform x with built-in logistic CDF
  logistic_x = plogis(x),
  # Transform x with manual logistic
  logistic_x_man = 1 / (1 + exp(-x))
) 

# Check that each logistic function gives the same results
all.equal(
  logistic_distn_cdf$logistic_x, 
  logistic_distn_cdf$logistic_x_man
)
```

2. Using the `logistic_distn_cdf` dataset, plot `logistic_x` versus `x` as a line plot.

```{r}
# From previous step
logistic_distn_cdf <- tibble(
  x = seq(-10, 10, 0.1),
  logistic_x = plogis(x),
  logistic_x_man = 1 / (1 + exp(-x))
)

# Using logistic_distn_cdf, plot logistic_x vs. x
ggplot(logistic_distn_cdf, aes(x, logistic_x)) +
  # Make it a line plot
  geom_line()
```

Super sigmoid plotting! The logistic distribution's cumulative distribution function has an “S” shape, known as a sigmoid curve.

## Inverse cumulative distribution function

<!--
LO: Knows what a logistic distribution inverse CDF looks like.
-->
The logistic function (logistic distribution CDF) has another important property: each x input value is transformed to a unique value. That means that the transformation can be reversed. The *logit function* is the name for the *inverse logistic function*, which is also the *logistic distribution inverse cumulative distribution function*. (All three terms mean exactly the same thing.)

The logit function takes values between zero and one, and returns values between minus infinity and infinity.

**Steps**

1. Create a tibble containing three columns.

    * `x` values as a sequence from minus `0.001` to `0.999` in steps of `0.001`.
    * `logit_p` made from `p` transformed with the logistic distribution inverse CDF.
    * `logit_p_man` made from `p` transformed with the equation <img src="https://render.githubusercontent.com/render/math?math={\color{red}log(\frac{p}{1 - p})}">.
    * Check that both logit transformations (`logit_p` and `logit_p_man`) have the same values with `all.equal()`.

```{r}
logistic_distn_inv_cdf <- tibble(
  # Make a seq from 0.001 to 0.999 in steps of 0.001
  p = seq(0.001, 0.999, 0.001),
  # Transform with built-in logistic inverse CDF
  logit_p = qlogis(p),
  # Transform with manual logit
  logit_p_man = log(p / (1 - p))
) 

# Check that each logistic function gives the same results
all.equal(
  logistic_distn_inv_cdf$logit_p,
  logistic_distn_inv_cdf$logit_p_man
)
```

2. Using the `logistic_distn_inv_cdf` dataset, plot `logit_p` versus `p` as a line plot.

```{r}
# From previous step
logistic_distn_inv_cdf <- tibble(
  p = seq(0.001, 0.999, 0.001),
  logit_p = qlogis(p),
  logit_p_man = log(p / (1 - p))
)

# Using logistic_distn_inv_cdf, plot logit_p vs. p
ggplot(logistic_distn_inv_cdf, aes(p, logit_p)) +
  # Make it a line plot
  geom_line()
```

Incredible inversing! The inverse CDF is the "opposite" transformation to the CDF. If you flip the x and y axes on this plot, you get the same plot you saw in the previous exercise.

## binomial family argument

<!--
LO: Understands the fmaily argument for glm(), in the context of logistic regression.
-->
The big difference between running a linear regression with `lm()` and running a logistic regression with `glm()` is that you have to set `glm()`'s `family` argument to `binomial`. `binomial()` is a function that returns a list of other functions that tell `glm()` how to perform calculations in the regression. The two most interesting functions are `linkinv` and `linkfun`, which are used for transforming variables from the whole number line (minus infinity to infinity) to probabilities (zero to one) and back again.

**Steps**

1. Examine the structure of the `binomial()` function. *Notice that it contains two elements that are functions,* `binomial()$linkinv`*, and* `binomial()$linkfun`.
2. Call `binomial()$linkinv()` on `x`, assigning to `linkinv_x`.
3. Check that `linkinv_x` and `plogis()` of `x` give the same results with `all.equal()`.
4. Call `binomial()$linkfun()` on `p`, assigning to `linkfun_p`.
5. Check that `linkfun_p` and `qlogis()` of `p` give the same results.

```{r}
# Create data
x <- seq(-10, 10, 0.2)
p <- seq(0.01, 0.99, 0.01)

# Look at the structure of binomial() function
str(binomial())

# Call the link inverse on x
linkinv_x <- binomial()$linkinv(x)

# Check linkinv_x and plogis() of x give same results 
all.equal(linkinv_x, plogis(x))

# Call the link fun on p
linkfun_p <- binomial()$linkfun(p)

# Check linkfun_p and qlogis() of p give same results 
all.equal(linkfun_p, qlogis(p))
```

Terrific transforming! The  `binomial()` family object contains `linkinv` and `linkfun` functions that correspond to the logistic distribution CDF and inverse CDF respectively. These are used to translate between numbers and probabilities.

## Logistic distribution parameters

<!--
LO:
-->
The logistic distribution CDF is not just a single curve. In the same way that the normal distribution has mean and standard deviation parameters that affect the CDF curve, the logistic distribution has *location* and *scale* parameters. Here, you'll visualize how changing those parameters changes the CDF curve.

> *Question*
> ---
> How do changes to the parameters change the CDF curve?<br>
> <br>
> ⬜ As `location` increases, the logistic CDF curve moves rightwards. As `scale` increases, the steepness of the slope increases.<br>
> ⬜ As `location` increases, the logistic CDF curve moves upwards. As `scale` increases, the steepness of the slope increases.<br>
> ✅ As `location` increases, the logistic CDF curve moves rightwards. As `scale` increases, the steepness of the slope decreases.<br>
> ⬜ As `location` increases, the logistic CDF curve moves upwards. As `scale` increases, the steepness of the slope decreases.<br>

Perfect parameter play! The logistic distribution consists of a whole family of curves specified by the location and scale parameters. This allows logistic model prediction curves to have different positions or steepnesses. 

## How logistic regression works

Theory. Coming soon ...


**1. How logistic regression works**

Let's see how logistic regression works. The principle is the same as for linear regression: choose a metric that measures how far the predicted responses are from the actual responses, and optimize that metric.

**2. Sum of squares doesn't work**

In the linear regression case, the metric to optimize was the sum of squares. That is, you calculated each predicted response minus the corresponding actual response, squared it, then took the sum.In the case of logistic regression, the actual response is always either zero or one, and the predicted response is between these two values.It turns out that the sum of squares metric optimizes poorly under these restrictions, and that there is a better metric.

**3. Likelihood**

This is the likelihood metric. Unlike the sum of squares, where the goal was to find the minimum possible value, with likelihood you want to find the maximum value.You take the product of the predicted and actual responses,

**4. Likelihood**

and add the product of one minus the predicted responses and one minus the actual responses.

**5. Likelihood**

Then you sum over all data points.Since the actual response only has two possible values, this equation simplifies in two different ways.When the actual response is one, the equation for each observation simplifies to the predicted response, y_pred. As y_pred increases, the metric increases too, and the maximum likelihood occurs when y_pred is one, the same as the actual value.When the actual response is zero, the equation simplifies to one minus the predicted response. As y_pred decreases, the metric increases, and the maximum likelihood occurs when y_pred is zero.In either case, you get a higher likelihood score when the predicted response is close to the actual response.

**6. Log-likelihood**

When calculating the likelihood, y_pred is often close to zero or one, which means you end up adding up lots of very small numbers, which introduces numerical error.It is more efficient to compute the log-likelihood. The only difference in this equation is that you take the logarithm of the predicted response terms.Optimizing to find the log-likelihood gives the same coefficients as optimizing to find the likelihood.

**7. Negative log-likelihood**

Since we want to maximize likelihood, but the optim function defaults to finding minimums, one tweak to make is to calculate the negative log-likelihood. That is, include a minus sign when you calculate the sum of each observation's likelihood contribution.

**8. Logistic regression algorithm**

Now we are set to write our logistic regression algorithm. The metric function takes a coefficients argument, you extract the intercept and slope from it, perform some more calculation that you'll see in the exercises, and then call optim to find the coefficients that minimize the metric.

**9. Let's practice!**

Almost there!

## Likelihood & log-likelihood

<!--
LO: Understands that logistic regression involves maximizing the likelihood metric. 
-->
Linear regression tries to optimize a "sum of squares" metric in order to find the best fit. That metric isn't applicable to logistic regression. Instead, logistic regression tries to optimize a metric called *likelihood*, or a related metric called *log-likelihood*.

The dashboard shows churn status versus line time since last purchase from the `churn` dataset. The blue dotted line is the logistic regression prediction line calculated by ggplot2's `geom_smooth()`. (That is, it's the "best fit" line.) The black solid line shows a prediction line calculated from the intercept and slope coefficients you specify as `plogis(intercept + slope * time_since_last_purchase)`.

Change the intercept and slope coefficients and watch how the likelihood and log-likelihood values change.

> *Question*
> ---
> As you get closer to the best fit line, what statement is true about likelihood and log-likelihood?<br>
> <br>
> ✅ Both likelihood and log-likelihood increase to a maximum value.<br>
> ⬜ Both likelihood and log-likelihood decrease to a minimum value.<br>
> ⬜ Likelihood increases to a maximum value; log-likelihood decreases to a minimum value.<br>
> ⬜ Likelihood decreases to a minimum value; log-likelihood increases to a maximum value.<br>

It seems likely that you know what you are doing! Logistic regression chooses the prediction line that gives you the maximum likelihood value. It also gives maximum log-likelihood. 

## Logistic regression algorithm

<!--
LO: Can write a simple logistic regression algorithm for a specific dataset. 
-->
Let's dig into the internals and implement a logistic regression algorithm. Since R's `glm()` function is very complex, you'll stick to implementing simple logistic regression for a single dataset.

Rather than using sum of squares as the metric, we want to use likelihood. However, log-likelihood is more computationally stable, so we'll use that instead. Actually, there is one more change: since we want to maximize log-likelihood, but `optim()` defaults to finding minimum values, it is easier to calculate the *negative* log-likelihood.

The log-likelihood value for each observation is

<img src="https://render.githubusercontent.com/render/math?math={\color{red}log(y_{pred}) * y_{actual} %2B log(1 - y_{pred}) * (1 - y_{actual})}">

The metric to calculate is minus the sum of these log-likelihood contributions.

The explanatory values (the `time_since_last_purchase` column of `churn`) are available as `x_actual`.
The response values (the `has_churned` column of `churn`) are available as `y_actual`.

**Steps**

1. Set the intercept to one.
2. Set the slope to `0.5`.
3. Calculate the predicted y-values as the intercept plus the slope times the actual x-values, all transformed with the logistic distribution CDF.
4. Calculate the log-likelihood for each term as the log of the predicted y-values times the actual y-values, plus the log of one minus the predicted y-values times one minus the actual y-values.
5. Calculate minus the sum of the log-likelihoods for each term.

```{r}
# Set the intercept to 1
intercept <- 1

# Set the slope to 0.5
slope <- 0.5

# Calculate the predicted y values
y_pred <- plogis(intercept + slope * x_actual)

# Calculate the log-likelihood for each term
log_likelihoods <- log(y_pred) * y_actual + log(1 - y_pred) * (1 - y_actual)

# Calculate minus the sum of the log-likelihoods for each term
-sum(log_likelihoods)
```

6. Complete the function body.

    * Get the intercept from the first element of `coeffs`.
    * Get the slope from the second element of `coeffs`.
    * Calculate the predicted y-values as the intercept plus the slope times the actual x-values, transformed with the logistic distribution CDF.
    * Calculate the log-likelihood for each term as the log of the predicted y-values times the actual y-values, plus the log of one minus the predicted y-values times one minus the actual y-values.
    * Calculate minus the sum of the log-likelihoods for each term.

```{r}
calc_neg_log_likelihood <- function(coeffs) {
  # Get the intercept coeff
  intercept <- coeffs[1]

  # Get the slope coeff
  slope <- coeffs[2]

  # Calculate the predicted y values
  y_pred <- plogis(intercept + slope * x_actual)

  # Calculate the log-likelihood for each term
  log_likelihoods <- log(y_pred) * y_actual + log(1 - y_pred) * (1 - y_actual)

  # Calculate minus the sum of the log-likelihoods for each term
  -sum(log_likelihoods)
}
```

7. Optimize the sum of squares metric.

    * Call an optimization function.
    * Initially guess that the intercept is zero and the slope is one.
    * Use `calc_neg_log_likelihood` as the optimization function.

```{r}
# From previous step
calc_neg_log_likelihood <- function(coeffs) {
  intercept <- coeffs[1]
  slope <- coeffs[2]
  y_pred <- plogis(intercept + slope * x_actual)
  log_likelihoods <- log(y_pred) * y_actual + log(1 - y_pred) * (1 - y_actual)
  -sum(log_likelihoods)
}

# Optimize the metric
optim(
  # Initially guess 0 intercept and 1 slope
  par = c(intercept = 0, slope = 1),
  # Use calc_neg_log_likelihood as the optimization fn 
  fn = calc_neg_log_likelihood
)

# Compare the coefficients to those calculated by glm()
glm(has_churned ~ time_since_last_purchase, data = churn, family = binomial)
```

Ace algorithmic implementation! To make a really simple version of `glm()`, you just needed a function to calculate the negative log-likelihood, and a general-purpose optimization function.

## Congratulations

Theory. Coming soon ...


**1. Congratulations**

Well done! You've learned linear and logistic regression!

**2. You learned things**

In Chapter 1 you fitted, visualized, predicted and assessed parallel slopes linear regression models.In Chapter 2 you explored interactions between explanatory variables, and tried to resolve Simpson's Paradox.In Chapter 3 you saw that although visualization gets tricky with more explanatory variables, modeling easily handles them.In Chapter 4 you ran logistic regression with multiple explanatory variables, and explored the logistic distribution.On top of that, you implemented algorithms for linear regression and logistic regression. That's a really advanced thing, so well done!

**3. There is more to learn**

In order to master regression, there are a few more things to learn that we didn't have time for in this course.Firstly, it is common practice to split your data into separate training and testing datasets. You calculate the model coefficients on the training set, then assess the performance of the model on the testing set. This helps prevent overfitting, where your model only works well for one specific dataset.Secondly, an even stronger remedy against overfitting is to use cross-validation. This involves fitting and assessing the model several times, against several random training-testing splits to give a set of models from which to choose.Thirdly, in order to help determine which explanatory variables should be included in the model, linear and logistic regression return the significance of each coefficient, which is a measure of whether the coefficient is really different from zero, or it just appears to be different from zero by chance.

**4. Advanced regression**

DataCamp has many more courses on advanced regression techniques. You're now ready to take them!

**5. Let's practice!**

Have fun regressing!

