```{r global-options, include=FALSE}
library(tidyverse)
library(tidymodels)
library(palmerpenguins)
library(patchwork)
library(marginaleffects)
library(kableExtra)

tidymodels_prefer()

penguins <- palmerpenguins::penguins

options(kable_styling_bootstrap_options = c("hover", "striped"))

theme_set(theme_bw(base_size = 14))

set.seed(123)
penguins_split <- initial_split(penguins)
penguins_train <- training(penguins_split)
penguins_test <- testing(penguins_split)
```

## Where We Stand

We've come a long way over the past few weeks. We've followed the trajectory below.

+ **Simple Linear Regression Models:** $\mathbb{E}\left[y\right] = \beta_0 + \beta_1 \cdot x_1$

  + *General Assumption*: We can predict the response, $y$, using knowledge of another variable, $x_1$. The relationship between $y$ and $x_1$ is a straight line relationship, with random noise (irreducible error).
  + *Model Form*: A simple linear regression model is a straight line.

+ **Multiple Linear Regression Models:** $\mathbb{E}\left[y\right] = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_k x_k$

  + *General Assumption*: We can predict the response, $y$, using knowledge of other variables, $x_1,~x_2,~\cdots,~x_k$.
    
    + For any of the $x_i$ which are numerical, the association between $y$ and $x_i$ is a straight line relationship, with slope ($\beta_i$) independent of all other predictors $x_j$.
    + For any of the $x_i$ which are dummy variables corresponding to a categorical predictor, the impact of including $x_i$ in the model is a vertical shift. That is, a change in *intercept* specific to that particular level of the categorical variable.
    + All other fluctuation in $y$ is random noise (irreducible error).
  
  + *Model Form*: A multiple linear regression model with no categorical predictors is a *flat* hyperplane, while a multiple linear regression model with categorical predictors is a set of parallel, flat hyperplanes.
  
We can see these model forms below. For ease of visualization, we'll use a simple linear regression model to predict penguin `body_mass_g` using `bill_length_mm` on the left. In the plot on the right, we'll extend that model by introducing dummy variables corresponding to the `species` variable. Note that a hyperplane in two dimensions (`bill_length_mm` and `body_mass_g`) is a straight line.

```{r echo = FALSE, message = FALSE, warning = FALSE}
simp_reg_spec <- linear_reg() %>%
  set_engine("lm")

simp_reg_rec <- recipe(body_mass_g ~ bill_length_mm, data = penguins_train)

simp_reg_wf <- workflow() %>%
  add_model(simp_reg_spec) %>%
  add_recipe(simp_reg_rec)

simp_reg_fit <- simp_reg_wf %>%
  fit(penguins_train)

multi_reg_spec <- linear_reg() %>%
  set_engine("lm")

multi_reg_rec <- recipe(body_mass_g ~ bill_length_mm + species, data = penguins_train) %>%
  step_dummy(species)

multi_reg_wf <- workflow() %>%
  add_model(multi_reg_spec) %>%
  add_recipe(multi_reg_rec)

multi_reg_fit <- multi_reg_wf %>%
  fit(penguins_train)

new_bill_lengths <- seq(min(penguins_train$bill_length_mm, na.rm = TRUE), 
                max(penguins_train$bill_length_mm, na.rm = TRUE),
                length.out = 250)

new_data <- crossing(bill_length_mm = new_bill_lengths, species = c("Adelie", "Chinstrap", "Gentoo"))

new_data <- simp_reg_fit %>%
  augment(new_data) %>%
  rename(.pred_simp = .pred)

new_data <- multi_reg_fit %>%
  augment(new_data) %>%
  rename(.pred_multi = .pred)

p1 <- penguins_train %>%
  ggplot() +
  geom_point(aes(x = bill_length_mm,
                 y = body_mass_g),
             alpha = 0.4) +
  geom_line(data = new_data,
            aes(x = bill_length_mm,
                  y = .pred_simp),
              color = "blue") +
  labs(title = "Predicting Body Mass",
       x = "Bill Length (mm)",
       y = "Body Mass (g)")

p2 <- penguins_train %>%
  ggplot() +
  geom_point(aes(x = bill_length_mm,
                 y = body_mass_g,
                 color = species),
             alpha = 0.4,
             show.legend = FALSE) +
  geom_line(data = new_data,
            aes(x = bill_length_mm,
                  y = .pred_multi,
                  color = species),
            show.legend = FALSE) +
  labs(title = "Predicting Body Mass",
       x = "Bill Length (mm)",
       y = "Body Mass (g)")

(p1 + p2)
```
  
