linear_model_extensions.qmd

# Extending the Linear Model {#sec-lm-extend}

![](img/chapter_gp_plots/gp_plot_5.svg){width=75%}

```{r}
#| include: false
#| label: r-data-import

df_reviews = read_csv('data/movie_reviews.csv')
df_happiness = read_csv('data/world_happiness_2018.csv')
df_happiness_all = read_csv('data/world_happiness_all_years.csv')

model_baseline = lm(rating ~ children_in_home + genre, data = df_reviews)
model_interaction = lm(rating ~ children_in_home * genre, data = df_reviews)
```

```{python}
#| include: false
#| label: py-data-import

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

np.set_printoptions(formatter={'float': lambda x: '{0:0.4f}'.format(x)})
pd.set_option('display.float_format', lambda x: '%.3f' % x)

df_reviews = pd.read_csv('data/movie_reviews.csv')
df_happiness = pd.read_csv('data/world_happiness_2018.csv')
```

With just linear and generalized linear models, we have a very solid foundation for modeling, and we've seen how there is a notable amount we can do with a conceptually simple approach. We've also seen how we can extend the linear model to handle different types of target distributions to help us understand and make some inferences about the relationships between our features and target.

In this chapter, we will extend our linear models with additional common and valuable modeling tools. These methods provide good examples of how we can think about our data and modeling approach in different ways, and they can serve as a foundation for exploring more advanced techniques in the future. A thread that binds these techniques together is the ability to use a linear model to investigate explicitly nonlinear relationships!


## Key Ideas {#sec-lm-extend-key-ideas}

- The standard and generalized linear models are great and powerful starting points for modeling, but there's even more we can do!
- Linear models can be used to model nonlinear feature-target relationships!
- While these seem like different approaches, we can still use our linear model concepts and approach at the core, take similar estimation steps, and even have similar interpretation. However, we'll have even more results to explore and interpret.


### Why this matters {#sec-lm-extend-why}

The linear model is a great starting point for modeling. It is a simple approach that can be used to model a wide variety of relationships between features and targets, and it's also a great way to get a feel for how to think about modeling. But linear and generalized models are just the beginning, and the models depicted here are common extensions used in a variety of disciplines and industries. More generally, the following techniques allow for nonlinear relationships while still employing a linear model approach. This is a very powerful tool to have in your toolkit, and it's a great way to start thinking about how to model more complex relationships in your data.


### Helpful context {#sec-lm-extend-good}

While these models are extensions of the linear model, they are not significantly more complicated in terms of how they are implemented or how they are interpreted. However, like anything new, it can take a bit more effort to understand. You likely want to be comfortable with standard linear models at least before you start to explore these extensions.


## Interactions {#sec-lm-extend-interactions}

Things can be quite complex in a typical model with multiple features, but just adding features may not be enough to capture the complexity of the relationships between features and target. Sometimes, we need to consider how features interact with each other to better understand how they correlate with the target. A common way to add complexity in linear models is through **interactions**. This is where we allow the effect of a feature to vary depending on the values of another feature, or even itself! 

As a conceptual example, we can think about a movie's rating being different for movies from different genres. For example, maybe by default ratings are higher for kids movies, and lower for horror movies. But, genre and season might work together in some way to affect rating, e.g., action movies get higher ratings in summer. Or maybe having kids in the home might also interact with genre ratings by naturally resulting in higher ratings for kids movies. As a different example, we might also consider that the length of a movie might positively relate to rating, but plateau or even have a negative effect on rating after a certain point. In other words, it would have a **curvilinear** effect where really long movies aren't as highly rated as those of shorter length. 

All of these are types of interactions we can explore. Interactions allow us to incorporate nonlinear relationships into the model, and so greatly extend the linear model's capabilities - we basically get to use a linear model in a nonlinear way!  

With that in mind, let's explore how we can add interactions to our models. Going with one of our examples, let's see how having kids impacts the relationship between genre and rating. We'll start with a standard linear model, and then add an interaction term. Using a formula approach makes it very straightforward to add an interaction term. We just need to add a `:` between the two features we want to interact, or a `*` to denote both main effects and the interaction. As elsewhere, we present simplified results in the next table.

:::{.panel-tabset}

##### R

```{r}
#| label: model-interaction-r
#| output: false
df_reviews = read_csv('https://tinyurl.com/moviereviewsdata')

model_baseline = lm(rating ~ children_in_home + genre, data = df_reviews)
model_interaction = lm(rating ~ children_in_home * genre, data = df_reviews)

summary(model_interaction)
```

##### Python

```{python}
#| echo: true
#| label: model-interaction-py
#| eval: false

import pandas as pd
import statsmodels.formula.api as smf

df_reviews = pd.read_csv('https://tinyurl.com/moviereviewsdata')

model_baseline = smf.ols(
    formula = 'rating ~ children_in_home + genre', 
    data = df_reviews
).fit()

model_interaction = smf.ols(
    formula = 'rating ~ children_in_home * genre', 
    data = df_reviews
).fit()

model_interaction.summary()
```

:::


```{r}
#| echo: false
#| label: model-interaction-get-reference-group
ref = model_baseline$xlevels$genre[1]
```

Here is a quick look at the model output for the interaction vs. no interaction interaction model. Starting with the baseline model, the coefficients look like what we've seen before, but we have several coefficients for genre. The reason is that genre is composed of several categories, and converted to a set of dummy variables (refer to @sec-lm-categorical-features and @sec-data-cat). In the baseline model, the intercept tells us what the mean is for the reference group, in this case `r ref`, and the genre coefficients tell us the difference between the mean for that genre and the reference genre. For example, the mean rating for `r ref` is `r round(coef(model_baseline)['(Intercept)'], 2)`, and the difference between that genre rating for the drama genre is `r round(coef(model_baseline)['genreDrama'], 2)`. Adding the two gives us the mean for drama movies `r round(coef(model_baseline)['(Intercept)'], 2)` + `r round(coef(model_baseline)['genreDrama'], 2)` = `r round(coef(model_baseline)['(Intercept)'] + coef(model_baseline)['genreDrama'], 2)`. We also have the coefficient for the number of children in the home, and this does not vary by genre in the baseline model.

```{r}
#| echo: false
#| label: tbl-model-interaction-output
#| tbl-cap: Model Coefficients with Interaction

library(broom)

tidy_base = tidy(model_baseline) |> 
    rename(
        feature = term,
        coef = estimate,
        se = std.error,
        t = statistic,
        p = p.value
    )

tidy_interaction = tidy(model_interaction) |> 
    rename(
        feature = term,
        coef = estimate,
        se = std.error,
        t = statistic,
        p = p.value
    )

nchildren_coef = round(coef(model_interaction)['children_in_home'], 2)
sign_nchildren_coef = sign(nchildren_coef)
genre_kids_coef = round(coef(model_interaction)['genreKids'], 2)
nkids_comedy_coef = round(coef(model_interaction)['children_in_home:genreComedy'], 2)
nkids_kids_coef = round(coef(model_interaction)['children_in_home:genreKids'], 2)
sign_nkids_kids_coef = sign(nkids_kids_coef)
# mean(df_reviews$children_in_home == 0)
# mean(df_reviews$genre == ref)

# Create a gt table from the tidy coefficients
# COLOR CODING WON'T WORK IN LATEX BECAUSE IT SUCKS
tidy_base |> 
    full_join(tidy_interaction, by = 'feature', suffix = c('_base', '_inter')) |>
    select(feature, starts_with('coef')) |>
    gt(decimals = 3) |> 
    sub_missing(
        missing_text = "",
    ) |> 
    # gt_highlight_rows(
    #   columns = vars(coef_base, coef_inter),
    #   rows = str_detect(feature, 'children_in_home$|children.*Kids'),
    #   fill = okabe_ito[5],
    #   alpha = .1
    # ) |>
    tab_caption('Model coefficients with and without an interaction') 
```

But in our other model we have an interaction between two features: 'children in the home' and 'genre'. So let's start with the coefficient for children. It is `r nchildren_coef`, which means that for every additional child, the rating for any movie `r ifelse(sign_nchildren_coef, 'increases', 'decreases')` by that amount. But because of the interaction, we now interpret that as the *effect of children when genre is the reference group `r ref`*. 

Now let's look at the interaction effect for children and the kids genre. It is `r nkids_kids_coef`, which means that for the kids genre, the *effect of having children in the home `r ifelse(sign_nkids_kids_coef, 'increases', 'decreases')`* by that amount. So our actual effect for an additional child for the kids genre is  `r nchildren_coef` + `r nkids_kids_coef` = `r round(nchildren_coef + nkids_kids_coef, 2)` `r ifelse(sign_nchildren_coef + sign_nkids_kids_coef, 'increase', 'decrease')` in the expected review rating. In other words, the effect of children in the home is stronger for kids movies than for other genres.

We can also interpret this interaction in the reverse fashion. It is also correct to say that the difference in rating between the kids genre and the reference group `r ref` when there are no kids in the home is `r genre_kids_coef`. This means kids movies are generally rated worse if there are no kids in the home. But, with an increase in children, the difference in rating between the kids genre and `r ref` `r ifelse(sign_nkids_kids_coef, 'increases', 'decreases')` by `r abs(nkids_kids_coef)`. In other words, it is a **difference in differences**[^diffindiff].

[^diffindiff]: Some models that employ an interaction that investigates categorical group differences like this actually call their model a difference-in-difference model.

When we talk about differences in coefficients across values of features, it can get a little bit hard to follow. To combat this, we believe that in every case that you employ an interaction, you should look at the interaction visually for interpretation. Here is a plot of the predictions from the interaction model. We can see that the effect of children in the home is stronger for kids movies than for other genres, which makes a whole lot of sense! In other genres, the effect of having children seems to produce little difference, and in others it still has a positive effect, but not as strong as for kids movies.

