--- title: "Simple linear regression: Formalization (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 = "Organize your files"} This **qmd file** is where you'll type notes, code, etc. Directions: - Save this file in the `inclass_activities` sub-folder of the `stat155` folder you created before today's class. Use a file name related to the activity number and/or today's date (eg: "activity_4" or "4_simple_linear_regression"). ::: ## Warm-up {-} **Guiding question** How is bikeshare ridership related to various weather factors? \ **Context** We'll explore data from [Capital Bikeshare](https://capitalbikeshare.com/), a company in Washington DC. Our main goal will be to explore daily ridership among *registered* users, as measured by `rides`. ```{r} # Load packages and import data # Shorten the "riders_registered" variable to "rides" # Keep only certain columns / variables library(tidyverse) bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% rename(rides = riders_registered) %>% select(date, day_of_week, weekend, temp_feel, humidity, windspeed, rides) ``` \ \ ### EXAMPLE 1: Get to know the data ```{r} # Check out the first 6 rows # What does a case represent, i.e. what's the unit of observation? # How much data do we have? ``` **NOTE** It's important to think about the "who, what, when, where, why, and how" of a dataset. In this example, the *where* and *when* are particularly relevant. These data were collected in D.C. in 2011-12. The patterns we observe may not extend to other places and/or time periods! \ \ \ \ ### EXAMPLE 2: Explore ridership Let's get acquainted with the `rides` variable. Remember: Code = communication. Pay special attention to formatting / readability. ```{r} # Calculate the average (central tendency) and standard deviation (variability) in daily ridership bikes %>% ___(mean(___), sd(___)) ``` ```{r} # Construct a plot of how ridership varies from day to day bikes %>% ___(aes(___)) ``` **Follow-up** Summarize, in words, what you learned from the plot and numerical summaries. Remember to comment on: *shape*, *central tendency*, *spread* / *variability*, and, if relevant, any *outliers*. \ \ \ \ ### EXAMPLE 3: Relationships Now that we better understand ridership (`rides`), natural and important follow-up questions in nearly every statistical analysis are: What *explains* why ridership varies from day to day, i.e. why there are more riders on some days than others?! For example, is ridership *related to* / can it be *predicted by* various weather factors such as `windspeed` in miles per hour (mph), or the "feels like" temperature in degrees Fahrenheit (`temp_feel`)? In this setting: - response variable Y = `rides` - potential predictors X = `windspeed`, `temp_feel` Let's start with the relationship of `rides` with `windspeed`. #### Part a: intuition Quick: Do you anticipate that `rides` is related to `windspeed`? If so, do you anticipate that the relationship is positive or negative? Strong, weak, or moderate? #### Part b: visualization It's always good to start our analysis with a picture! We can *visualize* the relationship between these 2 quantitative variables using a **scatterplot**. This will require some new code, but let's start with what we know. ```{r} bikes %>% ggplot(aes(___)) ``` Add the following two lines to your plot. ```{r} # Add a blue simple linear regression line (line of best fit) # and a orange *curve* of best fit PUT YOUR PLOT HERE + geom_smooth(method = "lm", se = FALSE) + geom_smooth(color = "orange", se = FALSE) ``` #### Part c: description Summarize the relationship. Remember to comment on: *shape*, *direction*, *strength*, and, if relevant, any *outliers*. \ \ \ \ ### EXAMPLE 4: Correlation Some of our above description is a bit vague, and somewhat subjective. We can more precisely *numerically measure* the *strength* and *direction* of a *linear* relationship using **correlation**. Guess the correlation between ridership and windspeed. (You'll calculate this in the exercises below.) \ \ \ \ ### EXAMPLE 5: Simple linear regression (Part 1) We can also more accurately summarize this relationship by finding a formula for the linear trend: ```{r} bikes %>% ggplot(aes(y = rides, x = windspeed)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` Consider the *sample* linear regression model of this relationship, i.e. the formula for the line in the scatterplot: E[ridership | windspeed] = $\hat{\beta}_0$ + $\hat{\beta}_1$ windspeed \ \ In the exercises, you'll show that: E[ridership | windspeed] = 4490.10 - 65.34 windspeed Consider the **intercept coefficient** 4490.10. Physically, this is the *y-intercept* of our model line (where it crosses the y-axis). What's the *contextual* meaning? a. A 1 mph increase in wind speed is associated with a 4490.10 increase in expected ridership. b. On days with wind speeds of 0 mph (no wind!), we expect 4490.10 riders. c. On average, there are 4490.10 riders per day. \ \ \ \ ### EXAMPLE 6: Simple linear regression (Part 2) Consider the **windspeed coefficient**, -65.34 Physically, this is the *slope* of our model line (rise over run). Pick the most effective / correct *interpretations* of this slope below. For each other option, indicate what's wrong with the interpretation. #### Part 1 a. Increasing X by 1 is associated with a 65.34 decrease in the average Y value. b. A 1 mph increase in windspeed is associated with a decrease of 65.34 riders. c. A 1 mph increase in windspeed is associated with a decrease of 65.34 in average ridership. d. Windier weather discourages ridership. A 1 mph increase in windspeed produces 65.34 fewer riders on average. e. If tomorrow is 1 mph warmer than today, we should expect 65.34 fewer riders. f. A windspeed change of 1 is associated with a 65.34 decrease in average ridership. g. A 1 rider increase is associated with a 65.34 mph decrease in average windspeed. #### Part 2 Keeping in mind the *units* and *range* of our variables, what do you think about the *magnitude* of the slope? Is a decrease of 65.34 in average ridership per 1-mph increase in windspeed *meaningful* in this context? \ \ \ \ ::: {.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 ::: ::: {.callout-note title = "Interpreting intercept and slope"} Language is fluid -- there are many ways to phrase a correct interpretation of the intercept and slope! Here's one general option, that can be translated into the context at hand: E[Y|X] = a + bX - Intercept: When X is 0 units, we expect Y to be a. - Slope: A 1-unit increase in X is associated with a b-unit (increase / decrease) in the average Y value. ::: \ \ \ \ ## Exercises {-} **DIRECTIONS** You'll work on these exercises in your groups. Collaboration is a key learning goal in this course. - Why? Discussion & collaboration deepens your understanding of a topic while also improving confidence, communication, community, & more. (eg: Deeply learning a new language requires more than working through Duolingo alone. You need to talk with and listen to others!) - How? You are expected to: - Use your group members' names & pronouns. It's ok to ask if you don’t remember! - Actively contribute to discussion. Don't work on your own. - Actively include all other group members in discussion. - Create a space where others feel comfortable sharing ideas & questions. - We won't discuss these exercises as a class. With that, when you get stuck: - Carefully re-read the problem. Make sure you didn't miss any directions -- it can be tempting to skip words and go straight to an R chunk, but don't :). - Discuss any questions you have with your group. - If the question is unresolved by the group, be sure to ask the instructor! - Remember that there are solutions in the online manual, at the bottom of the activity. \ \ \ \ ### Exercise 1: Correlation In the warm-up, you took a guess at the correlation between ridership and windspeed: ```{r} bikes %>% ggplot(aes(y = rides, x = windspeed)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` Let's check your intuition! Fill in the code below to calculate the correlation. ```{r} # HINT: What function is useful for calculating statistics (eg: means)? bikes %>% ___(cor = cor(rides, windspeed)) ``` \ \ \ \ ### Exercise 2: Build the model You were told above that our (sample) simple linear regression model was: E[ridership | windspeed] = 4490.10 - 65.34 windspeed Verify this using the code below. Step through code chunk slowly, make note of new code. ```{r} # Construct and save the model as bike_model_1 # What's the purpose of "rides ~ windspeed"? # What's the purpose of "data = bikes"? bike_model_1 <- lm(rides ~ windspeed, data = bikes) ``` ```{r} # A long summary of the model stored in bike_model_1 summary(bike_model_1) ``` ```{r} # A simplified model summary coef(summary(bike_model_1)) ``` \ \ \ \ ### Exercise 3: Predictions and residuals On March 8, 2012, the `windspeed` was 29.58472 mph and there were 4896 riders: ```{r} # Note the use of filter() and select() # You'll learn more about filter below! bikes %>% filter(date == "2012-03-08") %>% select(rides, windspeed) ``` #### Part a In the scatterplot below, locate the point that corresponds to March 8, 2012 (by eye, not using code). Is it close to the trend? Were there *more* riders than expected or *fewer* than expected? ```{r} bikes %>% ggplot(aes(y = rides, x = windspeed)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` #### Part b Using the model formula below, calculate how many riders we would have *expected* or *predicted* on March 8, 2012 based on its windspeed: E[ridership | windspeed] = 4490.10 - 65.34 windspeed ```{r} # Calculate the prediction IN THIS CHUNK # (Convince yourself that you can connect this to the scatterplot) 4490.10 - 65.34*___ ``` #### Part c Check your prediction using the `predict()` function. ```{r} # What is the purpose of including bike_model_1? # What is the purpose of newdata = ___??? predict(bike_model_1, newdata = data.frame(windspeed = 29.58472)) ``` #### Part d How far does the *observed* ridership on March 8, 2012 fall from its *model prediction*? Calculate the **residual** or **prediction error**: ```{r} # Recall: residual = observed Y - predicted Y ``` #### Part e Are *positive* residuals above or below the trend line? When we have *positive* residuals, does the model over- or under-estimate ridership? Are *negative* residuals above or below the trend line? When we have *negative* residuals, does the model over- or under-estimate ridership? \ \ \ \ ### Exercise 4: Ridership vs temp_feel (visualization) Your turn! Let's practice what you've learned in the preceding exercises to investigate the relationship between `rides` and `temp_feel`. #### Part a Below, you'll plot the relationship of `rides` and `temp_feel`. What do you *expect* this to look like? Will it be... - positive / negative? - linear? - strong / moderate / weak? #### Part b ```{r} # Construct a visualization of the relationship of rides with temp_feel # Include a representation of: # - the individual data points # - a simple linear regression line (line of best fit) # - a *curve* of best fit ``` #### Part c Summarize, in words, your findings from the plot. Remember: shape, direction, strength, and, if relevant, outliers. \ \ \ \ ### Exercise 5: Filtering our data The relationship of `rides` with `temp_feel` is nonlinear! Though we'll soon learn techniques for handling nonlinear relationships, right now let's focus on days with feels-like temperatures below 80 degrees (where the relationship is linear). In the `tidyverse` package, whereas `select()` subsets our data to include only certain *columns / variables* of interest, `filter()` subsets our data to include only certain *rows / cases* of interest. Let's practice: ```{r} bikes %>% filter(rides > 6911) ``` ```{r} bikes %>% filter(rides >= 6911) ``` ```{r} bikes %>% filter(rides == 6911) ``` Your turn: show only days where `temp_feel` is under 80 degrees. ```{r} # There are many such days! Show only the first 6. bikes %>% filter(temp_feel ___) %>% head() ``` Now that you've convinced yourself that the filter works, store the cases where `temp_feel` is under 80 degrees as a new dataset called `bikes_sub`: ```{r} # Don't use head() here! Otherwise we'd only store the first 6 dates. bikes_sub <- bikes %>% filter(temp_feel ___) ``` \ \ \ \ ### Exercise 6: Ridership vs temp_feel (model) Now that we've filtered the data, let's explore the linear relationship of `rides` with `temp_feel` on dates with temperatures below 80 degrees. NOTE: We'll be using the `bikes_sub` data now!!! ```{r} bikes_sub %>% ggplot(aes(x = temp_feel, y = rides)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` #### Part a ```{r} # Use lm() to construct a model of rides vs temp_feel # Save this as bike_model_2 # !!!Remember to use the bikes_sub data!!! ``` ```{r} # Get a short summary of bike_model_2 ``` #### Part b Write out a formula for the sample model trend. E[??? | ???] = ??????? #### Part c Interpret the intercept *and* the `temp_feel` coefficients. NOTE: Your interpretation of the intercept should be "correct" but won't make sense here since it **extrapolates** our model far beyond the range of our observed data. #### Part d How well would our model predict the ridership on October 29, 2012? Follow the steps below to calculate the residual on this date. Interpret this residual in words! ```{r} # Obtain the temp_feel and ridership on October 29, 2012 # (Once you obtain that info, identify the corresponding point on the scatterplot!) # Predict ridership on this date using bike_model_2 (hence the temperature info) # NOTE: Use predict() but also convince yourself you could do this from scratch # Calculate the residual ``` \ \ \ \ ### Exercise 7: Data drill (`filter`, `select`, `summarize`) This exercise is designed to help you practice your tidy data wrangling skills. We can't do statistics or really understand data without these skills! We'll work with a simpler set of 10 data points: ```{r} new_bikes <- bikes %>% select(date, temp_feel, humidity, rides, day_of_week) %>% head(10) ``` #### Part a Thus far, in the `tidyverse` grammar you've learned 3 **verbs** or action words: `summarize()`, `select()`, `filter()`. First, refresh your memory. Take note of the following code and the output it produces: ```{r} # First check out the data for a point of reference new_bikes ``` ```{r} # summarize() new_bikes %>% summarize(mean(temp_feel), mean(humidity)) ``` ```{r} # select() new_bikes %>% select(date, temp_feel) ``` ```{r} # select() again new_bikes %>% select(-date, -temp_feel) ``` ```{r} # filter() new_bikes %>% filter(rides > 850) ``` ```{r} # filter() again new_bikes %>% filter(day_of_week == "Sat") ``` ```{r} # filter() again new_bikes %>% filter(rides > 850, day_of_week == "Sat") ``` #### Part b Your turn. Use `tidyverse` verbs to complete each task below. ```{r} # Keep only information about the humidity and day of week # Keep only information about the humidity and day of week using a different approach # Calculate the maximum and minimum temperatures # Keep only information for Sundays # Keep only information for Sundays with temperatures below 50 # Calculate the median ridership on Sundays with temperatures below 50 ``` \ \ \ \ ### Reflection **Learning goals** Recall the learning goals for today's activity. Explore numerical and visual summaries for the relationship between a *quantitative* **response (aka outcome) variable** and a *quantitative* **predictor (aka explanatory variable)**: - Construct and interpret a **scatterplot visualization** of the relationship. - Compute and interpret the linear **correlation** between the 2 variables. - Analyze and apply a simple linear regression model of the relationship: - obtain and write a **model formula** - interpret the model's **intercept** and **slope** coefficients *in context* - compute and interpret **predictions** (aka **expected** or **fitted values**) and **residuals** *in context* - explain the connection between residuals and the **least squares criterion** Among these learning goals: - Which do you have a good handle on? - Which are struggling with? What feels challenging right now? - What are some wins from the day? - Statistics is a particular kind of language and collection of tools for channeling **curiosity** to improve our world. How do the concepts we practiced today facilitate curiosity? **Code** In addition to exploring the learning goals, you learned some new code. - If you haven't already, you're highly encouraged to start tracking and organizing new code in a cheat sheet (eg: a Google doc). This will be a handy reference for you, and the act of making it will help deepen your understanding and retention. - Reflect upon the `lm()` function. This function takes 2 arguments. What are they and what's their purpose? Mainly, what information does R need in order to build a model? - Reflect upon the `predict()` function. This function takes 2 arguments. What are they and what's their purpose? Mainly, what information does R need in order to make a prediction? \ \ \ \ ## Extra practice After class, you're encouraged to work through the optional exercises below. These explore the *bills* (beaks) of Chinstrap penguins: ```{r} data(penguins) chinstraps <- penguins %>% filter(species == "Chinstrap") ``` Specifically, let's explore the relationship of bill length with bill depth, both measured in mm: ```{r} # Visual summary chinstraps %>% ggplot(aes(y = bill_len, x = bill_dep)) + geom_point() + geom_smooth(method = "lm", se = FALSE) ``` ```{r} # Model summary bill_model <- lm(bill_len ~ bill_dep, data = chinstraps) coef(summary(bill_model)) ``` ### Exercise 8: model formulas Which of the following is the correct (sample) model formula of this relationship: a. bill_dep = 1.92 + 13.43 bill_len b. bill_len = 1.92 + 13.43 bill_dep c. E[bill_dep | bill_len] = 1.92 + 13.43 bill_len d. E[bill_len | bill_dep] = 1.92 + 13.43 bill_dep e. bill_dep = 13.43 + 1.92 bill_len f. bill_len = 13.43 + 1.92 bill_dep g. E[bill_dep | bill_len] = 13.43 + 1.92 bill_len h. E[bill_len | bill_dep] = 13.43 + 1.92 bill_dep \ \ \ \ ### Exercise 9: Coefficient interpretations #### Part a Interpret the slope: A 1 ___ increase in ___ is associated with a 1.92 ___ [increase/decrease] in ___ ___. #### Part b Does it make sense to interpret the intercept in context in this example? \ \ \ \ ### Example 10: More predictions **Challenge yourself to use just the model coefficients to answer the below questions.** #### Part a Predict the bill length of a penguin that has a bill depth of 0cm (which doesn't make sense :)). #### Part b You come across 2 groups of penguins. In group "A", all penguins have a bill depth of 19cm. In group "B", all penguins have a bill depth of 20cm. Do you expect that the average bill length in group B is greater or less than that of group A? How *much* greater / less (provide a specific number). #### Part c Among penguins with bill depths of 17mm, the average bill length is 46.07mm. What's the average bill length among penguins with bill depths of 18mm? \ \ \ \ ### Example 11: Correlation Using only the scatterplot, what do you think is the correlation between bill_len and bill_dep: -1, -0.95, -0.65, -0.15, 0, 0.15, 0.65, 0.95, or 1? ## 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 `inclass_activities` subfolder within your `stat155` folder and locate the HTML file -- you can open it again in your browser. - Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' pane. - Check the solutions in the course website, at the bottom of the corresponding chapter. - Work on homework and/or any extra practice exercises!