---
title: "Model evaluation (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')
```



::: {.callout-note title = "Learning goals"}

- Is it wrong?    
    Use residual plots to evaluate the correctness of a model
- Is it strong?   
    - Explain the rationale for the R-squared metric of model strength
    - Interpret the R-squared metric
- Is it fair?    
    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


![](https://ajohns24.github.io/155_Fall_2025/images/concept_map_4.png)

:::



\




::: {.callout-note title="Additional resources"}

**Required videos**

- [Model evaluation: is the model wrong?](https://youtu.be/uUf4XYF4x3Y) ([slides](https://drive.google.com/file/d/1H6DYqemLNp2x6YiQQ606w9FiO-bkgRpd/view?usp=drive_link))
- [Model evaluation: is the model strong?](https://youtu.be/ew9eCkcUWGw) ([slides](https://drive.google.com/file/d/1Pd84YX20P1X5xfM34ZvWsYFZoljwK_pK/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))


**Optional**

- 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.


:::


    

\
\
\
\


# References {-}

Just as we should always question the safety of a foraged mushroom, and just as medications come with warning labels, we should always evaluate the quality of our model before trusting, using, or reporting it.
To this end, there are 3 important questions to ask:


\

## Is it correct? (Is it wrong?) {-}


There are four assumptions. These spell out **LINE**:   

- Linearity   
    Residuals are balanced above & below 0 across the entire range. (The model accurately describes the "shape" of the relationship between Y and X.)

- Independence    
    Residuals are unrelated / independent of one another.

- Normality   
    Residuals are *normally* distributed or clustered around 0 across the entire range.

- Equal variance    
    The variability in residuals (thus how far observed outcomes fall from the model) is roughly equal across the entire range.




::: {.callout-note title="OPTIONAL Math Box: Model assumptions"}

In practice, the **LINE** assumptions for a (Normal) simple linear regression model are mathematically presented as follows:

$$Y_i = E[Y | X_i] + \varepsilon_i = \beta_0 + \beta_1 X_i + \varepsilon_i
\;\;\;\;
\text{ where } \;\;\;\;
\varepsilon_i \stackrel{ind}{\sim} N(0, \sigma^2)$$

Equivalently, 

$$Y_i \stackrel{ind}{\sim} N(\beta_0 + \beta_1 X_i, \sigma^2)$$
:::




\
\


## Is it strong? {-}

Calculate R-squared ($R^2$).    

- $0 \le R^2 \le 1$    

- R^2^ measures the proportion of the variability in the observed response values that's *explained* by the model. Thus 1 - R^2^ is the proportion of the variability in the observed response values that's left *UNexplained* by the model.

- The closer R^2^ is to 1, the better the model.


::: {.callout-note title="OPTIONAL Math Box: Calculating R-squared"}

Let $y = (y_1, y_2, ..., y_n)$ denote our sample of observed $Y$ values, 
$\hat{y} = (\hat{y}_1, \hat{y}_2, ..., \hat{y}_n)$ denote the corresponding model predictions,
and
$y - \hat{y} = (y_1-\hat{y}_1, y_2-\hat{y}_2, ..., y_n-\hat{y}_n)$ denote the corresponding model residuals.
Then we can *partition* the variance in the observed outcomes:

$$Var(y) = Var(\hat{y}) + Var(y - \hat{y})$$

It follows that

$$R^2 = \frac{Var(\hat{y})}{Var(y)} = \frac{\sum (\hat{y}_i - \bar{y})^2}{\sum (y_i - \bar{y})^2}$$

Equivalently,

$$R^2 = 1 - \frac{Var(y - \hat{y})}{Var(y)} = \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2}$$


:::



\
\


## Is it fair? {-}
    
There's no mathematical measure or assumption here!
Rather, we must evaluate the potential ethical implications of data collection and modeling.


\
\
\
\



# Warm-up {-}

Let's explore the `high_peaks` data, which includes information on hiking trails in the 46 "high peaks" in the Adirondack mountains of northern New York state.^[original data source = `HighPeaks` data in the `Stat2Data` package; original image source = https://commons.wikimedia.org/wiki/File:New_York_state_geographic_map-en.svg]


![](https://upload.wikimedia.org/wikipedia/commons/thumb/3/35/New_York_state_geographic_map-en.svg/1600px-New_York_state_geographic_map-en.svg.png)




\
\


**EXAMPLE 1: Modeling review**

Let's begin by modeling the `time` (in hours) that it takes to complete each hike by the hike's vertical `ascent` (in feet).

```{r}
# Load packages and data
library(tidyverse)
peaks <- read_csv("https://mac-stat.github.io/data/high_peaks.csv")
```

```{r}
# Check out the data