::: {.content-visible when-format='html'}
```{r }
#| echo: false
#| label: fig-model-interaction-plot
#| fig-cap: Interaction Plot

library(ggeffects)

p_init = plot(
    ggpredict(model_interaction, terms = c('children_in_home', 'genre')), 
    alpha = .05, 
    # line_size = I(rep(c(rep(.10, 7), 2), 4)) # values of children in home = 0:3, kids is last genre, but it's ignored anyway
    line_size = 1
    )$data |>
    as_tibble()

p_init |> 
    ggplot(aes(x = x, y = predicted, color = group, fill = group)) +
    geom_line(linewidth = 1.5) +
    # geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1, color = NA) +
    geom_borderline(linewidth = 2.5, data = p_init |> filter(group == 'Kids') ) +
    scale_color_manual(values = unname(okabe_ito), aesthetics = c('color', 'fill')) +
    # scale_color_brewer(palette = 'Set1', aesthetics = c('color', 'fill')) +
    scale_y_continuous(breaks = seq(0, 5, .25)) +
        labs(
        x = 'Children in home', 
        y = 'Predicted\nrating',
        caption = 'Bold line is the Kids genre'
    ) +
    theme_clean() +
    theme(
        panel.grid.major = element_line(color = alpha('black', .05)),    
    )


p_print = p_init |>     
    filter(group != 'Kids') |> 
    ggplot(aes(x = x, y = predicted)) +
    geom_line(aes(linetype = group), linewidth = 1) +
    geom_line(
        aes(linetype=group),
        linewidth = 3, 
        linetype = 'solid',
        data = p_init |> filter(group == 'Kids') |> droplevels()
    ) +
    scale_y_continuous(breaks = seq(0, 5, .25)) +
    labs(
        x = 'Children in home', 
        y = 'Predicted\nrating',
        caption = 'Bold line is the Kids genre'
    ) +
    theme(
        legend.key.size = unit(1.5, 'cm'),
        panel.grid.major = element_line(color = alpha('black', .05)),    
    )

# plot(modelbased::estimate_expectation(model_interaction), si_alpha = .01)
ggsave('img/lm-extend-interaction_plot.svg', p_print, width = 8, height = 6)
```
:::

