--- title: "Model Evaluation (Notes)" subtitle: "Stat 253" author: "Your Name" format: html: toc: true toc-depth: 2 embed-resources: true --- # Learning Goals {.unnumbered .smaller} - Practice building and evaluating linear regression models in RStudio, including: - Create and interpret residuals vs. fitted, residuals vs. predictor plots to identify improvements in modeling and address ethical concerns. - Calculate and interpret MSE, RMSE, MAE, and R-squared in a contextually meaningful way. - Review some STAT 155 concepts and RStudio tools. - Start using the `tidymodels` package to build and evaluate regression models. - Get used to the course workflow, e.g. taking notes. \ # Notes: Quarto {.smaller .unnumbered} If you've taken STAT 155 or COMP/STAT 112 recently, you're already familiar with Quarto. If not, here are a few highlights: - Replacement for R Markdown - Document structure is nearly identical (YAML header, code chunks, etc.) - File names end in `.qmd` instead of `.rmd` - `Render` instead of `Knit` - Structure of header is slightly different - More resources available in the *R Resources* Appendix Check out the QMD for today's class! Questions? Comments? \ # Small Group Discussion {-} ## Instructions {.unnumbered .smaller} - In small groups, please first introduce yourself (in whatever way you feel appropriate) and check in with each other as human beings. - When everyone is ready, glance through the summary of concepts covered in the video (see "Video Recap" below) and discuss the following prompts: - What vocabulary or notation was new to you? - What concepts were new to you? - What concepts are still unclear to you at this moment? - Prepare to share a few highlights from your group discussion with the entire class ## Video Recap {.unnumbered .smaller} ![](https://mac-stat.github.io/STAT253/images/MLdiagram1.jpg){width=90%} We are in the **regression** setting. We want to build a model of some **quantitative** output variable $y$ by some predictors $x$: $$y = f(x) + \varepsilon$$ There are many regression tools that we might use to build this model. We'll use a **linear regression** model which assumes that $y$ is a linear combination of the $x$'s: $$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots \beta_p x_p + \varepsilon$$ After *building* any model, it's important to *evaluate* it: **Is our _regression model_ a "_good_" model?** 1. Is the model wrong? 2. Is the model strong? 3. Does the model produce accurate predictions? 4. Is the model fair? We will review these concepts through today's exercises. A detailed overview is provided in the Unit 1 "Motivating Question" section on the course website. \ # Notes: R code {-} ## Intro to `tidymodels` {.unnumbered .smaller} Throughout the semester, we are going to use the **tidymodels** package in R. - Similar flavor to **tidyverse** structure - More general structure that allows us to fit many other types of models \ At first, it will seem like *a lot more code* (perhaps even *unnecessarily so*). \ For example, what you did in STAT 155 with ```{r} #| eval: false #| echo: true lm(y ~ x1 + x2, data = sample_data) ``` will now look like ```{r tidymodels-example} #| eval: false #| echo: true # STEP 1: model specification lm_spec <- linear_reg() %>% # we want a linear regression model set_mode("regression") %>% # this is a regression task (y is quantitative) set_engine("lm")# we'll estimate the model using the lm function # STEP 2: model estimation model_estimate <- lm_spec %>% fit(y ~ x1 + x2, data = sample_data) ``` \ For now, you'll need to *trust me* that this will be end up being useful. (It will be!) \ \ ## Quick Highlight {.unnumbered .smaller} A few useful functions to use on `model_estimate`: ```{r eval=FALSE,echo=TRUE} model_estimate %>% tidy() #gives you coefficients (and se, t-statistics) ``` ```{r eval=FALSE,echo=TRUE} model_estimate %>% augment(new_data = sample_data) # gives you predictions and residuals for sample_data ``` ```{r eval=FALSE,echo=TRUE} model_estimate %>% glance() #gives you some model evaluation metrics (is it strong?) ``` ```{r eval=FALSE,echo=TRUE} model_estimate %>% augment(new_data = sample_data) %>% # get predictions, and then... mae(truth = y, estimate = .pred) # calculate MAE to measure accuracy of predictions ``` More info, for future reference, below! \ \ ## Future Reference {.unnumbered .smaller} :::{.callout-note} This section is for **future** reference. It is a summary of code you'll learn below for building and evaluating regression models. You do not need to (and in fact should not) run the code in this section verbatim in R; it is example code and meant for future reference only. ::: **Building and evaluating regression models using `tidymodels`.** Throughout, suppose we wish to build and evaluate a linear regression model of `y` vs `x1` and `x2` using a dataset called `sample_data`. You'll need both of these packages: ```{r eval = FALSE} # Load packages library(tidyverse) library(tidymodels) ``` **Building a linear regression model** ```{r eval = FALSE} # STEP 1: specify the type of model to build lm_spec <- linear_reg() %>% # we want a linear regression model set_mode("regression") %>% # this is a regression task (y is quantitative) set_engine("lm") # we'll estimate the model using the lm function # STEP 2: estimate the specified model using sample data model_estimate <- lm_spec %>% fit(y ~ x1 + x2, data = sample_data) # Get the model coefficients model_estimate %>% tidy() ``` **Obtaining predictions (& residuals) for each observation** ```{r eval = FALSE} # Obtain y predictions and residuals for each observation in our sample_data # (We can replace sample_data with any data frame that includes y, x1, and x2) model_estimate %>% augment(new_data = sample_data) # Obtain y predictions (but not residuals) for some given x1, x2 values, when we haven't yet observed y # (We can replace the data.frame with any data frame that includes x1 and x2) model_estimate %>% augment(new_data = data.frame(x1 = ___, x2 = ___)) # Another approach using predict() model_estimate %>% predict(new_data = data.frame(x1 = ___, x2 = ___)) ``` **Evaluating the model** ```{r eval = FALSE} # Is it strong? (R^2) model_estimate %>% glance() # Does it produce accurate predictions? (MAE) model_estimate %>% augment(new_data = sample_data) %>% mae(truth = y, estimate = .pred) # Is it wrong? (residual plot) model_estimate %>% augment(new_data = sample_data) %>% ggplot(aes(x = .pred, y = .resid)) + geom_point() + geom_hline(yintercept = 0) ``` \ # Exercises {-} ## Instructions {.smaller .unnumbered} :::{.incremental} - Work through these exercises as a group, talking through your ideas, questions, and reasoning as you go and taking notes in your QMD. - **Focus on patterns in code.** Review, but do not try to *memorize* any provided code. Focus on the general steps and patterns. - If you're given some starter code with blanks (e.g. below), don't type in those chunks. Instead, copy, paste, and modify the starter code in the chunk below it. ```{r example-startercode} #| eval: false # Start small: rides vs temp ggplot(___, aes(y = ___, x = ___)) + geom___() ``` - **Ask questions!** We will not have time to discuss all exercises at the end of class. Talk through your questions as a group, and ask me questions as I walk around the room! - **Be kind to yourself/each other!** You will be rusty and make mistakes, and that's ok! Mistakes are important to learning. - If you're feeling rusty on STAT 155 material, in particular, check out the *STAT 155 Review* Appendix on our course website. - **Collaborate.** We're sitting in groups for a reason. Collaboration improves higher-level thinking, confidence, communication, community, and more. I expect you to: - Actively contribute to discussion (don't work on your own). - Actively include all group members in discussion. - Create a space where others feel comfortable making mistakes & sharing their ideas (remember that we all come to this class with different experiences, both personal and academic). - Stay in sync while respecting that everybody has different learning strategies, work styles, note taking strategies, etc. If some people are working on exercise 10 and others are on exercise 2, that's probably *not* a good collaboration. - Don't rush. You won't hand anything in and can finish up outside of class. ::: \ ## Questions {.smaller .unnumbered} Capital Bikeshare provides a bike-sharing service in the Washington DC area. Customers can pick up and drop off bikes at any station around the city. Of primary interest to the company is: **How many registered riders can we expect today?** To this end, you will build, evaluate, and compare 2 different **linear regression models** of ridership using the following Capital Bikeshare dataset (originally from the UCI Machine Learning Repository): - Model 1: `rides ~ windspeed + temp` - Model 2: `rides ~ windspeed + temp + weekend` \ \ 0. **Explore the data** Take a minute to get familiar with the data. ```{r} # Load packages we'll need to wrangle and plot the data library(tidyverse) # Load the data bikes <- read.csv("https://Mac-Stat.github.io/data/bike_share.csv") # Only keep / select some variables # And round some variables (just for our demo purposes) bikes <- bikes %>% rename(rides = riders_registered, temp = temp_feel) %>% mutate(windspeed = round(windspeed), temp = round(temp)) %>% select(rides, windspeed, temp, weekend) ``` ```{r} # Check out the dimensions dim(bikes) # Check out the first 3 rows head(bikes, 3) ``` This dataset contains the following information for a sample of different dates: | variable | description | |:-------------|:---------------------------------------------| | rides | count of daily rides by registered users | | windspeed | wind speed in miles per hour | | temp | what the temperature feels like in degrees Fahrenheit | | weekend | whether or not it falls on a weekend | \ \ 1. **Plot the relationships** First, let's *plot* these relationships. *REMEMBER: Don't write in any chunk with starter code. Copy, paste, and modify the code in the chunk below it.* ```{r q1a-starter} #| eval: false # Start small: rides vs temp ggplot(___, aes(y = ___, x = ___)) + geom___() ``` ```{r} ``` ```{r q1b-starter} #| eval: false # rides vs temp & windspeed ggplot(bikes, aes(y = ___, x = ___, ___ = windspeed)) + geom_point() ``` ```{r} ``` ```{r q1c-starter} #| eval: false # rides vs temp & windspeed & weekend ggplot(bikes, aes(y = ___, x = ___, ___ = windspeed)) + geom_point() + facet_wrap(~ ___) ``` ```{r} ``` \ \ 2. **tidymodels STEP 1** We'll build and evaluate our two models of ridership using the **tidymodels** package. This code is more complicated than the `lm()` function we used in STAT 155. BUT: - **tidymodels** is part of the broader **tidyverse** (what we use to plot and wrangle data), thus the syntax is more consistent - **tidymodels** generalizes to the other ML algorithms we'll survey in this course, thus will eventually *minimize* the unique syntax we need to learn ```{r} # Load package library(tidymodels) ``` The first step is to **specify what type of model we want to build**. We'll store this as `lm_spec`, our linear regression model (`lm`) specification (`spec`). ```{r} lm_spec <- linear_reg() %>% # we want a linear regression model set_mode("regression") %>% # this is a regression task (y is quantitative) set_engine("lm") # we'll estimate the model using the lm function ``` This code *specifies* but doesn't *build* any model -- we didn't even give it any data or specify the variables of interest! ```{r} # Check it out lm_spec ``` \ \ 3. **tidymodels STEP 2** We can now estimate or `fit` our two ridership models using the specified model structure (`lm_spec`) and our sample `bikes` data: ```{r} # Fit bike_model_1 bike_model_1 <- lm_spec %>% fit(rides ~ windspeed + temp, data = bikes) # Check out the coefficients bike_model_1 %>% tidy() ``` ```{r} # YOUR TURN # Fit bike_model_2 & check out the coefficients ``` \ \ 4. **Is it fair?** Now, let's *evaluate* our two models. First, do you have any concerns about the context in which the data were collected and analyzed? About the potential impact of this analysis? \ \ 5. **Is it strong?** We can measure and compare the strength of these models using $R^2$, *the proportion of variability in our response variable that's explained by the model*. Report which model is stronger and interpret its $R^2$. ```{r} # Obtain R^2 for bike_model_1 bike_model_1 %>% glance() ``` ```{r} # YOUR TURN # Obtain R^2 for bike_model_2 ``` \ \ 6. **Pause -- predictions and residuals** Our next model evaluation questions will focus on the models' predictions and prediction *errors*, or *residuals*. We can obtain this information by *augmenting* our models with our original `bikes` data. For example: ```{r} # Calculate predicted ridership (.pred) & corresponding residuals (.resid) using bike_model_1 # Just look at first 6 days bike_model_1 %>% augment(new_data = bikes) %>% head() ``` We can also predict outcomes for *new* observations using either `augment()` or `predict()`. Note the difference in the output: ```{r} # Predict ridership on a 60 degree day with 20 mph winds bike_model_1 %>% augment(new_data = data.frame(windspeed = 20, temp = 60)) ``` ```{r} # Predict ridership on a 60 degree day with 20 mph winds bike_model_1 %>% predict(new_data = data.frame(windspeed = 20, temp = 60)) ``` \ \ 7. **Does it produce accurate predictions?** Recall that the **mean absolute error** (MAE) measures the typical prediction error. Specifically, it is the *mean* of the *absolute* values of the residual *errors* for the days in our dataset. a. Use the residuals to calculate the MAE for the 2 models. HINT: `abs()` ```{r eval = FALSE} # DON'T TYPE IN THIS CHUNK # MAE for bike_model_1 bike_model_1 %>% augment(new_data = bikes) %>% summarize(mae = ___(___(___))) ``` ```{r} # MAE for bike_model_1 ``` ```{r} # YOUR TURN: MAE for bike_model_2 ``` b. Doing the calculation from scratch helps solidify your understanding of how MAE is calculated, thus interpreted. Check your calculations using a shortcut function. ```{r} # Calculate MAE for the first model bike_model_1 %>% augment(new_data = bikes) %>% mae(truth = rides, estimate = .pred) ``` ```{r} # YOUR TURN # Calculate MAE for the second model ``` c. Which model has more accurate predictions? d. Interpret the MAE for this model, in context, and comment on whether it's "large" or "small". NOTE: "large" or "small" is defined by the context (e.g. relative to the observed range of ridership, the consequences of a bad prediction, etc). \ \ 8. **Is it wrong?** To determine whether the linear regression assumptions behind `bike_model_1` and `bike_model_2` are reasonable, we can review residual plots, i.e. plots of the residuals vs predictions for each observation in our dataset. Run the code below and summarize your assessment of whether our models are wrong. RECALL: We want the points to appear *random* and *centered around 0* across the entire range of the model / predictions. ```{r} # Residual plot for bike_model_1 bike_model_1 %>% augment(new_data = bikes) %>% ggplot(aes(x = .pred, y = .resid)) + geom_point() + geom_hline(yintercept = 0) ``` ```{r} # YOUR TURN # Residual plot for bike_model_2 ``` \ \ 9. **Art vs science** Inspecting residual plots is more art than science.^[Stefanski, Leonard A. (2007). Residual (Sur)Realism. "The American Statistician," 61, pp 163-177.] It requires a lot of practice. Consider another example using simulated data. First, build a model that assumes all predictors are roughly linearly related: ```{r} # Import data simulated_data <- read.csv("https://Mac-Stat.github.io/data/simulated_data.csv") # Model y by the 6 input variables new_model <- lm_spec %>% fit(y ~ x1 + x2 + x3 + x4 + x5 + x6, simulated_data) ``` Next, check out a pairs plot. Is there anything here that makes you think that our model assumption is bad? (*NOTE: if you are not familiar with the `ggpairs` function we are using below, chat with your group and look for resources to figure out what this is doing!*) ```{r fig.height = 5, fig.width = 5} library(GGally) ggpairs(simulated_data) ``` Finally, check out a residual plot. Any concerns now? ```{r eval = FALSE} new_model %>% ___(new_data = ___) %>% ggplot(aes(x = ___, y = ___)) + geom_point(size = 0.1) + geom_hline(yintercept = 0) ``` ```{r} ``` \ \ 10. **Details -- communication & code style** Communication is a key machine learning skill, including written summaries, presentations, and **code**. Just like an essay, code must have structure, signposts, and grammar that will make it easier for others to follow. The below code runs, but it is "bad code". - Fix this code and add **comments** so that it is easier for yourself and others to follow. - Also pay attention to what this code *does*. ```{r} bikes%>%group_by(weekend)%>%summarize(median(rides)) ``` ```{r} mynewdatasetissmallerthantheoriginal<-bikes%>%filter(rides<=700,weekend==FALSE,temp>60) mynewdatasetissmallerthantheoriginal ``` ```{r} mynewdatasetusescelsius<-bikes%>%mutate(temp=(temp-32)*5/9) head(mynewdatasetusescelsius) ``` \ \ 11. **STAT 155 Review -- model interpretation & application** Let's **interpret** and **apply** `bike_model_2`. ```{r eval = FALSE} ___ %>% tidy() ``` ```{r} ``` a. How can we interpret the `temp` coefficient? - We expect roughly 54 more riders on warm days. - We expect roughly 54 more riders per every 1 degree increase in temperature. - When controlling for windspeed and weekend status, we expect roughly 54 more riders on warm days. - When controlling for windspeed and weekend status, we expect roughly 54 more riders per every 1 degree increase in temperature. b. How can we interpret the `weekendTRUE` coefficient? - We expect roughly 858 fewer riders on weekends. - We expect roughly 858 fewer riders per every extra weekend. - When controlling for windspeed and temperature, we expect roughly 858 fewer riders on weekends. - When controlling for windspeed and temperature, we expect roughly 858 fewer riders per every extra weekend. c. Reproduce the predicted ridership and corresponding residual for day 1 from scratch (how were these calculated?!): ```{r eval = FALSE} bike_model_2 %>% ___(new_data = bikes) %>% head(1) ``` ```{r} ``` \ \ 12. **STAT 155 Review -- data wrangling** Through the "Details: communication & code style" and elsewhere, you've reviewed the use of various `dplyr` data wrangling verbs: `filter()`, `mutate()`, `summarize()`, `group_by()`, `select()`, `arrange()`. Use these to complete the following tasks. a. Calculate the mean temperature across all days in the data set. ```{r} ``` b. Calculate the mean temperature on weekends vs weekdays. ```{r} ``` c. Print out the 3 days with the *highest* temperatures. HINT: `arrange()` or `arrange(desc())` ```{r} ``` d. Name and store a new data set which: - only includes the days that fall on a weekend *and* have temps below 80 degrees - has a new variable, `temp_above_freezing`, which calculates how far the temperature is above (or below) freezing (32 degrees F) - only includes the `windspeed`, `temp`, and `temp_above_freezing` variables. ```{r} ``` \ \ 13. **STAT 155 Review -- plots** Construct plots of the following relationships: ```{r} # rides vs temp ``` ```{r} # rides vs weekend ``` ```{r} # rides vs temp and weekend ``` \ \ # Done! - Render your notes. - Check the solutions in the course website (Solution drop downs). - If you finish all that during class, start your homework!