We can see that our models in both cases are straight lines, and that the models are parallel lines in the case where we included the `species` variable. We'll move beyond these restrictions in this notebook and the next one too.

+ In this notebook, we'll discuss methods for adding curvature to our models. This is sometimes called curvi-linear regression.
+ In our next notebook, we'll discuss methods for allowing interaction between terms. Among other things, interaction terms will allow for non-parallel hyperplanes (different rates of change / slopes for different categories).

## Objectives

In previous notebooks we learned how to build, assess, and interpret simple and multiple linear regression models, including those with categorical predictors. Those models all made the following assumptions:

+ The association between the response and any individual predictor is linear.
+ All predictor variables act independently in influencing/predicting the response.
+ The association between any numerical predictor and the response (the *slope* of the model) is the same, regardless of the class/category 

What happens if these assumptions aren't reasonable (or we want to explore whether relaxing these assumptions improves fit)? In this notebook, we'll see how to update our models to allow for curved relationships and predictors which do not act independently to influence our response variable. After reviewing this notebook, you should be able to:

+ Use exploratory data analysis to identify visual evidence for curved relationships and interaction between predictors.
+ Use `step_*()` functions to augment a `recipe()`, updating how predictor variables are utilized in a regression model.

  + These `step_*()` functions are feature engineering steps, which transform our underlying data set.

+ Use `step_poly()` to allow for polynomial terms in a model -- that is, terms whose corresponding predictor variable is raised to a positive integer power.
+ Fit, assess, reduce, interpret, and utilize models including polynomial terms.

***

## A New Data Set: Boston Housing

There's not much visual evidence to suggest curved relationships between the physiological features and penguin body mass in the `palmerpenguins` data set. Because of this, we'll temporarily shift to a new data set on median property values in Boston from the 1970's.

```{r}
boston <- read_csv("https://raw.githubusercontent.com/selva86/datasets/master/BostonHousing.csv")

boston %>%
  head() %>%
  kable() %>%
  kable_styling()
```