::: {.content-visible when-format='pdf'}
![Interaction Plot](img/lm-extend-interaction_plot.svg){#fig-model-interaction-plot}
:::

So we can see that interactions can allow a linear effect to vary depending on the values of another feature. But the real take home message from this is that *the general effect is actually not just a straight line*! The linear effect *changes* depending on the setting. Furthermore, the effect for children when the interaction is present only applies when 'genre' is at its default group, or when other features are at their default or zero.

So, when we have interactions, we can't talk about a feature's relationship with a target without considering the other features it interacts with. Some might see this as a downside, but it's actually how most feature relationships work - they don't exist in isolation. Interactions let us model these complex relationships, and they're used a lot in real-world situations.

### Summarizing Interactions

So what is *the* effect of children in the home? Or a particular genre, for that matter? This is a problematic question, because the effect of one feature depends on the setting of the other feature. We can summarize interactions, and we show two ways in which to do so. But we need to be very careful about summarizing a single feature's effects when we know it interacts with another feature.

### Average effects {#sec-lm-extend-avgeffect}

One thing we can do is get an average effect for a feature. In other words, we can say what the effect of a feature is *on average* across the settings of the other features. This is called the **average marginal effect**, something we've talked about in @sec-avg-marginal-effects[^ME]. Here is the average effect of children in the home across all genres.

[^ME]: These results are provided by the [marginaleffects]{.pack} R package, which is a great tool.

\small
```{r}
#| echo: false
#| label: tbl-model-interaction-avg-effect
#| tbl-cap: Average Marginal Effects of Children in the Home
library(marginaleffects)

me_children = marginaleffects::avg_slopes(model_interaction, variables = c("children_in_home"))

me_children |> 
    as_tibble() |> 
    select(-s.value) |> 
    gt() |> 
    fmt_number(
        columns = c(estimate), 
        decimals = 3
    )
```

\normalsize


So-called **marginal effects**, and related approaches such as SHAP values (see @sec-knowing-shap-values), attempt to boil down the effect of a feature to a single number. Here we see the average coefficient for children in the home is `r round(me_children$estimate, 2)`, with a range from `r round(me_children$conf.low, 2)` to `r round(me_children$conf.high, 2)`. 

But this is difficult even in the simpler GLM settings, and can be downright misleading in interaction models. We saw from @tbl-model-interaction-output that this average is slightly larger than what we would estimate in the baseline model, and we saw in @fig-model-interaction-plot it's actually near zero (flat) for some genres in the interaction model. So what is the average effect really telling us? Consider a more serious case of drug effects across demographic groups, where the effect of the drug is much stronger for some groups than others. Would you want your doctor to prescribe you a drug based on the average effect across all groups, or the specific group to which you belong?




### ANOVA {#sec-lm-extend-anova}

A common method for summarizing categorical effects in linear models is through **analysis of variance** or ANOVA, something we briefly mentioned in our chapter introducing the linear model @sec-lm-categorical-anova. ANOVA breaks down the variance in a target attributable to different features or their related effects such as interactions. It's a bit beyond the scope here to get into all the details, but we demonstrate it here, as it's also used to summarize interactions. It also summarizes the random effects and spline terms we'll see in the coming sections on mixed and generalized additive models.

:::{.panel-tabset}

##### R


```{r}
#| label: anova-r
#| eval: false

anova(model_interaction)
```

##### Python

```{python}
#| label: anova-py
#| eval: false
sm.stats.anova_lm(model_interaction)
```

:::



```{r}
#| echo: false
#| label: tbl-anova-r
#| tbl-cap: ANOVA Table for an Interaction Model
anova(model_interaction) |> 
    tidy() |> 
    rename(
        feature = term,
        df = df,
        sum_sq = sumsq,
        mean_sq = meansq,
        f = statistic,
        p = p.value
    ) |>
    gt() |> 
    sub_missing(
        missing_text = "",
    )
```

In this case, it doesn't appear that the general interaction effect is statistically significant if we use the typical .05 cut-off. We know the effect of children in the home varies across genres, but this result suggests maybe it's not as much as we might think. However, we also saw that the estimate for the children effect more than doubled for the kids genre, so maybe we don't want to ignore it. That's for you to decide.

The ANOVA approach can be generalized to provide a statistical test to compare models. For example, we can compare the baseline model to the interaction model to see if the interaction model is a better fit. However, it's entirely consistent with just looking at the interaction's statistical result in the ANOVA for the interaction model alone, so it doesn't provide any additional information in this case. Note that the only models that can be compared with ANOVA must be *nested*, i.e., one model is an explicit subset of the other.

:::{.callout-note title='ANOVA is not a model' collapse='true'}
It's worth noting that ANOVA is often confused with being a model itself. When people use it this way, it is just a linear regression with only categorical features, something that is usually only seen within strict experimental designs. It's pretty difficult to think of a linear regression setting where no continuous features would be of theoretical interest, but back when people were doing this stuff by hand, they just categorized everything to enable doing an ANOVA, which was tedious arithmetic but more manageable. It's a bit of a historical artifact, but might be useful for exploratory purposes of feature relationships. For model comparison however, other approaches are more general than ANOVA and generally preferred. An example is using AIC, or the other methods employed that we discuss the in machine learning chapters, like cross-validation error.
:::



### Interactions in Practice

When dealing with interactions in a model, it's best to consider how the feature-target relationship changes based on the values of other features it interacts with. Visualizing these effects is important in helping us understand how the relationships change. It's also helpful to consider what the predicted outcome is at important feature values, and how this changes with different feature values. This is the approach we've used with interactions, and it's a good strategy overall.



## Mixed Models {#sec-lm-extend-mixed-models}


### Knowing your data {#sec-lm-extend-mixed-knowing}

As much fun as modeling is, knowing your data is far more important. You can throw any model you want at your data, from simple to fancy, but you can count on disappointment if you don't fundamentally know the structure that lies within your data. Let's take a look at the following visualizations. In @fig-length-release-rating, we see a positive relationship between the length of the movie and ratings. 


:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-length-release-rating
#| fig-cap: Linear Relationship Between Length of Movie and Rating
ggplot(df_reviews, aes(length_minutes, rating)) +
    geom_point() +
    geom_smooth(method="lm", se = FALSE) +
    labs(
        x = 'Length of movie (minutes)',
        y = 'Rating',
        # caption = 'Data from movie reviews'
    )

ggsave("img/lm-extend-length_release_rating.svg", width = 8, height = 6)
```
:::
:::{.content-visible when-format='pdf'}
![Linear Relationship Between Length of Movie and Rating](img/lm-extend-length_release_rating.svg){#fig-length-release-rating}
:::



We could probably just stop there, but given what we just saw with interaction, we might think to ourselves to be ignoring something substantial within our data: genre. We might want to ask a question, "Does this relationship work the same way across the different genres?"


```{r}
#| echo: false
#| eval: false
#| label: fig-length-genre-rating
#| fig-cap: Genre Effects on Length and Rating
#| out-width: 100%

# p_mod = lm(rating ~ length_minutes*genre, 
#           df_reviews |> mutate(genre = recode(genre, "Action/Adventure" = "Act/Adv")))

# plot(modelbased::estimate_expectation(p_mod, trend = 'length_minutes', at = 'genre')) +
#   scale_x_discrete(guide = guide_axis(n.dodge=2)) +
#   labs(title = '', y = '') +
#   theme_clean()

mods = df_reviews |> 
    nest_by(genre) |>
    mutate(model = list(lm(rating ~ length_minutes, data = data))) |>
    reframe(broom::augment(model))

p_slopes = ggplot(mods, aes(length_minutes, .fitted, color = genre)) +
    geom_line(aes(linetype = genre), lwd = .5) +
    scale_y_continuous(limits = c(1, 5)) +
    labs(title = '', y = '', x = '') 

p_mod = lm(rating ~ genre, df_reviews |> mutate(genre = recode(genre, "Action/Adventure" = "Act/Adv")))

p_dat = plot(modelbased::estimate_means(p_mod))$data

p_group_means = p_dat |> 
    ggplot(aes(x = genre, y = Mean)) +
    geom_dotplot(
        aes(x = genre, y= rating),
        binaxis='y', 
        stackdir='center',
        dotsize=.1,
        color = 'black',
        data = df_reviews |> 
        mutate(genre = recode(genre, "Action/Adventure" = "Act/Adv")), 
        position = position_jitter(width = .0, height = .15),
        alpha = .25
    ) +
    geom_pointrange(
        aes(ymin = CI_low, ymax = CI_high, color = genre),
        fatten = 5,
        alpha = 1,
        show.legend = FALSE
    ) +
    scale_x_discrete(guide = guide_axis(n.dodge=2)) +
    labs(title = '', y = '', x = '') +
    scale_y_continuous(limits = c(1, 5))

see::plots(
    p_group_means,
    p_slopes
) + 
    plot_layout(guides = 'collect') &
    # scale_color_brewer(palette = 'Set1') &
    plot_annotation(
        title = 'Genre Effects on Length and Rating',
        caption = 'Left Plot: Group Means by Genre\nRight Plot: Length Slopes by Genre'
    ) &
    theme(
        panel.grid.major.y = element_line(color = scales::alpha('black', .1)),
        legend.key.size = unit(1.25, "cm"),
    )

ggsave("img/lm-extend-length_genre_rating.svg", width = 8, height = 6)
```


![Genre Effects on Length and Rating](img/lm-extend-length_genre_rating.svg){#fig-length-genre-rating}

A very quick examination of @fig-length-genre-rating might suggest that the rating varies by genre, and that the relationship between length and rating varies significantly over the different genres. The group means in the right panel show variability across genres. In addition, on the left panel, some genres show a strong positive relationship, some show less of a positive relationship, a couple even show a negative relationship, and one even looks flat. We can also see that they would have different intercepts. This is a very important thing to know about your data! If we had just run a model with length as a feature and nothing else, we would have missed this important information.

Now consider something a bit more complicated. Here is a plot of the relationship between the length of a movie and the rating, but across release year.  Again we might think there is notable variability in the effect across years, as some slopes are positive, some very strongly positive, and some are even negative.  How can we account for this when there are so many group effects to consider?

:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-length-release-rating-year
#| fig-cap: Release Year Effects on Length and Rating

# glimpse(df_reviews)
mods_release = df_reviews |> 
    nest_by(release_year) |>
    mutate(model = list(lm(rating ~ length_minutes, data = data))) |>
    reframe(broom::augment(model))

p_slopes_release = ggplot(mods_release, aes(length_minutes, .fitted, color = release_year)) +
    geom_line(aes(group= release_year), lwd = 1) +
    scale_color_scico(palette = 'oslo', end = .85, breaks = seq(1985, 2020, 10)) +
    scale_y_continuous(limits = c(1, 5)) +
    labs(title = '', y = '', x = '') 

p_slopes_release

ggsave("img/lm-extend-length_release_rating_year.svg", width = 8, height = 6)
```
:::
:::{.content-visible when-format='pdf'}
![Release Year Effects on Length and Rating](img/lm-extend-length_release_rating_year.svg){#fig-length-release-rating-year}
:::


### Overview of mixed models {#sec-mixed-models-overview}

What we've just seen might initially bring to mind an interaction effect, and that's the right way to think about it! A **mixed model** can be used to incorporate that type of relationship into our model, which we can think of as a group interaction, without much hassle and additional explainability. But it's actually a quite flexible class that can also allow for more complicated but related types.

Before going too much further, the term *mixed model* is as vanilla as we can possibly make it, but you might have heard of different names such as *hierarchical linear models*, or *multilevel models*, or maybe *mixed-effects models*. Maybe you've even been exposed to ideas like *random effects* or *random slopes*. These are in fact all instances of what we're calling a *mixed model*. 

What makes a model a *mixed* model? The mixed model is characterized by the idea that a model can have **fixed effects** and **random effects**. Fortunately, you've already encountered *fixed* effects -- those are the features that we have been using in all of our models so far! We are assuming a single true parameter (e.g., coefficient/weight) to estimate for each of those features, and that parameter is *fixed*.

In mixed models, a *random effect* instead comes from a specific distribution, and this is almost always a normal distribution. This random effect adds a unique source of variance in the target variable. This distribution of effects can be based on a grouping variable (such as genre), where we let those parameters, i.e., coefficients (or weights), vary across the groups, creating a distribution of values.

Let's take our initial example with movie length and genre. Formally, we might specify something like this:

$$
\text{Rating} = b_{\text{int[genre]}} + b_\text{length}*\text{length} + \epsilon
$$

In this formula, we are saying that genre has its own unique effect for this model in the form of specific intercepts for each genre[^mixedeps]. This means that whenever an observation belongs to a specific genre, it will have an intercept that reflects that genre, and that means that two observations with the same movie length but from different genres would have different predictions.

[^mixedeps]: The error term $\epsilon$ is still assumed normal with mean zero and variance $\sigma^2$. 

We also posit, and this a key point, that those group effects come from a *random* distribution. We can specify that as:

$$b_{\text{int[genre]}} \sim \text{N}(b_\text{intercept}, \sigma_\text{int\_genre})$$

This means that the random intercepts will be normally distributed and the overall intercept is just the mean of those random intercepts, and with its own variance. It also means our model will have to estimate that variance along with our residual variance. Another very common depiction is as follows, that makes explicit the addition of a random effect to the intercept:

$$\text{re}_{[\text{int\_genre}]} \sim \text{N}(0, \sigma_\text{int\_genre})$$

$$b_{\text{int[genre]}} = b_\text{intercept} +\text{re}_{[\text{int\_genre}]}$$

The same approach would apply with a random slope, where we would have a random slope for each group, and that random slope would be normally distributed with its own variance. Compared to our first formulation it would look like this:

$$b_{\text{length[genre]}} \sim \text{N}(b_\text{length}, \sigma_\text{length\_genre})$$

A simple random intercept and slope is just the start. As an example, we can let the intercepts and slopes correlate with one another, and we could have multiple grouping factors, as well as allowing multiple features and even interactions themselves to vary by group! So mixed models can get quite complex, but the basic idea is still the same: we are allowing parameters to vary across groups, and we are estimating the variance of those parameters.


### Using a mixed model  {#sec-mixed-models-using}

One of the key advantages of a mixed is that we can use it when the observations within a group are not independent. This is a very common situation in many fields, and it's a good idea to consider this when you have grouped data. As an example we'll use the happiness data for all available years, and we'll consider the country as a grouping variable. In this case, observations within a country are likely to be more similar to each other than to observations from other countries. Such **longitudinal data** is a classic example of when to use a mixed model. This is also a case where we wouldn't just throw in 'country' as a feature like any other factor, since there are `r n_distinct(df_happiness_all$country)` countries in the data. We need an easier way to handle so many groups!

In general, to use mixed models we have to specify a random effect pertaining to the categorical feature of focus, but that's the primary difference from our previous approaches used for linear or generalized linear models. For our example, we'll look at a model with a random intercept for the country feature, and one that adds a random coefficient for the yearly trend across countries. This means that we are allowing the intercepts and slopes to vary across countries, and the intercepts and slopes can correlate with one another. Furthermore, by recoding year to start at zero, the intercept will represent the happiness score at the start of the data. In addition, to see a more reasonable effect, we also divide the yearly trend by 10, so the coefficient provides the change in happiness score per decade.


:::{.panel-tabset}

##### R

We'll use the [lme4]{.pack} package in R which is the most widely used package for mixed models. 

```{r}
#| eval: false
#| label: r_mixed_model

library(lme4)

df_happiness_all = read_csv("https://tinyurl.com/worldhappinessallyears")

df_happiness_all = df_happiness_all |> 
    mutate(decade_0 = (year - min(year))/10)

# random intercepts are specified by a 1
model_ran_int = lmer(
    happiness_score ~ decade_0 + (1| country), 
    df_happiness_all
)

model_ran_slope = lmer(
    happiness_score ~ decade_0 + (1 + decade_0 | country), 
    df_happiness_all
)

# not shown
summary(model_ran_int)
summary(model_ran_slope)
```

##### Python

As with our recommendation with GAMs later, R is the better tool for mixed models, as the functionality is overwhelmingly better there for modeling and post-processing. However, you can use [statsmodels]{.pack} in Python to fit them as well[^mixedpy].

[^mixedpy]: One of your authors worked for several years with the key developer of the mixed models functionality in [statsmodels]{.pack}. As such, we can say there is zero doubt about the expertise going into its development, as there are few in the world with such knowledge. Even so, you probably won't find the functionality is not as mature or as expansive as what you get in R.


```{python}
#| eval: true
#| label: python_mixed_model
#| output: false
import statsmodels.api as sm

df_happiness_all = pd.read_csv("https://tinyurl.com/worldhappinessallyears")

df_happiness_all = (
    df_happiness_all
    .assign(decade_0 = lambda x: (x['year']- x['year'].min())/10)
)

model_ran_int = sm.MixedLM.from_formula(
    "happiness_score ~ decade_0", 
    df_happiness_all, 
    re_formula='1',
    groups=df_happiness_all["country"]  
).fit()

model_ran_slope = sm.MixedLM.from_formula(
    "happiness_score ~ decade_0", 
    df_happiness_all, 
    re_formula='1 + decade_0',
    groups=df_happiness_all["country"]  
).fit()

# not shown
# model_ran_int.summary()
# model_ran_slope.summary()
```

:::

@tbl-mixed-model shows some typical output from a mixed model, focusing on the random slope model. The fixed effect part (Fixed) is your basic GLM result and interpreted the same way. Nothing new there, and we can see a slight positive decade trend in happiness, though maybe not a strong one. But the random effects (Random) are where the action is! We can see the standard deviation of the random effects, i.e., the intercepts and slopes. We can also see the correlation between the random intercepts and random slopes. And finally we also have the residual (observation level) standard deviation, which is interpreted the same as with standard linear regression. Depending on your modeling tool, the default result may be in terms of variances and covariances rather than standard deviations and correlations, but generally similar.

\small
```{r}
#| echo: false
#| label: tbl-mixed-model
#| tbl-cap: Mixed Model Results
#| out.width: 90%

library(lme4)

df_happiness_all = read_csv("data/world_happiness_all_years.csv")
df_happiness_all_show = df_happiness_all |> 
    mutate(decade_0 = (year - min(year))/10)

model_ran_int = lmer(
    happiness_score ~ decade_0 + (1| country), 
    df_happiness_all_show
)

model_ran_slope = lmer(
    happiness_score ~ decade_0 + (1 + decade_0 | country), 
    df_happiness_all_show
)

vc = model_ran_slope |> mixedup::extract_vc(ci_level = 0)
vcor = model_ran_slope |> mixedup::extract_vc(ci_level = 0, show_cor = TRUE)
vcor = vcor$Cor$country
cor_sign = ifelse(sign(vcor[2]) > 0, 'positive', 'negative')


# gtsummary::tbl_regression(model_ran_slope, tidy_fun = broom.mixed::tidy) # doesn't include intercept and default is 'include everything'.
# broom.mixed::tidy(model_ran_slope) |>
#     mutate(
#       effect = ifelse(str_detect(effect, 'ran_pars'), 'Random', 'Fixed'),
#       term = ifelse(str_detect(term, 'Observation'), '', term),
#       term = str_remove_all(term, '\\(|\\)'),
#     ) |>
#     group_by(effect) |>
#     gt() |> 
#     sub_missing(
#       missing_text = "",
#     )
parameters::parameters(model_ran_slope) |> 
    as_tibble() |>
    select(-CI) |> 
    mutate(
        across(where(is.numeric), \(x) as.character(round(x, 2))),
        Effects = str_to_title(Effects)
    ) |>
    select(-c(t:p)) |> 
    group_by(Effects) |>
    gt() |> 
    sub_missing(
        missing_text = "",
    )
```
\normalsize

```{r}
#| echo: false
#| label: get-happiness-range
happy_range = round(range(df_happiness_all$happiness_score))
```

In this case, we can see notable variability attributable to the random effects. How do we know this? Well, if our happiness score is on a `r happy_range[1]` - `r happy_range[2]`  scale, and we have a standard deviation of happiness of `r round(sd(df_happiness_all$happiness_score), 2)` before accounting for anything else. So, we might surmise that having an effect of that size (roughly `r round(vc$sd[1], 2)`) is a relatively notable amount, as it suggests we would move around about a standard deviation of happiness just going from country to country (on average). 

The trend variability is also notable. Although the overall effect is small, the coefficient standard deviation shows that it ranges from a notable negative trend to a notable positive trend ($\pm$ `r round(vc$sd[2], 2)`) depending on which country we're looking at.

In addition, we can also see that the correlation between the random intercepts and random slopes is  `r cor_sign`, which means that the groups with higher starting points have `r ifelse(cor_sign == 'positive', 'more positive', 'more negative')` slopes. 

Now let's look at the estimates for the random effects for the model with both intercepts and slopes[^mixedup].

[^mixedup]: MC provides a package for mixed models in R called [mixedup]{.pack}. It provides a nice way to extract random effects and summarize such models ([link](https://github.com/m-clark/mixedup)).


:::{.panel-tabset}

##### R

```{r}
#| eval: true
#| label: r_mixed_model_random_effects
estimated_RE = ranef(model_ran_slope)
# mixedup::extract_random_effects(model_ran_slope) # prettier version
```

##### Python

```{python}
#| eval: true
#| label: python_mixed_model_random_effects
estimated_RE = pd.DataFrame(model_ran_slope.random_effects)
```

:::


```{r}
#| echo: false
#| eval: false
#| label: fig-random-effects
#| fig-cap: Estimated Random Effects
#| out-width: 100%
#| fig-height: 8

model_ran_int_code = lmer(
    happiness_score ~ decade_0 + (1| cntry_code), 
    df_happiness_all_show
)

model_ran_slope_code = lmer(
    happiness_score ~ decade_0 + (1 + decade_0 | cntry_code), 
    df_happiness_all_show
)

x = visibly::plot_coefficients(model_ran_slope_code, ranef = TRUE, which_ranef = 'cntry_code') 

all_re = mixedup::extract_random_effects(model_ran_slope_code)

re_usa = all_re  |> 
    filter(group == 'USA')

ran_coef_usa = mixedup::extract_random_coefficients(model_ran_slope_code)  |> 
    filter(group == 'USA')
    
ran_ef_extreme_pos =  all_re |>
    arrange(desc(value)) |> 
    group_by(effect) |> 
    slice(1:3)
ran_ef_extreme_neg =  all_re |>
    arrange(value) |> 
    group_by(effect) |> 
    slice(1:3)


p = (
    (x[[1]] + labs(x = '', y = 'Intercept')) /
    (x[[2]] + labs(x = 'Country', y = 'Decade\nTrend')) + scale_linewidth_manual(values = rep(10, 1000)) 
) &
    scale_x_discrete(guide = guide_axis(n.dodge=2)) &
    plot_layout(guides = 'collect') &
    plot_annotation(
        subtitle = 'Estimated random effects for country',
        caption = 'Each point represents the random effect for a country.\nBold points are significantly different from 0.'
    ) &
    theme(
        axis.ticks.x = element_line(linewidth = .5),
        axis.text.x = element_text(size = 4, angle = 90, hjust = 1, vjust = .5),
        panel.grid.major = element_blank()
        # axis.text.x = element_blank(),
    )

p

ggsave("img/lm-extend-random_effects.svg", width = 8, height = 8)
```

![Estimated Random Effects](img/lm-extend-random_effects.svg){#fig-random-effects width=100%}


How do we interpret these *deviations*? For starters, they are deviations from the fixed effect for the intercept and decade trend coefficient. So here that means anything negative is an intercept or slope below the corresponding fixed effect value, and anything positive is above that value. If we want the specific effect for a country, we just add its random effect to the fixed effect, and we can refer to those as **random coefficients**. 

For example, if we wanted to know the effects for the US, we would add its random effects to the population level fixed effects. This is exactly what we did in the previous section in our interpretation with the interaction model. However, you can typically get these from the package functionality directly. The result shows that the US starts at a very high happiness score, but actually has a negative trend over time[^lifeexpusa].

[^lifeexpusa]: Life expectancy in the US has actually [declined recently](https://www.hsph.harvard.edu/news/hsph-in-the-news/whats-behind-shocking-u-s-life-expectancy-decline-and-what-to-do-about-it/).

:::{.panel-tabset}

##### R

```{r}
#| eval: true
#| label: r_mixed_model_random_effects_usa

# coef(model_ran_slope) stored here
ranef_usa = estimated_RE$country
ranef_usa = ranef_usa |> 
    rownames_to_column('country') |> 
    filter(country == 'United States') 

ranef_usa[1, c('(Intercept)',  'decade_0')] + fixef(model_ran_slope) 
```

##### Python

```{python}
#| eval: true
#| label: python_mixed_model_random_effects_usa

ranef_usa = estimated_RE['United States'].rename({'Group': 'Intercept'})
ranef_usa + model_ran_slope.fe_params
```
:::


:::{.callout-note title='Averages in Mixed Models' collapse='true'}
In the linear mixed effect model setting with a random intercept, the fixed effects can be seen as the (population) average effects, but this is not exactly what you are getting from the mixed model. To make the distinction clear, consider family groups and a gender effect for males vs. females. The linear regression and some other types of models (e.g., estimated via generalized estimating equations) would give you the average effect male-female difference across all families. The mixed model actually tells you the male-female difference as if they were in the same family (e.g., siblings). Again, in the simplest mixed model setting these are the same. Beyond that that, when we start dealing with random slopes and non-gaussian distributions, they are not. 

In general, if we set the random effect to 0 to get a prediction, that tells us what the prediction would be for a typical group, in this case, a typical country. Often we want to get something more like the average slope or prediction across countries that we would have with linear regression. This gets us back to the idea of the **marginal effect** we discussed earlier. While the mechanics are not straightforward for mixed models, the tool use generally takes no additional effort.
:::

Let's plot those random coefficients together to see how they relate to each other. 

```{r}
#| echo: false
#| eval: false
#| label: fig-random-effects-cor-plot
#| fig-cap: Random Effects for Mixed Model
#| out-width: 100%
all_re_wide =  mixedup::extract_random_coefficients(model_ran_slope_code)  |> 
    select(group, effect, value) |>
    pivot_wider(names_from = effect, values_from = value)

p_re_cor = all_re_wide |>
    ggplot(aes(x = Intercept, y = `decade_0`)) +
    geom_hline(yintercept = 0, color = 'gray75', linewidth = 2) +
    geom_point(aes(), 
        data = all_re_wide |> filter(group %in% c(ran_ef_extreme_neg$group, ran_ef_extreme_pos$group)), 
        color = okabe_ito[2],
        size =3
    ) +
    geom_point(aes()) +
    geom_text_repel(
        aes(label = group, x = Intercept, y = `decade_0`), 
        segment.color = 'gray75',
        size = 2,
        max.overlaps = 20
    ) +
    labs(
        x = 'Intercept', 
        y = 'Decade\nTrend',
        subtitle = 'Estimated random coefficients for country',
        caption = 'Some more extreme country points highlighted.'
    )

p_re_cor

ggsave("img/lm-extend-random_effects_cor_plot.svg", width = 8, height = 6)
```

![Random Effects for Mixed Model](img/lm-extend-random_effects_cor_plot.svg){#fig-random-effects-cor-plot width=100%}

From this plot, we can sense why the estimated random effect correlation was negative. For individual country results, we can see that recently war-torn regions like Syria and Afghanistan have declined over time even while they started poorly as well. Some like Guinea and Togo started poorly but have improved remarkably over time.  Many western countries started high and mostly stayed that way, though generally with a slight decline. Perhaps there's only one direction to go when you're already starting off well!

:::{.callout-note title='Always scale features for mixed models' collapse='true'}

Your authors have run a lot of these models. Save yourself some trouble and standardize or otherwise scale your features before fitting the model. Just trust us, or at least, don't be surprised when your model doesn't converge.

:::


### Mixed model summary {#sec-mixed-models-summary}

For a model with just one feature, we certainly had a lot to talk about! And this is just a glimpse of what mixed models have to offer, and the approach can be even richer than what we've just seen.  But you might be asking- Why don't I just put genre country the model like other categorical features? In the case of genre for movie reviews where there are few groups, that's okay. But doing that with features with a lot of levels would typically result in estimation issues due to having so many parameters to estimate. In general mixed models provide several advantages for the data scientist:

- Any coefficient can be allowed to vary by groups, including other random effects. It actually is just an interaction in the end as far as the linear predictor and conceptual model is concerned.
- The group-specific effects are *penalized*, which shrinks them toward the overall mean, and makes this a different approach from just adding a 'mere interaction'. This helps avoid overfitting, and that penalty is related to the variance estimate of the random effect. In other words, you can think of it as running a penalized linear model where the penalty is applied to the group-specific effects (see @sec-estim-penalty).
- Unlike standard interaction approaches, we can estimate the covariance of the random effects, and specify different covariance structures for observations within groups. Standard interactions implicitly assume independence.
- Because of the way they are estimated, mixed models can account for lack of independence of observations[^indieobs], which is a common issue in many datasets. This is especially important when you have repeated measures, or when you have a hierarchical structure in your data, such as students within schools, or patients within hospitals.
- Standard modeling approaches can estimate these difficult models very efficiently, even with thousands of groups and millions of observations.
- The group effects are like a very simplified **embedding** (@sec-data-cat), where we have taken a categorical feature and turned it into a numeric one, like those shown in @fig-random-effects. This may help you understand other embedding techniques that are used in other places like deep learning.
- When you start to think about random effects and/or distributions for effects, you're already thinking like a Bayesian (@sec-unc-bayes), who is always thinking about the distributions for various effects. Mixed models are the perfect segue from standard linear model estimation to Bayesian estimation, where everything is random.
- The random effect is akin to a latent variable of 'unspecified group causes'. This is a very powerful idea that can be used in many different ways, but importantly, you might want to start thinking about how you can figure out what those 'unspecified' causes may be!
- Group effects will almost always improve your model's predictive performance relative to not having them, especially if you weren't including those groups in your model because of how many groups there were.

[^indieobs]: Independence of observations is a key assumption in linear regression models, and when it's violated, the standard errors of the coefficients are biased, which can lead to incorrect inferences. Rather than hacking a model (so-called 'fixed effects' models) or 'correcting' the standard error (e.g., with some 'sandwich' or estimator), mixed models can account for this lack of independence through the model itself.

In short, mixed models are a fun way to incorporate additional interpretive color to your model, while also getting several additional benefits to help you understand your data!




## Generalized Additive Models {#sec-gam}

<!-- :::{.content-visible when-format='html'}
> Wiggle, wiggle, wiggle, yeah!
> -- LMFAO
:::
 -->

<!-- TODO: Chelsea suggested diving a bit deeper here, and we should explain the function code a bit more. But getting into a spline demo seems a bit much. -->


### When straight lines aren't enough {#sec-gam-curve}

Linear models, as their name implies, generally assume a linear relationship between the features and the target by default. But fitting a line through your data isn't always going to be the best approach. Not every relationship is linear and not every relationship is monotonic. Sometimes, you need to be able to model a relationship that has a fair amount of nonlinearity -- they can appear as slight curves, waves, and any other type of wiggle that you can imagine.


In other words, we can go from the straight line here:

:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-regular-linear-effect
#| fig-cap: A Standard Linear Relationship
df_happiness |> 
    filter(year == 2018) |>
    ggplot(aes(healthy_life_expectancy_at_birth, happiness_score)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = 'Life Expectancy', y = 'Happiness\nScore')

ggsave("img/lm-extend-linear_effect.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![A Standard Linear Relationship](img/lm-extend-linear_effect.svg){#fig-regular-linear-effect}
:::



To the curve seen here:

:::{.content-visible when-format='html'}  
```{r}
#| echo: false
#| label: fig-gam-model-effect
#| fig-cap: A Curvilinear Relationship

# df_happiness |> 
    read_csv('https://tinyurl.com/worldhappiness2018') |> 
    ggplot(aes(healthy_life_expectancy_at_birth, happiness_score)) +
    geom_point() +
    geom_smooth(method = 'gam', se = FALSE) +
    labs(x = 'Life Expectancy', y = 'Happiness\nScore')

ggsave("img/lm-extend-nonlinear_effect.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![A Curvilinear Relationship](img/lm-extend-nonlinear_effect.svg){#fig-gam-model-effect}
:::


There are many ways to go about this, and common approaches include using polynomial regression, using feature transformations, or using a nonlinear regression model. The curved line in @fig-gam-model-effect is the result of what's called a **spline**. It is created by a feature and expanding it to multiple columns, each of which is a function of the original feature. We then fit a linear model to that data as usual. What this ultimately means is that we can use a linear model to fit a curve through the data. While this might not give us the same tidy explanation that a typical linear approach would offer, we will certainly get better predictions if it's appropriate, and a better understanding of the reality and complexity of the true relationship. It's also useful for exploratory purposes, and visualization tools can make it easy. Here is some demonstration code (result not shown).

<!-- [^ggplotly]: [Plotly]{.pack} is directly available in R and Python, and [plotnine]{.pack} is the [ggplot]{.pack} equivalent in Python. -->


:::{.panel-tabset}

##### R

```{r}
#| eval: false
#| label: r_gam_ggplot

x = rnorm(1000)
y = sin(x)

tibble(x, y) |> 
    ggplot(aes(x = x, y = y)) +
    geom_smooth(method = 'gam', se = FALSE) 
```

##### Python

```{python}
#| eval: false
#| label: python_model_gam

import plotly.graph_objects as go
import numpy as np

x = np.random.normal(size = 1000)
y = np.sin(x)

fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x = x, 
        y = y,
        line_shape = 'spline'
    )
)
```

:::

Models incorporating this type of effect belong to a broad group of *generalized additive models* (**GAM**s). When we explored interactions and mixed models, we explored how the feature-target relationship varies with another feature. There we focused on our feature and its relationship to the target at different values of other features. When we use a GAM, our initial focus is just on a specific feature and how its relationship with the target changes at different of that feature's values. How are we going to do this, you might ask? Conceptually, we will have a model that looks like this:

$$
y = f(X) + \epsilon
$$

This isn't actually any different than what you've seen - it really isn't! It's just shorthand for the input $X$ being fed into a function $f()$  of some kind, exactly as we saw in @sec-models-expression-math. That function is very flexible, and can be anything you want. 

What we're going to do with the function now is expand the feature $x$ in some way, which on a practical level means it will actually become multiple columns of input instead of just one. Some approaches to creating the expansion can be quite complex, with the goal of tackling spatial, temporal, or other aspects of the data. But, in the end, it's just extra columns in the **model matrix** that go into the model fitting function like any other feature. This enables the model to explore a more complex representation space, ultimately helping us capture non-linear patterns in our data.

At this point, you might be asking yourself, "Why couldn't I just use some type of polynomial regression or even a nonlinear regression?". Of course you could, but both have limitations relative to a GAM. If you are familiar with polynomial regression, where we add columns that are squares, cubes, etc. of the original feature, you can think of GAMs as a more general approach, and very similar in spirit. But that polynomial approach assumes a specific form of nonlinearity, and has no regularization. This means that it tends to overfit the data you currently have, and *you* are forcing curves to fit through the data, rather than exploring what naturally may arise.

Another approach besides polynomials is to use a **nonlinear regression model**. In this setting, you need to know what the underlying functional form is. An example is a [logistic growth curve model for population growth](https://en.wikipedia.org/wiki/Population_growth#Logistic_equation). Without taking extra steps, such models can also overfit. Furthermore, outside of well-known physical, chemical, or biological processes, it's rarely clear what the underlying functional form should be. At the very least, we wouldn't know a formula for life expectancy and happiness!

GAMs are better in such settings because they fit the data well without needing to know the underlying form. They also prevent overfitting in smaller data and/or more complex settings by using a penalized estimation approach. We can use them for multiple features at once, and even include interactions between features.  We also can use different types of splines, which you can think of as different functions to apply to features, to capture different types of nonlinearities. Here is another formal definition of a GAM that makes it more clear we can deal with multiple features.

$$
\hat{y} = \sum \mathbf{X_j\beta_j}
$$

In this case, each $X_j$ is a matrix of the feature and its **basis expansion**, which is created by the spline function. The $\beta_j$ are the coefficients for each of those basis expansion columns. But a specific $X_j$ could also just be a single feature and it's coefficient to model a linear relationship. You can pick and choose, which features to use it on, and how to use it.

The nice thing is that you don't have to worry about the details of the basis expansion - the package you choose will take care of that for you. You'll have different options, and often the default is fine, but sometimes you'll want to adjust the technique and how 'wiggly' you want the curve to be.


### A standard GAM {#sec-gam-standard}

Now that you have some background, let's give this a shot! In most respects, we can use the same sort of approach as we did with our other linear model examples. For our example here, we'll use the model that was depicted in @fig-gam-model-effect, which looks at the relationship between life expectancy and happiness score from the world happiness data (2018). Results are simplified in the subsequent table.

:::{.panel-tabset}

##### R

We'll use the very powerful `mgcv` package in R. The `s` function will allow us to use a spline approach to capture the nonlinearity. The `bs` argument specifies the type of spline, and here we use a B-spline, but there are many options. The `k` argument provides a means to control the complexity of the spine, historically thought of as 'knots' where bends in the curve can occur.
```{r}
#| eval: false
#| label: r_gam_demo
library(mgcv)

df_happiness_2018 = read_csv('https://tinyurl.com/worldhappiness2018')

model_gam = gam(
    happiness_score ~ s(healthy_life_expectancy_at_birth, bs = 'bs'), 
    data = df_happiness_2018
)

summary(model_gam)
```

##### Python

We can use the [statsmodels]{.pack} package in Python to fit a GAM, or alternatively, [pygam]{.pack}, and for consistency with previous models, we'll choose the former. Honestly though, you should use R's [mgcv]{.pack}, as these require notably more work without having some of the basic functionality available there. In addition, there is an ecosystem of R packages to further extend [mgcv]{.pack}'s capabilities, while these are not as mature in Python. Here we must create the spline explicitly first, a so-called 'B-spline' and add it to the model.

```{python, python_gam_demo}  
#| eval: false
#| label: python_gam_demo
from statsmodels.gam.api import GLMGam, BSplines

df_happiness_2018 = pd.read_csv('https://tinyurl.com/worldhappiness2018')

bs = BSplines(
    df_happiness_2018['healthy_life_expectancy_at_birth'],
    df = 9,
    degree = 3
)

gam_happiness = GLMGam.from_formula(
    'happiness_score ~ healthy_life_expectancy_at_birth', 
    smoother = bs,
    data = df_happiness_2018
)

model_gam = gam_happiness.fit()

model_gam.summary()
```

:::

\footnotesize
```{r}
#| echo: false
#| label: tbl-gam-model-output
#| tbl-cap: GAM Model Results
# maybe revisit if there is a better table approach in the future.
library(mgcv)

df_happiness_2018 = read_csv('data/world_happiness_2018.csv')

model_gam = gam(
    happiness_score ~ s(healthy_life_expectancy_at_birth, bs = "bs"), 
    data = df_happiness_2018
)

# library(flextable)
# looks best but not easy to manipulate at all, docs for models are non-existent, the only thing it says about signif is to remove the R global option which then proceeds to f up the footer
# this is awful to have to do and not at all generalizable, but it works for now.

d_0 = flextable::as_flextable(model_gam)$body$dataset  |>
    mutate(
        Term = ifelse(Term == '(Intercept)', 'Intercept', Term),
        Term = ifelse(Term == 's(healthy_life_expectancy_at_birth)', 's(Life Exp.)', Term),
        Component = str_remove(Component, 'A. |B. '),
        across(where(is.numeric), \(x) round(x, 2))
    ) |> 
    select(-signif)
header_smooth = c(' ',' ', 'EDF','Ref.df','F.value','p.value')
d_0[2,] = header_smooth

d_0 |> 
    gt() |> 
    tab_style(locations = cells_body(rows = 2), style = cell_text(size = 'small', color = 'gray50'))  |> # px or words but not actual font size specification; ignored by latex anyway
    tab_style(locations = cells_column_labels(), style = cell_text(color = 'gray50'))
```
\normalsize



When you look at the model output, what you get will depend a lot on the tool you use, and the details are mostly beyond the scope we want to present here (check out Michael's [-@clark_generalized_2022] document on GAMs for a lot more). But in general, the following information will be provided as part of the summary ANOVA table or as an attribute of the model object:

- **coefficients**: The coefficients for each of the features in the model. For a GAM, these are the coefficients for the basis expansion columns, as well as standard linear feature effects.

- **global test of a feature**: Some tools will provide a statistical test of the significance of the entire feature's basis expansion, as opposed to just the individual coefficients. In the table we have the intercept and the summarized smooth term.

- **edf**/**EDoF**: Effective degrees of freedom. This is a measure of 'wiggle' in the relationship between the feature and the target. The higher the value, the more wiggle you have. If you have a value close to 1, then you have a linear relationship. With our current result, we can be pretty confident that a nonlinear relationship gives a better idea about the relationship between a country's life expectancy and happiness than a linear one.

- **R-squared**: Adjusted/Pseudo $R^2$ or *deviance explained*. This is a measure of how much of the variance in the target is explained by the model. The higher the value, the better the model. Deviance explained is an analog to the unadjusted $R^2$ value for a Gaussian model that is used in the GLM setting. It's fine as a general assessment of prediction-target correspondence, and in this case, we might be feeling pretty good about the model.

Far more important than any of these is the visual interpretation, and we can get plots from GAMs easily enough.

:::{.panel-tabset}

##### R

```{r}
#| eval: false
#| label: r_gam_plot
# not shown
plot(model_gam)
```

##### Python

```{python}
#| eval: false
#| label: python_gam_plot
# not shown
model_gam.plot_partial(0, cpr=True)
```

:::


Unfortunately the default package plots are not pretty, and sadly also aren't provided in the same way we'd expect for interpretation. But they're fine for a quick look at your wiggly result. We provide a better looking one here[^gamplot]. The main interpretation is that there is not much relationship between life expectancy and happiness score until you get to about 60 years of life expectancy, and then it increases at a faster rate.  Various tools are available to easily plot the derivatives for more understanding.

[^gamplot]: We used the [see]{.pack} package in R for a quick plot. We also recommend the functionality via the [gratia]{.pack} package to visualize the derivatives, which will show where the feature effect is changing most.

:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-gam-plot
#| fig-cap: Visualizing a GAM
#| out.width: 100%
model_gam = mgcv::gam(
    happiness_score ~ s(healthy_life_expectancy_at_birth, bs = "bs"), 
    data = df_happiness
)

library(see)

plot(
    modelbased::visualisation_recipe(
        modelbased::estimate_expectation(model_gam),
        ribbon = list(alpha = 0.1),
        point  = list(color = okabe_ito[['orange']]),
        line   = list(color = okabe_ito[['darkblue']])
    ),
    alpha = 0.5
) +
    labs(
        x = 'Life Expectancy',
        y = 'Happiness\nScore',
        title = 'GAM Results'
    )

ggsave("img/lm-extend-gam_plot.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![Visualizing a GAM](img/lm-extend-gam_plot.svg){#fig-gam-plot}
:::

### GAM Summary

To summarize, we can use a GAM to model nonlinear relationships with a linear model approach. We can use splines to capture those nonlinearities, and we can use a penalized approach to control the amount of wiggle in our model. What's more, we can interact the wiggle with other categorical and numeric features to capture even more complexity in our data. This allows us to model spatial, temporal, and other types of data that have complex relationships.

GAMs are a *very* powerful modeling tool that take us a step toward more complex models, but without the need to go all the way to a neural network or similar. Plus they still provide standard statistical inference information. In short, they're a great tool for modeling!

:::{.callout-note title='GAMs are random effects models' collapse='true'}

It turns out that [GAMs have a very close relationship to mixed models](https://m-clark.github.io/generalized-additive-models/appendix.html#a-comparison-to-mixed-models), [where splines can be thought of as random effects](https://tjmahr.github.io/random-effects-penalized-splines-same-thing/) (see [Gavin Simpson's post also](https://fromthebottomoftheheap.net/2021/02/02/random-effects-in-gams/)), and GAMs can even incorporate the usual random effects for categorical variables. So you can think of GAMMs, or generalized additive mixed models, as a way to combine the best of both GAM and mixed model worlds into one extremely powerful modeling that covers a lot of ground. It's also a great baseline and sanity check model for when you're trying boosting or deep learning models, where a well-specified GAMM can be hard to beat for tabular data and still be highly interpretable, while coming with built-in uncertainty estimates and better visualization tools.

:::

## Quantile Regression {#sec-lm-extend-quantile}

<!-- :::{.content-visible when-format='html'}
> Oh, you think the median is your ally. But you merely adopted the median. I was born in it, molded by it. I didn't see the mean until I was already a data scientist, and by then it was nothing to me but misleading.
~ Bane (probably)
::: -->


People generally understand the concept of the arithmetic mean, or 'average'. You see it some time during elementary school, it gets tossed around in daily language, and it is statistically important. After all, so many distributions depend on it! Why, though, do we feel so tied to it from a regression modeling perspective? Yes, it has handy features, but  it can limit the relationships we can otherwise model effectively. Here we'll show you what to do when the mean betrays you -- and trust us, the mean will betray you at some point!


### When the mean breaks down {#sec-quantile-break}

In a perfect data world, we like to assume the mean is equal to the middle observation of the data: the *median*. But that is only when things are symmetric though, and usually our data comes loaded with challenges. Skewness and even just a few extreme scores in your data may cause a rift between the median and the mean. 

Let's say we take the integers between 1 and 10, and find the mean. 

$$\frac{1+2+3+4+5+6+7+8+9+10}{10} =  5.5$$

The middle value in that vector of numbers would also be 5.5. 

What happens if we replace the 1 with a more extreme value, like -10?

$$\frac{-10+2+3+4+5+6+7+8+9+10}{10} =  4.5$$

With just one dramatic change, our mean went down by a whole point. The median observation, though, is still 5.5. In short, the median is invariant to wild swings out in the tails of your numbers.  

You might be saying to yourself, "Why should I care about this central tendency chicanery?" Let us tell you why - the least squares approach to the standard linear model dictates that the regression line needs to be fit through the means of the variables. If you have extreme scores that influence the mean, then your regression line will also be influenced by those extreme scores.

Consider the following regression line:

:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-linear-line-no-extremes
#| fig-cap: Linear Relationship without Extreme Scores
set.seed(123)
N = 1000 
k = 2   
X = matrix(rnorm(N * k), ncol = k)  
y = -.5 + .2*X[, 1] + .1*X[, 2] + rnorm(N, sd = .5)  

dfXy = data.frame(X, y)

ggplot(dfXy, aes(x = X1, y = y)) +
    geom_point(alpha = .25) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = "x", y = "y") + 
    lims(y = c(-3, 3))

ggsave("img/lm-extend-linear_line_no_extremes.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![Linear Relationship without Extreme Scores](img/lm-extend-linear_line_no_extremes.svg){#fig-linear-line-no-extremes}
:::


Now, what would happen if we replaced a few of our observations with extreme scores?

:::{.content-visible when-format='html'}
```{r linear_line_extremes}
#| echo: false
#| label: fig-linear-line-extremes
#| fig-cap: Linear Relationship with Extreme Scores
new_df = dfXy

new_df$y[dfXy$y > quantile(dfXy$y, .95) & dfXy$X1 > quantile(dfXy$X1, .95)] = rnorm(10, 2.5, .1)

new_df$X1[dfXy$y > quantile(dfXy$y, .95) & dfXy$X1 > quantile(dfXy$X1, .95)] = rnorm(10, 2.5, .1)

new_df = new_df |> mutate(extreme = X1 != dfXy$X1 | y != dfXy$y)

ggplot(new_df, aes(x = X1, y = y)) +
    geom_point(aes(
        # color = extreme, 
        size = I(extreme + 2), 
        alpha = I(extreme + .25))
    ) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = "x", y = "y") + 
    lims(y = c(-3, 3))

ggsave("img/lm-extend-linear_line_extremes.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![Linear Relationship with Extreme Scores](img/lm-extend-linear_line_extremes.svg){#fig-linear-line-extremes}
:::


With just a casual glance, it doesn't look like our two regression lines are that different. They both look like they have a similar positive slope, so all should be good. To offer a bit more clarity, let's put those lines in the same space:

:::{.content-visible when-format='html'}
```{r}
#| echo: false
#| label: fig-both-linear-lines
#| fig-cap: Linear Relationships with and without Extreme Scores

dfXy |> 
    mutate(extreme = 'No extreme scores') |>
    bind_rows(new_df |> mutate(extreme = 'Extreme scores')) |>
    ggplot() +
    geom_smooth(
        aes(
            X1, 
            y, 
            color = extreme,
            linetype = extreme,
        ), 
        method = 'lm',
        se = FALSE
    ) +
    see::scale_color_okabeito(order = 5:6) +
    scale_linetype() +
    labs(x = "x", y = "y") +
    # lims(y = c(-3, 3)) +
    theme(
        panel.grid.major = element_line(color = alpha('black', .1)),
        # panel.grid.minor = element_line(color = alpha('black', .1)),
    )

ggsave("img/lm-extend-both_linear_lines.svg", width = 8, height = 6)
```
:::

:::{.content-visible when-format='pdf'}
![Linear Relationships with and without Extreme Scores](img/lm-extend-both_linear_lines.svg){#fig-both-linear-lines}
:::


With 1000 observations, we see that having just 10 relatively extreme scores is enough to change the regression line, even if just a little. But that little bit can mean a huge difference for predictions and the conclusions we come to.

There are a few approaches we could take here, with common approaches being dropping those observations, [winsorizing](https://en.wikipedia.org/wiki/Winsorizing) them, or doing some transformation. Throwing away data because you don't like the way it behaves is almost statistical abuse, and winsorization is just replacing those extreme values with numbers that you like a little bit better. Transformations done independently of the model rarely work for this purpose either. So let's try something else!





### A standard quantile regression {#sec-quantile-standard}

A better answer to this challenge might be to try a median-based approach instead. This is where a model like **quantile regression** becomes handy. Quantile regression is a type of regression that allows us to model the relationship between the features and the target at different quantiles of the target. For example, we can examine models at the 10th, 25th, 50th, 75th, and 90th percentiles of the target. This is very cool, as it allows us to model the relationship between the features and the target in a way that is robust to outliers and extreme scores. It's also a way to understand a type of nonlinearity that is not captured by a standard linear model, as the feature target relationship may change at different quantiles of the target.

To demonstrate this type of model, let's use our movie reviews data. Let's say that we are curious about the relationship between the `word_count` and `rating`  to keep things simple. To make it even more straightforward, we will use the standardized (scaled) version of the feature. In our default approach, we will start with a median regression.



:::{.panel-tabset}

##### R

```{r}
#| echo: true
#| label: r_quantreg
#| results: hide
library(quantreg)

model_median = rq(rating ~ word_count_sc, tau = .5, data = df_reviews)

summary(model_median)
```

```{r}
#| echo: false
#| eval: false
#| label: quantiles-vs-prediction-intervals
ul = fitted(rq(rating ~ word_count_sc, tau = .975, data = df_reviews))
ll = fitted(rq(rating ~ word_count_sc, tau = .025, data = df_reviews))
main = predict(lm(rating ~ word_count_sc, data = df_reviews), interval = 'prediction')

c(mean(ll - main[,'lwr']), mean(ul - main[,'upr']))
plot(ll, main[,'lwr'])
plot(ul, main[,'upr'])

# does it statistically work for other models? e.g., the lgbm for the mean might not have the same feature space for the q 975 https://developer.ibm.com/articles/prediction-intervals-explained-a-lightgbm-tutorial/
# https://jmarkhou.com/lgbqr/
```

##### Python

```{python} 
#| eval: false
#| label: py_quantreg
#| results: 'hide'
model_median = smf.quantreg('rating ~ word_count_sc',  data = df_reviews)
model_median = model_median.fit(q = .5)

model_median.summary()     
```

:::


```{r}
#| echo: false
#| label: tbl-quantile-model-output
#| tbl-cap: Quantile Regression Model Results

broom::tidy(model_median) |> 
    rename(
        feature = term,
        coef = estimate
    ) |> 
    select(-tau) |> 
    gt()
```

Fortunately, our interpretation of this result isn't all that different from a standard linear model - the rating should decrease by `r round(coef(model_median)['word_count_sc'],2)` for every bump in a standard deviation for the number of words, which in this case is about `r round(sd(df_reviews$word_count))` words. However, this concerns the expected *median* rating, not the mean, as would be the case with standard linear regression. 

Quantile regression is not a one-trick-pony though - being able to compute a median regression is just the default. As mentioned, we can also explore different quantiles, and this gives us the ability to answer brand new questions - does the relationship between word count and their ratings change at different quantiles of rating? Very cool!

Let's now examine the trends within 5 different quantiles of the data - .1, .3 .5, .7, and .9[^morequants].  We aren't limited to just those quantiles though, and you can examine any of them that you might find interesting.  Here is a plot of the results of these models.

[^morequants]: The R function can take a vector of quantiles, while the Python function can only take a single quantile, so you would need to loop through the quantiles.


```{r}
#| echo: false
#| eval: false
#| label: fig-quantile-lines
#| fig-cap: Quantile Regression Lines
tau_values = c(.1, .3, .5, .7, .9)

quant_values = purrr::map_df(
    tau_values, 
    ~{
        result = coef(rq(rating ~ word_count_sc, tau = .x, data = df_reviews))
        result$tau = .x
        result
    }
)

colnames(quant_values) = c('intercept', 'word_count_sc', 'tau')

quant_values$tau = as.factor(quant_values$tau)

ggplot() +
    geom_point(
        data = df_reviews, 
        mapping = aes(word_count_sc, rating), 
        position = position_jitter(width = .05),
        alpha = .1
    ) +
    geom_abline(
        aes(
            slope = word_count_sc, 
            intercept = intercept, 
            color = tau,
            linetype = tau
        ),
        linewidth = 1,
        data = quant_values
    ) +
    labs(
        x = 'Word Count (Standardized)',
        y = 'Rating',
        caption = 'Points jittered for visibility\nLines represent different quantile regression models'
    ) +
    scale_color_manual(values = unname(okabe_ito)[c(-1,-4)]) +
    theme(
        # legend.key.height = unit(1, 'cm'),
        legend.key.width = unit(.75, 'cm'),
        # legend.position = 'right'
    )

ggsave('img/lm-extend-quantile_lines.svg', width = 8, height = 6)
```

![Quantile Regression Lines](img/lm-extend-quantile_lines.svg){#fig-quantile-lines}


To interpret our visualization, we could start by saying that all of the model results suggest a negative relationship. The 10th and 90th quantiles seem weakest, while those in the middle show a notably stronger relationship. We can also see that the 90th percentile model is better able to capture those values that would otherwise be deemed as outliers using other standard techniques. The following table shows the estimated coefficients for each of the quantiles, and suggests that all word count relationships are statistically significant, since the confidence intervals do not include zero.

\small
```{r}
#| echo: false
#| label: tbl-quantile-model-output-multi-quants
#| tbl-cap: Quantile Regression Model Results
tau_values = c(.1, .3, .5, .7, .9)

rq(
    rating ~ word_count_sc, 
    tau = tau_values, 
    data = df_reviews
    ) |> 
    parameters::parameters() |>
    rename(
        feature = Parameter,
        coef = Coefficient,
        quantile = Component
    ) |>
    select(-c(CI, t:p)) |>
    gt() |>
    tab_style(
        style = cell_text(weight = "bold"),
        locations = cells_body(
            columns = vars(coef),
            rows = feature == "word_count_sc"
        )
    ) |>
    tab_footnote(
        "95% confidence intervals are shown."
    ) |>
    tab_options(
        column_labels.font.size = 12,
        footnotes.font.size = 8
    )
```

\normalsize



### The quantile loss function {#sec-quantile-loss}

Formally, the objective function for the quantile regression model can be expressed as:


:::{.content-visible when-format='html'}
![Quantile Loss Function](img/eq-quantile-reg.png){#eq-quantile-loss-function}
:::

:::{.content-visible when-format='pdf'}
$$
\text{Value} =  \Sigma \left((\pmb{\tau} - 1)\sum_{\pmb{y_{i} < \hat{y}}}(y_{i} - \hat{y}) + \pmb{\tau} \sum_{\pmb{y_{i}\geq \hat{y}}}(y_{i} - \hat{y}) \right)
$${#eq-quantile-loss-function}
:::

Compared to a standard linear regression, we are given an extra parameter for the model: $\tau$. It's a  number between 0 and 1 representing the desired quantile (e.g., 0.5 for the median)[^leastabs]. The objective function treats positive residuals differently than negative residuals. If the residual is positive, then we multiply it by $\tau$. If the residual is negative, then we multiply it by $\tau -1$. The increased penalty for over-predictions for ensures that the estimated quantile $\hat{y}$ is such that $\tau$ proportion of the data falls below it, balancing the total loss and providing a robust estimate of the desired quantile.


[^leastabs]: This is equivalent to using the least absolute deviation objective.


### Rolling our own

Given how relatively simple the objective function is, let's demystify this model further by creating our own quantile regression model and see if we can get the same results. We'll start by creating a loss function that we can use to fit our model.

:::{.panel-tabset}

##### R

```{r}
#| echo: true
#| label: quantile-loss-r
quantile_loss = function(par, X, y, tau) {
    y_hat = X %*% par
    
    residual = y - y_hat
    
    loss = ifelse(
        residual < 0, 
        (tau - 1) * residual, 
        tau * residual
    )
    
    sum(loss)
}
```

##### Python

```{python}
#| echo: true
#| label: quantile-loss-py
def quantile_loss(par, X, y, tau):
    y_hat = X.dot(par)
    
    residual = y - y_hat
    
    loss = np.where(
        residual < 0, 
        (tau-1)*residual, 
        tau*residual
    )
    
    return sum(loss)
```

:::

This code is just the embodiment of the formula shown for the quantile objective function. Compared to some of our other approaches demonstrated, we add an argument for $\tau$, but otherwise proceed very similarly. We calculate the residuals or errors in prediction, and then we sum the loss function based on those errors.


Now that we have our data and our loss function, we can fit the model almost exactly like our standard linear model. Again, note the  $\tau$ value, which we've set to .5 to represent the median.

:::{.panel-tabset}

##### R

```{r}
#| label: quantile-optim-r

X = cbind(1, df_reviews$word_count_sc)
y = df_reviews$rating


optim(
    par = c(intercept = 0, word_count_sc = 0),
    fn  = quantile_loss,
    X   = X,
    y   = y,
    tau = .5
)$par
```

##### Python

```{python}
#| label: quantile-optim-py
from scipy.optimize import minimize
import numpy as np

X = pd.DataFrame(
    {'intercept': 1, 
    'word_count_sc': df_reviews['word_count_sc']}
)
y = df_reviews['rating']

minimize(quantile_loss, x0 = np.array([0, 0]), args = (X, y, .5)).x
```

:::


Let's compare this to our previous result, and the OLS results as well. As usual, our simple code does what we need it to do! We also see that the linear regression model would produce a relatively smaller coefficient.

```{r}
#| echo: false
#| label: tbl-quantile-model-output-comparison
#| tbl-cap: Comparison of Quantile Regression Models

main_mod = coef(model_median)

our_mod = optim(
    par = c(intercept = 0, word_count_sc = 0),
    fn  = quantile_loss,
    X   = X,
    y   = y,
    tau = .5
)$par

ols_mod = coef(lm(rating ~ word_count_sc, data = df_reviews))

tibble(
    feature = c("intercept", "word_count_sc"),
    OLS = c(ols_mod[1], ols_mod[2]),
    QReg = c(main_mod[1], main_mod[2]),
    Ours = c(our_mod[1], our_mod[2])
) |> 
    gt(decimals = 3)
```


::: {.callout-note title='Yet Another Interaction' collapse='true'}
One way to interpret this result is that we have a nonlinear relationship between the word count and the rating, in the same way we had an interaction previously. In this case, our effect of word count interacts with the target! In other words, we have a different word count effect for different ratings. This is a bit of a mind bender to process, but it's another good example of how a linear approach can be used to model quirky relationships!
:::




## Wrapping Up {#sec-lm-extend-wrap}

The standard linear model is useful across many different data situations. It does, unfortunately, have some issues when data becomes a little bit more "real". When you have extreme scores or relationships that a standard model might miss, you don't necessarily need to abandon your linear model in favor of something more exotic. Instead, you might just need to think about how you are actually fitting the line through your data. 


### The common thread {#sec-lm-extend-thread}

The models discussed in this chapter are all linear models, but they add flexibility in how they model the relationship between the features and the target, and provide a non-linear aspect to the otherwise linear model. Furthermore, with tools like mixed models, GAMs, and quantile regression, we generalize our GLMs to handle even more complex data settings.

### Choose your own adventure {#sec-lm-extend-choose}


No matter how much we cover in this book, there is always more to learn. Hopefully you've got a good grip on linear models and related topics, so feel free to try out some machine learning in @sec-ml-core-concepts!



### Additional resources {#sec-lm-extend-resources}

There is no shortage of great references for mixed effects models. If you are looking for a great introduction to mixed models, we would recommend starting with yet another tutorial by one of your fearless authors! Michael Clark's *Mixed Models with R* [-@clark_mixed_2023] is a great introduction to mixed models and is [freely available](https://m-clark.github.io/mixed-models-with-R/). For a more comprehensive treatment, you can't go wrong with Gelman & Hill's *Data Analysis Using Regression and Multilevel/Hierarchical Models* [-@gelman_data_2006], and their more recent efforts in *Regression and Other Stories* [-@gelman_regression_2020], which will soon have an added component for mixed models - *Advanced Regression and Multilevel Models* [-@gelman_advanced_2024].


If you want to dive more into the GAM world, we would recommend that you start with the **Moving Beyond Linearity** chapter in *An Introduction to Statistical Learning* [@james_introduction_2021]. Not only do they have versions for both R and Python, but both have been made [available online](https://www.statlearning.com/). If you are wanting something more technical, you can't beat Simon Wood's book, *Generalized Additive Models: An Introduction with R* [-@wood_generalized_2017], or a more digestible covering of the same content by one of your own humble authors [@clark_generalized_2022]. 

For absolute depth on quantile regression, we will happily point you to the OG of quantile regression, Roger Koenker. His book, *Quantile Regression* [-@koenker_quantile_2005] is a must read for anyone wanting to dive deep into quantile regression, or just play around with his R package `quantreg`. *Galton, Edgeworth, Frisch, and Prospects for Quantile Regression in Econometrics* is another of his. And finally, you might also consider @fahrmeir_regression_2021, which takes an even more generalized view of GLMs, GAMs, mixed models, quantile regression, and more (very underrated).




## Guided Exploration {#sec-lm-extend-exercise}

These models are so much fun, you should feel comfortable just swapping any feature(s) in and out of the models. For example, for the mixed model, try using gdp per capita or life expectancy instead of (just a) trend effect. For the GAM, try using several features with nonlinear effects and see what shakes out. For quantile regression, try a different feature like movie length and use different quantiles.

:::{.panel-tabset}

##### R

```{r}
#| eval: false
#| label: r_lm_extend_exercise

library(lme4)
library(mgcv)
library(quantreg)

df_happiness = read_csv('https://tinyurl.com/worldhappinessallyears')

# Mixed Model
model_mixed = lmer(
    happiness_score ~ ??? + (1 | country), 
    data = df_happiness_2018
)

model_gam = gam(
    happiness_score ~ s(???)..., 
    data = df_happiness
)

model_quant = rq(
    happiness_score ~ ???, 
    tau = ???, 
    data = df_happiness
)
```

##### Python

```{python}
#| eval: false
#| label: python_lm_extend_exercise
import pandas as pd

import statsmodels.formula.api as smf
from statsmodels.gam.api import GLMGam, BSplines

df_happiness = pd.read_csv('https://tinyurl.com/worldhappinessallyears')

# Mixed Model
model_mixed = smf.mixedlm(
    'happiness_score ~ ???', 
    data = df_happiness, 
    groups = df_happiness['country']
).fit()

bs = BSplines(
    df_happiness['???'],
    df = 9,  # fiddle with these if you like
    degree = 3
)

gam_happiness = GLMGam.from_formula(
    'happiness_score ~ ???', 
    smoother = bs,
    data = df_happiness
).fit()

model_quant = smf.quantreg(
    'happiness_score ~ ???',
    data = df_happiness['???']
).fit(q = ???) 
```
:::