```

```{r}
# Fit a model
peaks_model <- ___(___, ___)

# Check out a model summary table
# (the big one, not the simplified version!)

```

**Follow-up**

Write out the model formula.


\
\
\
\


::: {.callout-note title = "Interpreting model coefficients"}

When interpreting model coefficients, remember to:

1. interpret in context
2. use non-causal language
3. include units
4. talk about averages rather than individual cases

:::

\
\
\
\


**EXAMPLE 2: More modeling review**


Interpret the 2 model coefficients.


\
\
\
\


**EXAMPLE 3: Is it correct?**

Let's *evaluate* the quality of this model.
First, is it correct (i.e. not wrong)?
Comment on the **LINE** assumptions, using the below **residual plot** where relevant.

```{r}
# Residual plot
# NOTE: We start our code with the MODEL not the DATA!!!
# That's where the residuals (`.resid`) & predictions (`.fitted`) are stored
# We'll RARELY do this.
peaks_model %>% 
  ___(___(___ = .fitted, ___ = .resid)) + 
  ___() + 
  geom_hline(yintercept = 0)
```






\
\
\
\



**EXAMPLE 4: Plot the data**

Let's go back and plot the data, specifically the relationship of `time` with `ascent` (which we should always do before building the model!).
Include a representation of the linear regression model.

```{r}
# NOTE: "lm" stands for "linear model"
___ %>% 
  ___(aes(y = ___, x = ___)) +
  geom___() + 
  geom___(method = "lm", se = FALSE)
```



\
\
\
\



**EXAMPLE 5: Is it strong?**

We can measure the *strength* of this model using R-squared ($R^2$).
Report and interpret the R-squared value.
HINT: No new R code needed! Just check out the model summary table.

```{r}
summary(peaks_model)
```


\
\
\
\


**EXAMPLE 6: Is the model fair?**




\
\
\
\


# Exercises {-}


You'll work with a variety of datasets & scenarios in these exercises.
The goal is to explore the different ways that a model can be bad (or good!), not to work through 1 cohesive analysis.

\
\


**DIRECTIONS**

- Work with your group:   
    Be kind to yourself and others. Engage yourself and engage others.
    
- When you get stuck:   
    Re-read the problem, discuss with your group, ask the instructor, check solutions in the online manual.
    
    

\
\


## Exercise 1: Is the model wrong? {-}

Let's check out some weather data from Hobart, Australia, including the `maxtemp` (maximum temperature in degrees Celsius) and `day_of_year` (1--365):

```{r}
# Import data
# Take note of the filter code!
weather <- read_csv("https://mac-stat.github.io/data/weather_3_locations.csv") %>% 
  filter(location == "Hobart") %>% 
  mutate(day_of_year = yday(date)) %>% 
  select(maxtemp, day_of_year)

# Check it out
head(weather)
```

Let's use this data to model `maxtemp` by `day_of_year`: 

```{r}
# Fit a linear model of 
weather_model <- lm(maxtemp ~ day_of_year, data = weather)

# Check it out
summary(weather_model)
```

### Part a {-}

The estimated model formula is

E[`maxtemp` | `day_of_year`] = 20.105 - 0.012 `day_of_year`

Discuss the implications of the slope coefficient -0.012 in this context.
Does this make sense (what do you think is the actual relationship between these 2 variables)?



### Part b {-}

Now plot the relationship of `maxtemp` with `day_of_year` with both a curved and linear trend line.
This plot indicates that the simple linear regression model (blue) is wrong.
Specify which LINE assumption it violates.
  
```{r}
# Plot maxtemp vs day_of_year with...
# a linear regression model trend AND
# a smooth red curve (not necessarily linear)
___ %>% 
  ___(aes(y = ___, x = ___)) + 
  geom___() +
  geom___(method = "lm", se = FALSE) + 
  geom___(se = FALSE, color = "red") 
```


### Part c {-}

Check out the residual plot for this model.
Convince yourself that you could come to the same conclusion from Part a using this plot alone.

```{r}
# Check out the residual plot for weather_model
___ %>% 
  ggplot(aes(x = .fitted, y = .resid)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE)