This Boston housing dataset is quite famous ([and problematic](https://fairlearn.org/main/user_guide/datasets/boston_housing_data.html)), and includes features on each neighborhood and the corresponding median home value in that neighborhood. You can see a data dictionary [here](https://www.kaggle.com/code/prasadperera/the-boston-housing-dataset). The data set has many interesting features and even allows us some ability to explore structural racism in property valuation in 1970s Boston.

We'll assume that `medv`, the median home value in thousands of dollars, is our response variable and we'll focus on the `age`, `lstat`, and `chas` predictors.

## Motivating Interaction and Higher-Order Terms

Let's take a look at the relationships between some of the available predictors in the `boston` data set and the median home values (`medv`). As usual, we'll split off the *training* and *test* sets so that the *test* data remains hidden from ourselves and the model.

```{r}
boston_split <- initial_split(boston)
boston_train <- training(boston_split)
boston_test <- testing(boston_split)

boston_test %>%
  head() %>%
  kable() %>%
  kable_styling()
```

Now let's do some visualization.

```{r}
p1 <- boston_train %>%
  ggplot() + 
  geom_hex(aes(x = lstat, y = medv)) +
  labs(x = "Proportion Residents with Low Socioeconomic Status",
       y = "Median Home Values (000s)")

p2 <- boston_train %>%
  ggplot() +
  geom_point(aes(x = age, y = medv),
             alpha = 0.5) + 
  labs(x = "Proportion of Homes Built Before 1940",
       y = "Median Home Values (000s)")

p3 <- boston_train %>%
  mutate(chas = ifelse(chas == 1, "on", "off")) %>%
  ggplot() + 
  geom_boxplot(aes(x = chas, y = medv)) +
  labs(x = "On the Charles River",
       y = "Median Home Values (000s)") +
  coord_flip()

(p1 + p2) / p3
```

From the plots above, we see that the association between the proportion of residents having low socioeconomic status (`lstat`) and the median home values (`medv`) have a nonlinear association. Additionally, it may be reasonable to suspect that the "slope" of the model with respect to `age` (the proportion of homes built prior to 1940) is different for neighborhoods on the Charles River versus those which are away from the river. We'll explore all of these things, and more, in the coming sections of this notebook.

## New Terminology

Before we start building models, its worth developing some helpful terminology. Consider the model below:

$\begin{align}\mathbb{E}\left[\text{medv}\right] = &\beta_0 + \beta_1\cdot\left(\text{age}\right) + \beta_2\cdot\left(\text{lstat}\right) + \beta_3\cdot\left(\text{chas}\right) + \beta_4\cdot\left(\text{age}^2\right) +\\ &\beta_5\cdot\left(\text{chas}\cdot\text{age}\right) + \beta_6\cdot\left(\text{chas}\cdot\text{age}^2\right) + \beta_7\cdot\left(\text{chas}\cdot\text{lstat}\right) + \beta_8\cdot\left(\text{age}\cdot\text{lstat}\right) + \beta_9\cdot\left(\text{age}^2\cdot \text{lstat}\right)\end{align}$

+ The **order** of a model is the highest degree of any term in the model. Only numeric predictors contribute to the degree of a term. 

  + The order of the proposed model above is *third* order, because of the $\text{age}^2\cdot\text{lstat}$ term.
+ Models including interaction terms include **main effects** and **mixed effects** terms. The main effects terms contain a single predictor, while the mixed effects terms contain multiple predictors.

$\begin{align}\mathbb{E}\left[\text{medv}\right] = &\overbrace{\beta_0 + \beta_1\cdot\left(\text{age}\right) + \beta_2\cdot\left(\text{lstat}\right) + \beta_3\cdot\left(\text{chas}\right) + \beta_4\cdot\left(\text{age}^2\right)}^{\text{Main Effects}} +\\ &\underbrace{\beta_5\cdot\left(\text{chas}\cdot\text{age}\right) + \beta_6\cdot\left(\text{chas}\cdot\text{age}^2\right) + \beta_7\cdot\left(\text{chas}\cdot\text{lstat}\right) + \beta_8\cdot\left(\text{age}\cdot\text{lstat}\right) + \beta_9\cdot\left(\text{age}^2\cdot \text{lstat}\right)}_{\text{Mixed Effects}}\end{align}$

These *mixed effects* are interaction terms in our model, so we'll wait until the next notebook to work with them.

## Starting with a Main-Effects Model

Let's start out with a basic model which uses `age`, `lstat`, and `chas` as straight-line predictors of median home values. 

```{r}
main_reg_spec <- linear_reg() %>%
  set_engine("lm")

main_reg_rec <- recipe(medv ~ age + lstat + chas, data = boston_train)

main_reg_wf <- workflow() %>%
  add_model(main_reg_spec) %>%
  add_recipe(main_reg_rec)

main_reg_fit <- main_reg_wf %>%
  fit(boston_train)

main_reg_fit %>%
  extract_fit_engine() %>%
  tidy() %>%
  kable() %>%
  kable_styling()
```

Let's compute our performance metrics for our main effects model so that we have a benchmark to try improving from.

```{r}
my_metrics <- metric_set(rsq, rmse)

(main_reg_fit %>%
    augment(boston_train) %>%
    my_metrics(medv, .pred) %>%
    mutate(type = "train")
    ) %>%
  bind_rows(
    main_reg_fit %>%
      augment(boston_test) %>%
      my_metrics(medv, .pred) %>%
      mutate(type = "test")
    ) %>%
  pivot_wider(id_cols = .metric,
              names_from = type,
              values_from = .estimate) %>%
  kable() %>%
  kable_styling()
```

Okay, now that we've got those benchmark metrics, let's see if we can improve on them by using higher-order terms.

## Building a Model Including Curvature (`step_poly()`)

In order to build a model including curvature, we'll need to include some *polynomial* terms. These are simply terms which include a predictor raised to an integer power greater than $1$. We can do this by adding a *feature engineering* step to our `recipe()`, just like we did with `step_dummy()`. In order to add these higher-order terms that introduce curvature, we'll use `step_poly()`. Let's see it in action! We'll include second-order terms associated with `age` and `lstat`, but not for the `chas` variable. This is because `chas` is a dummy variable, which can only take the values $0$ or $1$ -- what happens when you square these values? Including polynomial terms with dummy variables doesn't actually do anything!

Before we build our model, there are just a few comments about `step_poly()` that we should highlight.

+ If you'd like to build polynomial terms associated with multiple variables, you can do this all at once by including each of the variables as the first few arguments to `step_poly()`, separating the variables by commas.

+ Set the `degree` of the polynomial terms (the amount of *wiggling* you'd like to allow) with the `degree` argument.

  + The default degree for `step_poly()` is degree $2$, allowing a single bend in the association between your predictor(s) and the response. You can override this default by setting `degree = 6`, for example. This would allow up to five bends in the association.

+ Maintain your ability to interpret the relationship between predictor and response by setting `options = list(raw = TRUE)`

  + There is some *linear algebra* that can be done to try and reduce the correlation between our polynomial terms and their corresponding lower-order terms. You can think of this as squeezing all of the predictive juice out of the predictor. The `step_poly()` function does this by default, but results in a model which is more difficult to interpret. I recommend overriding this any time you intend to interpret your models. You can do this by setting `options = list(raw = TRUE)` within `step_poly()`.

```{r}
poly_reg_spec <- linear_reg() %>%
  set_engine("lm")

poly_reg_rec <- recipe(medv ~ age + lstat + chas, data = boston_train) %>%
  step_poly(age, degree = 2, options = list( raw = TRUE)) %>%
  step_poly(lstat, degree = 2, options = list(raw = TRUE))

poly_reg_wf <- workflow() %>%
  add_model(poly_reg_spec) %>%
  add_recipe(poly_reg_rec)

poly_reg_fit <- poly_reg_wf %>%
  fit(boston_train)

poly_reg_fit %>%
  glance() %>%
  kable() %>%
  kable_styling()

poly_reg_fit %>%
  extract_fit_engine() %>%
  tidy()  %>%
  kable() %>%
  kable_styling()
```

Note that the term for `age_poly_2` is not statistically significant. We'll drop it from the model and re-run the regression.

```{r}
poly_reg_spec <- linear_reg() %>%
  set_engine("lm")

poly_reg_rec <- recipe(medv ~ age + lstat + chas, data = boston_train) %>%
  step_poly(lstat, degree = 2, options = list(raw = TRUE))

poly_reg_wf <- workflow() %>%
  add_model(poly_reg_spec) %>%
  add_recipe(poly_reg_rec)

poly_reg_fit <- poly_reg_wf %>%
  fit(boston_train)

poly_reg_fit %>%
  glance() %>%
  kable() %>%
  kable_styling()

poly_reg_fit %>%
  extract_fit_engine() %>%
  tidy()  %>%
  kable() %>%
  kable_styling()
```

There; now all of our model terms are statistically significant. Let's see what that model looks like!

```{r}
poly_reg_fit %>%
  augment(boston_train) %>%
  ggplot() +
  geom_point(aes(x = lstat, y = medv),
             alpha = 0.5) + 
  geom_point(aes(x = lstat, y = .pred), 
             color = "red",
             alpha = 0.4) + 
  labs(x = "Proportion of Residents with Low Socioeconomic Status",
       y = "Median Neighborhood Home Value (000s)",
       title = "Median Home Values and Socioeconomic Status",
       subtitle = "Actuals in Black, Predicted Values in Red")
```

We can also look at our model for a few `age` thresholds.

```{r}
new_data <- crossing(age = c(10, 25, 50),
                     chas = c(0, 1),
                     lstat = seq(min(boston_train$lstat, na.rm = TRUE),
                                max(boston_train$lstat, na.rm = TRUE,
                                    by = 1)))

poly_reg_fit %>%
  augment(new_data) %>%
  ggplot() + 
  geom_line(aes(x = lstat, 
                y = .pred, 
                color = as.factor(age),
                linetype = as.factor(chas)),
            lwd = 1) +
  labs(x = "Proportion of Residents with Low Socioeconomic Status",
       y = "Predicted Median Home Value",
       title = "Estimated Median Home Values",
       subtitle = "(by prevalence of low socioeconomic status)",
       color = "% Homes Built \nPrior to 1940",
       linetype = "On Charles River")
```



This is a model we can interpret and utilize. Before we do, including higher-order terms in a model requires some additional care.

**Note:** Since `lstat_poly_2`, the term corresponding to `lstat`$^2$, is statistically significant, we must keep the lower order term for `lstat` in the model. In general, if we have a statistically significant higher-order, or interaction term in a model, then we *must* retain all of the associated lower-order terms in the model, regardless of whether or not they are statistically significant. For example, if we had a term for `lstat`$\cdot$`age`$^2$ in our model, then we must keep the individual terms for `lstat` and `age`, the term for `age`$^2$ and the term for `lstat`$\cdot$`age`, regardless of whether the $p$-values associated with those terms are small.

### Interpreting a Model with Curvature

```{r}
poly_reg_fit %>%
  extract_fit_engine() %>%
  tidy() %>% 
  kable() %>%
  kable_styling()
```

Now that we have estimated coefficients, we can write down the estimated model form:
$$\mathbb{E}\left[\text{medv}\right] \approx 40.61 + 0.08\cdot\text{age} + 4.03\cdot\text{chas} - 2.68\cdot\text{lstat} + 0.05\cdot\text{lstat}^2$$

From here, we see that older neighborhoods tend to have slightly higher median home values. Neighborhoods along the Charles River have higher median home values. Before moving to interpret the `lstat` variable, let's formalize the interpretations associated with `age` and `chas`.

+ Holding all variables other than `age` constant (controlling for them), an increase by $1$ percentage point in the proportion of homes in a neighborhood built before 1940 is associated with higher median home values by about $\$80$, on average.

  + Alternatively, given two otherwise similar neighborhoods (`lstat` and `chas`), the neighborhood with a larger percentage of homes built before 1940 is expected to have a higher median home value by about $\$80$ per percentage point, on average.

+ Holding all variables other than `chas` constant (controlling for them), median home values for neighborhoods on the Charles River are expected to be about $\$4,030$ higher, on average.

  + Alternatively, given two otherwise similar (`age` and `lstat`) neighborhoods, the one on the Charles River is expected to have approximately $\$4,030$ higher median home values, on average, than the neighborhood away from the Charles River.

Because the effect of `lstat` is split across multiple terms, we need some way to aggregate that effect. As a reminder, we are looking for *rate of change* in `medv` with respect to `lstat`. Those of you with a calculus background will recognize that we are interested in finding the partial derivative of our model with respect to `lstat`. For those of you without a calculus background, you'll just use the expression below:

\begin{align} &\text{Rate of Change in medv With Respect to lstat: }\\ &~~~~\beta_{\text{lstat}} + 2\beta_{\text{lstat}^2}\cdot \text{lstat} = -2.68 + 2\left(0.05\right)\cdot\text{lstat}\end{align}

The rate of change in median home values (`medv`) with respect to the percentage of residents with low socioeconomic status now depends on the level of the percentage of residents with low socioeconomic status. The associated drop in median home values is large for percentage point increases from relatively low prevalence of residents with low socioeconomic status, but then for neighborhoods with already high levels of residents with low socioeconomic status, the decrease in median home values is less. For example, the difference in median home values between a neighborhood with 10% of residents having low socioeconomic status and a neighborhood with 11% of residents having low socioeconomic status is expected to be about $-2.68 + 0.1\cdot 10$. That is, we expect lower median home values by about $\$1,680$, on average. The difference in expected median home values between a neighborhood with 20% of residents having low socioeconomic status and a neighborhood with 21% though is expected to be about $-2.68 + 0.1\cdot 20$, or about $\$680$ less, on average.

The inclusion of these polynomial terms allows us to model more complex relationships and potentially build models with better predictive value! There are a few things to beware of though -- these models are more complicated to interpret and we risk tricking ourselves into believing our models are much more accurate than they will be in production (More on this later in our course).

#### Interpretation with `{marginaleffects}`

As we saw in an earlier notebook, we can use the `{marginaleffects}` package to estimate the marginal effect of a unit increase in `lstat` on median home values at different levels of the `lstat` variable. We can also obtain confidence interval estimates in these effects as well.

```{r}
counterfactual_lstat_df <- tibble(lstat = seq(0, 40, length.out = 500),
                                  age = 50,
                                  chas = 1)

mfx <- poly_reg_fit %>% 
  slopes(newdata = counterfactual_lstat_df,
         variables = "lstat") %>%
  tibble() %>%
  #Not sure why, but slopes() is obliterating lstat values
  mutate(lstat = seq(0, 40, length.out = 500))

mfx %>%
  select(term, lstat, estimate, conf.low, conf.high, std.error) %>%
  head(n = 10) %>%
  kable() %>%
  kable_styling()

mfx %>%
  select(lstat, estimate, conf.low, conf.high) %>%
  mutate(lstat = seq(0, 40, length.out = 500)) %>%
  ggplot() +
  geom_line(aes(x = lstat, y = estimate), 
            color = "purple", lty = "dashed", lwd = 1.5) +
  geom_ribbon(aes(x = lstat, ymin = conf.low, ymax = conf.high),
              fill = "grey", alpha = 0.5) +
  labs(x = "Pct Low Socioeconomic Status",
       y = "Marginal Effect",
       title = "Marginal Effects of Unit Increase in lstat")
```

As we can see from the plot above, the marginal effect of a unit increase in the percentage of residents with low socioeconomic status is expected to be a decrease in median home values up until low socioeconomic status residents make up about 27% of the neighborhood population. After this point, we expect the marginal effect to be an increase in median home values. We see this also from our earlier plot of estimated median home values by prevalence of low socioeconomic status, plotted at different thresholds of the `age` and `chas` variables.

### Assessing Model Performance

Before we move on, let's assess the performance of our new model. As a reminder, we had the following benchmarks from the *main effects* model.

```{r}
my_metrics <- metric_set(rsq, rmse)

(main_reg_fit %>%
    augment(boston_train) %>%
    my_metrics(medv, .pred) %>%
    mutate(type = "train")
    ) %>%
  bind_rows(
    main_reg_fit %>%
      augment(boston_test) %>%
      my_metrics(medv, .pred) %>%
      mutate(type = "test")
    ) %>%
  pivot_wider(id_cols = .metric,
              names_from = type,
              values_from = .estimate) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

We can compute the same metrics for our new model.

```{r}
(poly_reg_fit %>%
    augment(boston_train) %>%
    my_metrics(medv, .pred) %>%
    mutate(type = "train")
    ) %>%
  bind_rows(
    poly_reg_fit %>%
      augment(boston_test) %>%
      my_metrics(medv, .pred) %>%
      mutate(type = "test")
    ) %>%
  pivot_wider(id_cols = .metric,
              names_from = type,
              values_from = .estimate) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

We are doing quite a bit better now! The R-Squared value for our new model on the training data is up to 70% (from 58% with the main effects model). We've also observed an improvement in R-Squared on the test data -- 62% versus 52%. Much more of the variation in median home values is being explained by our curvi-linear model than was explained with our main effects model. Similarly, we've observed improvements in the RMSE (*root mean squared error*). The test RMSE for our main effects model was approximately $6.4$, indicating that we could expect our predictions to be accurate to within about $\$12,800$ on average. The test RMSE for our curvi-linear model is approximately 5.8, indicating that we should expect our predictions to be accurate to within about $\$11,600$ on average.

*** 

## Summary

Okay, let's recap. We worked with a new data set today -- the `boston` housing data set. For this particular application, we were interested in predicting median home values (`medv`) for neighborhoods in 1970's Boston, given 

+ `age` -- the proportion of homes in the neighborhood built prior to $1940$.
+ `lstat` -- the proportion of neighborhood residents having low socioeconomic status.
+ `chas` -- an indicator (dummy) variable for whether the neighborhood is along the Charles River or not.

In our initial exploratory analysis, we saw visual evidence that the association between `medv` and `lstat` was not a *straight-line* association. We suspected that a model allowing for a curved relationship between `lstat` and `medv` would outperform a model that assumed a straight-line relationship. We went ahead and built a model assuming straight line relationships between `medv` and both `age` and `lstat`, but which also accounted for the location of the neighborhood being on the Charles River or away from the Charles River. The terms in that model were statistically significant, so we obtained some performance metrics for this model and used those as a *benchmark*.

We then built a model allowing for curved relationships between `medv` and both `age` and `lstat`. There wasn't statistically significant evidence to keep the curvature between `medv` and `age`, so we dropped it and returned to a linear relationship between `medv` and `age`. There was, however, statistically significant evidence for a curved relationship between `medv` and `lstat`. We interpreted our new model and calculated performance metrics for it, noting that the model allowing for curvature out-performed our original model on the unseen test data.

In the next notebook, we'll return to the `penguins` data and see what interaction terms allow for in a model.