---
title: "Simple linear regression: model evaluation"
subtitle: "Notes and in-class exercises"
format: 
  html:
    embed-resources: true
    toc: true
---

# Notes

## Learning goals

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

- Use residual plots to evaluate the correctness of a model
- Explain the rationale for the R-squared metric of model strength
- Interpret the R-squared metric
- Think about ethical implications of modeling by examining the impacts of biased data, power dynamics, the role of categorization, and the role of emotion and lived experience

## Readings and videos

Choose either the reading or the videos to go through **before** class.

- Reading: Sections 1.7, 3.7, and 3.8 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)
    - Note: You do *not* need to focus on the "Ladder of Power" in Section 3.8. Transformations in general will be the focus of the *next* activity we do.
- Videos:
    - [Model evaluation: is the model wrong?](https://youtu.be/12Y4r57GuKc) ([slides](https://drive.google.com/file/d/1fkrCKhmUrorMaJEqKQBu3BFV6Y8T1A-Y/view?usp=drive_link))
    - [Model evaluation: is the model strong?](https://youtu.be/7ovhiVjXwM0) ([slides](https://drive.google.com/file/d/1QIYkfAVeSw-dFVBRJvj22zxf3cohd-SZ/view?usp=drive_link))
    - [Model evaluation: is the model fair?](https://youtu.be/40PWKLns0B4) ([slides](https://drive.google.com/file/d/1fhVyieXzUf4zqRj-0H5gu4ubJkIkl8RX/view?usp=drive_link))
    - [R Code for Evaluating and Using a Linear Model](https://voicethread.com/share/14931045/)


## Independence Assumption

One way to think about model evaluation is to consider whether or not underlying *assumptions* of our regression models are being met (or not). Asking ourselves if our models are "wrong", "strong", and "fair" approaches this from one perspective. Another perspective, that is more classical, is to consider the following four *assumptions* of linear regression:

1. Linearity
2. Independence
3. Normality
4. Equal Variance

Note that they spell "LINE" (how convenient!).

By *assumptions*, we mean that the above four "things" are needed **mathematically** in order for linear regression to "work". We won't get into the details of most of these (because it turns out that two of the above assumptions aren't really needed in most cases... ask your professor for more details), but one important one to note that was not covered in the videos/course notes is the Independence assumption.

What we mean by independence, is that the *residuals* in our model *do not depend* on one another. This may seem like an unsatisfying definition, so here are some examples:

-   Suppose I want to understand the association between a person's high school GPA and their college GPA. I collect data from *every* graduating senior, at three different high schools. If I have college GPA as my outcome, and high school GPA as my predictor, are my residuals independent? Probably not! It is reasonable to believe that students *from the same high school* may have similar GPAs, due to resources their high school may have had available, or specific teachers grading differently at one school or another. This is an example of *clustering*, where we have clusters of students within schools. The independence assumption of our linear regression model would be violated. One way to address this would be to include which high school they went to as an additional covariate in our regression model (we'll get to this with multiple linear regression), and more advanced methods are covered in a course on Correlated Data.

-   Suppose I want to understand the association between a mouse's weight and their water consumption across time. I collect data for 365 days for ten different mice, recording their weight and water consumption each day of the year. If I have weight as my predictor and water consumption as my outcome, are my residuals independent? Nope! This is an example of correlated data that is *longitudinal* in nature: I have *multiple observations* per individual (mouse) across time. A mouse's weight one day is certainly not independent of it's weight the following day. The independence assumption of our linear regression model would again be violated. One way to address this would be to include "Mouse ID" as a predictor in our regression model (again, we'll get to this with multiple linear regression).

All types of data that will violate the independence assumption of linear regression will have some sort of **correlation structure** (within individual, across time, across space, etc.). Think about clusters. If your observations fall neatly into specific clusters, your data may violate the independence assumption of linear regression.


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

# Exercises

We'll start by revisiting the Capital Bikeshare data from a previous activity. This data contains information on daily ridership for a bike sharing service in Washington DC alongside weather and calendar information. Read in the data and load required packages:

```{r warning=FALSE, message=FALSE}
# Load packages and import data
library(readr)
library(ggplot2)
library(dplyr)

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

## Exercise 1: Is the model correct?

We previously explored modeling daily ridership among registered users as a function of temperature. Create a plot of this relationship with both a curved and linear trend line. Based on this plot, do you think a linear model is correct?
  
```{r}
# Plot temp_feel vs riders_registered with a model trend
___(___, aes(x = ___, y = ___)) + 
    geom___() + 
    geom___(se = FALSE) +
    geom___(method = "lm", se = FALSE, color = "red")
```

> Your response here
   
## Exercise 2: Changing the model

In this course, we have and will continue to build "linear" regression models. "Linear" means we have a *linear combination of predictors*. It does *not* mean that the models themselves must look linear! It's possible to include "transformations" in our models in order to better match the trend. Below we create a squared version of temperature and visualize the predictions from 2 models: (1) with just temperature as a linear term and (2) with both temperature and squared temperature. (Don't worry about the new syntax in `geom_smooth()`.)

How does the correctness of the two models compare?

```{r}
bikes <- bikes %>% 
    mutate(temp_feel_squared = temp_feel^2)

# Plot the model WITHOUT squared temperature
ggplot(bikes, aes(x = temp_feel, y = riders_registered)) + 
    geom_point() + 
    geom_smooth(se = FALSE) +
    geom_smooth(method = "lm", se = FALSE, color = "red")
# Plot the model WITH squared temperature
ggplot(bikes, aes(x = temp_feel, y = riders_registered)) + 
    geom_point() + 
    geom_smooth(se = FALSE) +
    geom_smooth(method = "lm", formula = y ~ x + I(x^2), se = FALSE, color = "red")
```

> Your response here

## Exercise 3: Residual plots

Plotting the **residuals** vs the **predictions** (also called "fitted values") for each case can help us assess how wrong our model is. This will be a particularly important tool when evaluating models with multiple predictors. Construct the residual plots for our two temperature models. Do these suggest that `bike_mod1` is wrong? What about `bike_mod2`? Explain.
    
Note: Information about the residuals (`.resid`) and predictions (`.fitted`) are stored within our *model*, thus we start our `ggplot()` with the model name as opposed to the raw dataset. We will *rarely* start `ggplot()` with a model instead of the data.

```{r}
# Fit a linear model
bike_mod1 <- lm(riders_registered ~ temp_feel, data = bikes)
# Fit a quadratic model
bike_mod2 <- lm(riders_registered ~ temp_feel + temp_feel_squared, data = bikes)

# Check out the residual plot for bike_mod1 (the incorrect model)
ggplot(bike_mod1, aes(x = .fitted, y = .resid)) + 
    geom_point() + 
    geom_hline(yintercept = 0) +
    geom_smooth(se = FALSE)

# Construct the residual plot for bike_mod2 (the good model)

```

> Your response here

## Exercise 4: Another example of an incorrect model

The [Bechdel test](https://en.wikipedia.org/wiki/Bechdel_test) is a test applied to films to assess the quality of women's presence in the film. We'll be looking at movies that have had the Bechdel test applied available in the `fivethirtyeight` R package. In the **Console**, enter `install.packages("fivethirtyeight")` to install this package.

We'll be using the `bechdel` data. You can view the codebook by entering `?bechdel` in the Console.

Let's examine the relationship between international earnings (`intgross`) and movie budget (`budget`) for **films made in 1997**.

```{r, warning = FALSE, message = FALSE}
# Import the data
library(fivethirtyeight)
data(bechdel)

# Get only 1997 movies
movies_1997 <- bechdel %>% 
    filter(year == 1997)

# Construct the model
bechdel_model <- lm(intgross ~ budget, movies_1997)
```


a. Construct two plots:

```{r}
# Scatterplot of earnings and budget with linear and curved trend lines


# Residual plot for bechdel_model

```

b. These two plots confirm that our model is wrong. *What* is wrong and how might we fix it?

> Your response here

c. Identify which movie is causing the problem. Hint: `filter()` according to `budget`. Also, note that we could, but won't, take that film out of the data set and re-build the model.
    
```{r}
# NOTE: Many numbers could go in the ___
movies_1997 %>%
    filter(budget > ___)
```

> Your response here

## Exercise 5: Is the model strong? Developing R-squared intuition

The R-squared metric is a way to quantify the strength of a model. It measures how much variation in the outcome/response variable can be explained by the variation in the predictors.

Strong models have residuals that don't deviate far from 0.  So the smaller the variance in the residuals (thus larger the variance in the predictions), the stronger the model. Take a look at the picture below and write a few sentences addressing the following:

- The two rows of plots show a stronger and a weaker model. Just by looking at the blue trend line and the dispersion of the points about the line, which row corresponds to the stronger model? How can you tell? Which row would you expect to have a higher correlation?

> Your response here

- What is different about the variance of the residuals from the first to the second row?

> Your response here

![](https://mac-stat.github.io/images/155/rsquared.png)


Putting this together, the R-squared compares the *overall* variance in outcome (SSTO) to the variance in the residuals (SSE):

::: {.callout-note collapse="true"}
## R-squared

$$
R^2 = 1 - \frac{SSE}{SSTO} = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2}
$$
where $y_i$ are our observed outcomes, $i = 1, \dots, n$, $\hat{y}_i$ are our fitted values/predictions, and $\bar{y}$ is our observed *average* outcome.
:::

## Exercise 6: R-squared Interpretations

Recall `bikemod1` from Exercise 3, where we predicted registered riders by what the temperature felt like on a given day. Use the `summary` function to look out the model output for `bikemod1`, and **interpret** the $R^2$ value for this model, in the context of the problem. (NOTE: $R^2$ is reported in output here as "Multiple R-squared"). 

```{r eval = FALSE}
# Get R-squared
summary(bike_mod1)
```

> Your response here

## Exercise 7: Further exploring R-squared

In this exercise, we'll look at data from a synthetic dataset called Anscombe's quartet. Load the data in as follows, and look at the first few rows:

```{r, warning = FALSE, message = FALSE}
data(anscombe)

# Look at the first few rows

```

> Pause: Did you add code to the look at the first few rows, in the code chunk above?

The `anscombe` data is actually 4 datasets in one: `x1` and `y1` go together, and so forth. Examine the coefficient estimates (in the "Estimate" column of the "Coefficients:" part) and the "Multiple R-squared" value on the second to last line. What do you notice? How do these models compare?

```{r, eval = FALSE}
anscombe_mod1 <- lm(y1 ~ x1, data = anscombe)
anscombe_mod2 <- lm(y2 ~ x2, data = anscombe)
anscombe_mod3 <- lm(y3 ~ x3, data = anscombe)
anscombe_mod4 <- lm(y4 ~ x4, data = anscombe)

summary(anscombe_mod1)
summary(anscombe_mod2)
summary(anscombe_mod3)
summary(anscombe_mod4)
```

> Your response here

Now take a look at the following scatterplots of the 4 pairs of variables. What do you notice? What takeaway can we draw from this exercise?

```{r eval = FALSE}
ggplot(anscombe, aes(x = x1, y = y1)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x2, y = y2)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x3, y = y3)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)