```



\
\
\
\




## Exercise 2: What's incorrect about this model? {-}

Consider another example.
The `mammals` data includes data on the average brain weight (g) and body weight (kg) for a variety of mammals:

```{r}
# Import the data
mammals <- read_csv("https://mac-stat.github.io/data/mammals.csv")

# Check it out
head(mammals)
```

Fit a model of `brain` vs `body` weight:

```{r}
# Construct the model
mammal_model <- lm(brain ~ body, mammals)

# Check it out
summary(mammal_model)
```

### Part a {-}

Construct two plots that will help us evaluate `mammal_model`:

```{r}
# Scatterplot of brain weight (y) vs body weight (x)
# Include a model trend line (i.e. a representation of mammal_model)

```

```{r}
# Residual plot for mammal_model

```

### Part b {-}

These two plots confirm that our model is wrong.
*What* is wrong?
Specifically, beyond the relationship maybe not being linear, which of the *other* LINE assumptions is violated?
(NOTE: We again can *fix* this model by "transforming" one or both of the `brain` and `body` variables. We'll discuss this idea in our next class.)


### Part c {-}

Just for fun, let's dig into the mammals.
Discuss what you observe.

```{r}
# Label the points by the animal name!
# Discuss: What 2 things are new in this code?
mammals %>% 
  ggplot(aes(x = body, y = brain, label = animal)) + 
  geom_text() + 
  geom_smooth(method = "lm", se = FALSE) 
```

```{r}
# Zoom in
# Discuss: What 1 thing is new in this code?
mammals %>% 
  ggplot(aes(x = body, y = brain, label = animal)) + 
  geom_text() + 
  lims(y = c(0, 1500), x = c(0, 600))
```

```{r}
# Zoom in more
mammals %>% 
  ggplot(aes(x = body, y = brain, label = animal)) + 
  geom_text() + 
  lims(y = c(0, 500), x = c(0, 200))
```




\
\
\
\


## Exercise 3: Is the model strong? {-}

Let's revisit the `bikes` data, but focus just on *MONDAY* ridership:

```{r}
# Take a peak at this code for review!
mondays <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% 
  rename(rides = riders_registered) %>% 
  filter(day_of_week == "Mon")

# Check it out
head(mondays)
```

Check out the plot of `rides` vs `temp_feel`:

```{r}
mondays %>% 
  ggplot(aes(y = rides, x = temp_feel)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```

And the plot of `rides` vs `windspeed`:

```{r}
mondays %>% 
  ggplot(aes(y = rides, x = windspeed)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
```


### Part a {-}

Obtain, interpret, and compare the R-squared values for the 2 models of interest.

```{r}
# Fit the models
bike_model_1 <- lm(rides ~ temp_feel, data = mondays)
bike_model_2 <- lm(rides ~ windspeed, data = mondays)

# Obtain R-squared
___(bike_model_1)
___(bike_model_2)
```

### Part b {-}

Since the relationships modeled by `bike_model_1` and `bike_model_2` both include a single quantitative predictor, we can calculate their correlations.
Do that below *and* demonstrate how these correlations can be used to calculate the R-squared values from Part a.

```{r}
# Calculate correlation
mondays %>% 
  ___(___(rides, temp_feel), ___(rides, windspeed))

# Use the correlations to calculate the model R-squared values
# (within rounding)

```



\
\
\
\



## Exercise 4: Calculating R-squared (Part 1) {-}

Let's dig into where the R-squared metric comes from!

### Part a {-}

Check out the `model_1_info` table below that, for each day in the dataset, includes info on the *observed* ridership, the *predicted* ridership using `bike_model_1`, and the *residual*.
(You made similar calculations in the previous activity.)
*Discuss with each other how these 3 columns are related!*
For example, how can we get the `rides` column from the `predicted` and `residual` columns?

```{r}
# Don't worry about this code!
model_1_info <- mondays %>% 
  select(rides) %>% 
  mutate(predicted = bike_model_1$fitted.values, residual = bike_model_1$residuals)

# Check it out
head(model_1_info)
```


### Part b {-}

In Part a we observed that the observed `rides`, `predicted` values, and `residual` values all vary from day to day.
Here, calculate the *variance* of the observed `rides`, `predicted` values, and `residual` values.
*Identify and discuss with each other how these 3 values are related!*
For example, how can we calculate `var(rides)` using `var(predicted)` and `var(residual)`?

```{r}
model_1_info %>% 
  ___(var(rides), var(predicted), var(residual))
```
        
        
### Part c {-}

Recall that `bike_model_1` had an R-squared of 0.3137.
Challenge: How can you calculate this using the `var(rides)` and `var(predicted)` values from Part b?
Hint: Remember that R-squared is the proportion of the variability in a response variable Y that is explained by its model (thus its predictions).

```{r}
# Calculate R-squared
```

     



\
\
\
\


## Exercise 5: Calculating R-squared (Part 2) {-}

<!-- Recommended: Check out the rendered version of this exercise.  -->

Your work for `bike_model_1` is summarized below, along with additional calculations from `bike_model_2`.  
    
\

Model           Var(observed)  Var(predicted)  Var(residuals) R-squared
-------------- -------------- --------------- --------------- ----------
`bike_model_1`        2262666          709695         1552970 0.31
`bike_model_2`        2262666           65864         2196801 0.02

\

These results reflect 2 important facts:

1. We can **partition the variance** of the observed response values into the variability that's explained by the model (the variance of the predictions) and the variability that's left unexplained by the model (the variance of the residuals):   
    
$$\text{Var(observed) = Var(predicted) + Var(residuals)}$$

2. It follows that R-squared, the *proportion* of the variance in the observed response values that's explained the the model, can be calculated by:

$$R^2 = \frac{\text{Var(predicted)}}{\text{Var(observed)}} = 1 - \frac{\text{Var(residuals)}}{\text{Var(observed)}}$$


\

These properties are also reflected in the pictures below:

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



\
\


**FOLLOW-UP QUESTIONS:**

Reflecting on the above table of variances, R-squared formula, and plot...

- In terms of obtaining higher R-squared / stronger model, do we want the variance of the residuals to be close to 0 or close to the variance of the observed Y values? Why does this make sense? 

- In terms of obtaining higher R-squared / stronger model, do we want the variance of the prediction to be close to 0 or close to the variance of the observed Y values? Why does this make sense? 





\
\
\
\



## Exercise 6: Biased data, biased results {-}

In the above exercises, we focused on exploring our first 2 model evaluation questions: Is it wrong? Is it strong?
In the next 3 exercises, let's explore the third question: Is it fair?

**Data are not neutral**.
Data can reflect personal biases, institutional biases, power dynamics, societal biases, the limits of our knowledge, etc.
And biased data can lead to biased analyses.
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 hire-ability or potential for success at the company.
You can think of this model as something like:            

$$E[\text{potential | résumé features}] = \beta_0 + \beta_1 (\text{résumé features})$$

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.

- What are the potential implications, personal or societal, of the results produced from this biased data? Hence why is this model *not fair*?




\
\
\
\



## Exercise 7: Rigid data collection systems {-}

When working with categorical variables, we've observed 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 observe with these questions?

- What impact might this type of data collection have on a subsequent *analysis* of the census responses and the policies it might inform? That is, why might a resulting model *not be fair*?

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

**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 8: Data drill {-}

### Part a {-}

Let's end by practicing some data wrangling with our `peaks` data.
In addition to some basic functions (e.g. `head()`), use the `tidyverse` functions we've been accumulating: `select()`, `summarize()`, `filter()`.

```{r}
# How many hikes are in the data set?


# Show just the peak's name and its elevation for the first 6 hikes
___ %>% 
  ___ %>% 
  head()

# Calculate the average hiking time


# Show data for the hikes that are *at least* 17 miles long


# Calculate the average hiking time for hikes that are at least 17 miles long
# HINT: You'll need to use 2 tidyverse functions!
___ %>% 
  ___ %>% 
  ___

# Calculate the average hiking time for hikes with an "easy" rating
# HINT: You'll need to use 2 tidyverse functions again!


```

### Part b {-}

Let's practice some visualization!

```{r}
# Construct a visualization of hike rating


# Construct a visualization of hiking time


# Construct a visualization of the relationship between hiking time and length (distance)

```




\
\
\
\




# Reflection {-}

**Learning goals**

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


**Code**

- What new code did you learn today?

- Remember that you're highly encouraged to start tracking and organizing new code in a cheat sheet (eg: a Google doc).


## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