ggplot(anscombe, aes(x = x4, y = y4)) +
    geom_point() + 
    geom_smooth(method = "lm", color = "red", se = FALSE)
```

> Your response here

## Exercise 8: Biased data, biased results: example 1

DATA ARE NOT NEUTRAL. Data can reflect personal biases, institutional biases, power dynamics, societal biases, the limits of our knowledge, and so on. In turn, biased data can lead to biased analyses. Consider an example.    

- Do a Google image search for "statistics professor." What do you observe?

> Your response here

- These search results are produced by a search algorithm / model. Explain why the data used by this model are not neutral.

> Your response here

- What are the potential implications, personal or societal, of the search results produced from this biased data?

> Your response here



## Exercise 9: Biased data, biased results: example 2

Consider the example of a large company that developed a model / algorithm to review the résumés of applicants for software developer & other tech positions. The model then gave each applicant a score indicating their hireability or potential for success at the company. You can think of this model as something like:            

$$\text{potential for success } = \beta_0 + \beta_1 (\text{features from the résumé})$$

Skim this [Reuter's article](https://www.reuters.com/article/us-amazon-com-jobs-automation-insight/amazon-scraps-secret-ai-recruiting-tool-that-showed-bias-against-women-idUSKCN1MK08G) about the company's résumé model.    

- Explain why the data used by this model are not neutral.

> Your response here

- What are the potential implications, personal or societal, of the results produced from this biased data?

> Your response here



## Exercise 10: Rigid data collection systems

When working with categorical variables, we've seen that our units of observation fall into neat groups. Reality isn't so discrete. For example, check out questions 6 and 9 on [page 2 of the 2020 US Census](https://www2.census.gov/programs-surveys/decennial/2020/technical-documentation/questionnaires-and-instructions/questionnaires/2020-informational-questionnaire-english_DI-Q1.pdf). With your group, discuss the following:    

- What are a couple of issues you see with these questions?

> Your response here

- What impact might this type of data collection have on a subsequent *analysis* of the census responses and the policies it might inform?

> Your response here

- Can you think of a better way to write these questions while still preserving the privacy of respondents?

> Your response here

**FOR A DEEPER DISCUSSION:** Read [Chapter 4 of Data Feminism](https://data-feminism.mitpress.mit.edu/pub/h1w0nbqp/release/3) on "What gets counted counts". 



## Exercise 11: Presenting data: "Elevating emotion and embodiment"

Note: The following example highlights work done by W.E.B. Du Bois in the late 1800s / early 1900s. His work uses language common to that time period and addresses the topic of slavery.

The types of visualizations we've been learning in this course are standard practice, hence widely understood. Yet these standard visualizations can also suppress the lived experiences of people represented in the data, hence can miss the larger point. W.E.B. Du Bois (1868--1963), a "sociologist, socialist, historian, civil rights activist, Pan-Africanist, author, writer, and editor"^[https://en.wikipedia.org/wiki/W._E._B._Du_Bois], was a pioneer in elevating emotion and embodiment in data visualization. For the Paris World Fair of 1900, Du Bois and his team of students from Atlanta University presented 60 data visualizations of the Black experience in America, less than 50 years after the abolishment of slavery. To this end, Du Bois noted that *"I wanted to set down its aim and method in some outstanding way which would bring my work to notice by the thinking world."* That is, he wanted to increase the impact of his work by partnering technical visualizations with design that better connects to lived experiences. Check out:
    
- An [article by Allen Hillery](https://dataliteracy.com/web-du-bois-story-of-resilience/) (@AlDatavizguy).

- A [complete set of the data visualizations provided by Anthony Starks](https://speakerdeck.com/ajstarks/du-bois-visualizations-originals) (@ajstarks).

Discuss your observations. In what ways do you think the W.E.B. Du Bois visualizations might have been more effective at sharing his work than, say, plainer bar charts?

> Your response here

**FOR A DEEPER DISCUSSION AND MORE MODERN EXAMPLES:** Read [Chapter 3 of Data Feminism](https://data-feminism.mitpress.mit.edu/pub/5evfe9yd/release/5) on the principle of elevating emotion and embodiment, i.e. the value of "multiple forms of knowledge, including the knowledge that comes from people as living, feeling bodies in the world."



## Reflection

What has stuck with you most in our exploration of model evaluation? Why?

> **Response:** Put your response here.



## Render your work

- Click the "Render" button in the menu bar for this pane (blue arrow pointing right). This will create an HTML file containing all of the directions, code, and responses from this activity. A preview of the HTML will appear in the browser.
- Scroll through and inspect the document to check that your work translated to the HTML format correctly.
- Close the browser tab.
- Go to the "Background Jobs" pane in RStudio and click the Stop button to end the rendering process.
- Navigate to your "Activities" subfolder within your "STAT155" folder and locate the HTML file. You can open it again in your browser to double check.

