--- title: "Simple Linear Regression - Categorical Predictor" subtitle: "Notes and in-class exercises" format: html: embed-resources: true toc: 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') ``` You can download the .qmd file for this activity [here](../activity_templates/07-slr-cat-predictor.qmd) and open in R-studio. The rendered version is posted in the [course website](https://mutasim221b.github.io/Mac-STAT-155-Sp-26/) (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins. # Notes ## Learning goals By the end of this lesson, you should be able to: - Write a **model formula** for a simple linear regression model with a categorical predictor using **indicator variables** - Interpret the coefficients in a simple linear regression model with a categorical predictor ## Readings and videos Complete *both* the reading and the videos to go through **before** class. - Reading: Section 3.9 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/) only up through section 3.9.1 Indicator Variables - Videos: - [Simple linear regression: categorical predictor](https://youtu.be/QPpp6Kpk2To) ([slides](https://drive.google.com/file/d/1fbSJW4DhVWxR8rykowPUoAZ6WfzuZ857/view?usp=sharing)) - [R Code for Categorical Predictors](https://voicethread.com/share/14980390/) **File organization:** Save this file in the "Activities" subfolder of your "STAT155" folder. # Exercises **Context:** Today we'll explore data on thousands of diamonds to understand how physical characteristics relate to price. Read in the data below. ```{r warning=FALSE, message=FALSE} # Load packages and import data library(tidyverse) data(diamonds) # A little bit of data wrangling code - let's not focus on this for now diamonds <- diamonds %>% mutate( cut = factor(cut, ordered = FALSE), color = factor(color, ordered = FALSE), clarity = factor(clarity, ordered = FALSE) ) ``` For the first several exercises, our focus will be on the relationship between diamond price and cut. ## Exercise 1: Get to know the data Write R code to answer the following: - a. How many cases and variables do we have? What does a case represent? - b. What do the first few rows of the data look like? - c. Construct and interpret two different visualizations of the `price` variable. - d. Construct and interpret a visualization of the `cut` variable. ```{r eval = TRUE} dim(diamonds) head(diamonds) # Visualize price (outcome variable) ggplot(diamonds, aes(x = price)) + geom_histogram() ggplot(diamonds, aes(y = price)) + geom_boxplot() diamonds %>% summarize(mean(price), median(price), sd(price)) # Visualize cut (predictor variable) ggplot(diamonds, aes(x = cut)) + geom_bar() diamonds %>% count(cut) ``` ## Exercise 2: Visualizations Start by visualizing this relationship of interest, that between `price` and `cut`. a. The appropriate plot depends upon the type of variables we're plotting. When exploring the relationship between a quantitative response and a quantitative predictor, a scatterplot was an effective choice. After running the code below, explain why a scatterplot is not effective for exploring the relationship between the outcome price and *categorical* `cut` predictor. ```{r} # Try a scatterplot ggplot(diamonds, aes(y = price, x = cut)) + geom_point() ``` > **Response:** Put your response here. b. Separately run each chunk below, with two plots. Comment (`#`) on what changes in the code / output. b.1. ```{r} # Univariate boxplot ggplot(diamonds, aes(y = price)) + geom_boxplot() ``` > **Response:** Put your response here. b.2. ```{r} # ??? ggplot(diamonds, aes(y = price, x = cut)) + geom_boxplot() ``` > **Response:** Put your response here. b.3. ```{r} # Univariate density plot ggplot(diamonds, aes(x = price)) + geom_density() ``` > **Response:** Put your response here. b.4. ```{r} # Comparisons in density plot ggplot(diamonds, aes(x = price, color = cut)) + geom_density() ``` > **Response:** Put your response here. b.5. ```{r} # Univariate histogram ggplot(diamonds, aes(x = price)) + geom_histogram() ``` > **Response:** Put your response here. b.6. What’s the difference between this and b4 density plot comparison? Can we now interpret the relationship between price and cut as we could in density plot? Why not? ```{r} # ??? ggplot(diamonds, aes(x = price)) + geom_histogram() + facet_wrap(~ cut) ``` > **Response:** Put your response here. c. Do you notice anything interesting about the relationship between price and cut? What do you think might be happening here? > **Response:** Put your response here. ## Exercise 3: Numerical summaries Let's follow up our plots with some numerical summaries. a. To warm up, first calculate the mean `price` across all diamonds. ```{r} diamonds %>% summarize(mean(price)) ``` b. To summarize the trends we observed in the grouped plots above, we can calculate the mean `price` for each type of `cut`. This requires the inclusion of the `group_by()` function: ```{r} # Calculate mean price by cut diamonds %>% group_by(cut) %>% summarize(mean(price)) ``` c. Examine the group mean measurements. can you match these numbers up with what you see in the plots? d. Based on the results above, we can see that, on average, diamonds with a "Fair" cut tend to cost more than higher-quality cuts. Let's construct a new variable named `cutFair`, using on the following criteria: - `cutFair` = 1 if the diamond is of Fair cut - `cutFair` = 0 otherwise (any other value of cut (Good, Very Good, Premium, Ideal)) The `ifelse` function allows to create a new variable from an existing one, based on whether or not the values in that variable meet a certain "condition" (remember, you can always look up function documentation in R by typing `?ifelse` in the Console, and hitting enter!). Fill in the following code to create `cutFair`. The condition was given to you already. Try to use this to complete the code. ```{r} # In the first blank, put what value cutFair should have if the condition is "met", or TRUE # In the second blank, put what value cutFair should have if the condition is "not met", or FALSE diamonds <- diamonds %>% mutate(cutFair=ifelse(cut == "Fair", 1, 0)) ``` Variables like `cutFair` that are coded as 0/1 to numerically indicate if a categorical variable is at a particular state are known as an **indicator variable**. You will sometimes see these referred to as a "binary variable" or "dichotomous variable"; you may also encounter the term "dummy variable" in older statistical literature. e. Now, let's calculate the group means based on the new `cutFair` indicator variable: ```{r} diamonds %>% group_by(cutFair) %>% summarize(mean(price)) ``` ## Exercise 4: Modeling trend using a categorical predictor with exactly 2 categories Next, let's *model* the trend in the relationship between the `cutFair` and `price` variables using a simple linear regression model: ```{r} # Construct the model diamond_mod0 <- lm(price ~ cutFair, data = diamonds) # Summarize the model summary(diamond_mod0) ``` Compare these results to the output of exercise 3e. What do you notice? How do you interpret the intercept and `cutFair` coefficient terms from this model? ## Exercise 5: Modeling trend using a categorical predictor with >2 categories Using a single binary predictor like the `cutFair` indicator variable is useful when there are two clearly delineated categories. However, the `cut` variable actually contains 5 categories! Because we've collapsed all non-Fair classifications into a single category (i.e. `cutFair = 0`), the model above can't tell us anything about the difference in expected price between, say, *Premium* and *Ideal* cuts. The good news is that it is very straightforward to model categorical predictors with >2 categories. We can do this by using the `cut` variable as our predictor: ```{r} # Construct the model diamond_mod <- lm(price ~ cut, data = diamonds) # Summarize the model summary(diamond_mod) ``` a. Even though we specified a single predictor variable in the model, we are seeing 4 coefficient estimates--why do you think this is the case? We are seeing 4 coefficient estimates because each category is being assigned to a separate indicator variable--`cutGood` = 1 when `cut == "Good"` and 0 otherwise, `cutVery Good` = 1 when `cut == "Very Good" and 0 otherwise, and so on. NOTE: We see 4 indicator variables (for Good, Very Good, Premium, and Ideal), but we do not see `cutFair` in the model output. This is because `Fair` is the **reference level** of the `cut` variable (it's first alphabetically). b. After examining the summary table output from the code chunk above, complete the model formula: \ E[price | cut] = 4358.7578 - 429.8933 cutGood - 376.9979 cutVery Good + 225.4999 cutPremium - 901.2158 cutIdeal \ ## Exercise 6: Making sense of the model Recall our model: E[price | cut] = 4358.7578 - 429.8933 cutGood - 376.9979 cutVery Good + 225.4999 cutPremium - 901.2158 cutIdeal a. Use the model formula to calculate the expected/typical price for diamonds of *Good* cut. - Expected/typical price for diamonds of *Good* cut: E[price | cut] = 4358.7578 - 429.8933 * 1 - 376.9979 * 0 + 225.4999 * 0 - 901.2158 * 0 = 4358.7578 - 429.8933 = \$3928.865 ```{r eval = TRUE} predict(diamond_mod, newdata = data.frame(cut = "Good")) ``` b. Similarly, calculate the expected/typical price for diamonds of *Fair* cut. - Expected/typical price for diamonds of *Fair* cut: E[price | cut] = 4358.7578 - 429.8933 * 0 - 376.9979 * 0 + 225.4999 * 0 - 901.2158 * 0 = \$4358.7578 ```{r eval = TRUE} predict(diamond_mod, newdata = data.frame(cut = "Fair")) ``` c. Re-examine these 2 calculations. Where have you seen these numbers before?! These come from our group mean calculations in Exercise 3b! The predicted value for diamonds of *Fair* cut is also the same as what we obtained using the SLR model in exercise 4 with only a single `cutFair` indicator variable. ## Exercise 7: Interpreting coefficients Recall that our model formula is *not a formula for a line*. Thus we can't interpret the coefficients as "slopes" as we have before. Taking this into account and reflecting upon your calculations above... a. Interpret the intercept coefficient (`4358.7578`) in terms of the data context. *Make sure to use non-causal language, include units, and talk about averages rather than individual cases.* > The average price of a Fair cut diamonds is \$4358.7578. b. Interpret the `cutGood` and `cutVery Good` coefficients (`-429.8933` and `-376.9979`) in terms of the data context. Hint: where did you use these value in the prediction calculations above? - Interpretation of `cutGood` coefficient: On average, Good cut diamonds are worth \$429.89 less than Fair cut diamonds. - Interpretation of `cutVery Good` coefficient: On average, Very Good cut diamonds are worth \$377.00 less than Fair cut diamonds. ## Exercise 8: Modeling choices (CHALLENGE) Why do we fit this model in this way (using 4 indicator variables `cutGood`, `cutVery Good`, `cutPremium`, `cutIdeal`)? Instead, suppose that we created a single variable `cutCat` that gave each category a numerical value: 0 for Fair, 1 for Good, 2 for Very Good, 3 for Premium, and 4 for Ideal. How would this change things? What are the pros and cons of each approach? ## 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. # Additional Practice ## Exercise 9: The least squares criterion The coefficient estimates `diamond_mod` were selected in the same way as when our predictor is quantitative: by minimizing the sum of the squared residuals. Use this model to calculate the residual for the first diamond in the dataset: ```{r} # Observed data on the first diamond diamonds %>% select(price, cut) %>% head(1) ``` ## Exercise 10: Diamond color Consider modeling `price` by `color`. - Before creating a visualization that shows the relationship between price and color, write down what you expect the plot to look like. Then construct and interpret an appropriate plot. - Compute the average price for each color. - Fit an appropriate linear model with `lm()` and display a short summary of the model. - Write out the model formula from the above summary. - Which color is the reference level? How can you tell from the model summary? - Interpret the intercept and two other coefficients from the model in terms of the data context. ## Exercise 11: Diamond clarity Repeat the steps from the previous exercise for the `clarity` variable. ## Exercise 12: Evaluating model strength Let's study some penguin data! ```{r} data(penguins) penguins <- penguins %>% filter(!is.na(sex), !is.na(species)) ``` We'll focus on 3 variables with the goal of predicting `flipper_len`: - `flipper_len` = the length of the penguin's flippers (arms) in mm - `species` = Adelie, Chinstrap, or Gentoo - `sex` = female, male And build 2 models of `flipper_len`: ```{r} # Model of flipper_len by sex flipper_model_1 <- lm(flipper_len ~ sex, data = penguins) flipper_model_2 <- lm(flipper_len ~ species, data = penguins) ``` Just as with models that use *quantitative* predictors, it's important to *evaluate* our (eventual) models of `flipper_len` by `species` and `sex`: are they correct? strong? fair? Let's start with strength. How strong are these models? What's the best predictor? Let's explore. a. Based on the boxplots below, what is the stronger predictor of `flipper_len`: `sex` or `species`? ```{r} penguins %>% ggplot(aes(y = flipper_len, x = sex)) + geom_boxplot() ``` ```{r} penguins %>% ggplot(aes(y = flipper_len, x = species)) + geom_boxplot() ``` b. We used 2 approaches to measuring the strength of the simple linear regression models in previous activities: correlation and R-squared. Unfortunately, we get errors when we try to calculate the correlation between `flipper_len` and `sex`, and between `flipper_len` and `species`. Why? ```{r} penguins %>% summarize(cor(sex, flipper_len)) penguins %>% summarize(cor(species, flipper_len)) ``` c. Luckily, R-squared works no matter whether a predictor is quantitative or categorical! Interpret and compare the R-squared values for our separate models of `flipper_len` by `sex` and `species`. Which is the stronger predictor? Does this match your answer in part a? ```{r} # Model of flipper_len by sex summary(flipper_model_1) # Model of flipper_len by species summary(flipper_model_2) ``` ## Exercise 13: Evaluating model correctness Let's just focus on the stronger of our 2 models, that of `flipper_len` by `species` (`flipper_model_2`), and ask: is it correct (not wrong)? Recall that when our predictor was *quantitative*, residual plots provided some insight. Check out the residual plots below. They look a little goofy! Explain the goofiness (what's happening here) and describe what you learn about the model's "correctness". Is the model "correct"? ```{r} # Residual plot flipper_model_2 %>% ggplot(aes(x = .fitted, y = .resid)) + geom_point() + geom_hline(yintercept = 0) ``` ```{r} # Residual plot using boxes! flipper_model_2 %>% ggplot(aes(x = .fitted, y = .resid, group = .fitted)) + geom_boxplot() + geom_hline(yintercept = 0) ``` ## Exercise 14: Really understand how the coefficients work Next, consider the model of `flipper_len` by `island`. The average `flipper_len` of penguins on each `island` is calculated below: ```{r} penguins %>% group_by(island) %>% summarize(mean(flipper_len)) ``` a. Using just the above averages and your understanding of categorical predictors, fill in the model coefficients and indicator variables below. Do NOT use `lm()` yet!!! E[`flipper_len` | `island`] = ___ +/- ___ island??? +/- ___ island??? b. Check your work to Part a. ```{r} flipper_model_3 <- lm(flipper_len ~ island, penguins) summary(flipper_model_3) ``` \ \ \ \ # Solutions ## Exercise 1: Get to know the data - A case represents a single diamond. - The distribution of price is right skewed with considerable high outliers. The right skew is evidenced by the mean price (\$3932) being much higher than the median price (\$2401). - Most diamonds in this data are of Good cut or better. Ideal cut diamonds are the most common with each succesive grade being the next most common. ```{r eval = TRUE} dim(diamonds) head(diamonds) # Visualize price (outcome variable) ggplot(diamonds, aes(x = price)) + geom_histogram() ggplot(diamonds, aes(y = price)) + geom_boxplot() diamonds %>% summarize(mean(price), median(price), sd(price)) # Visualize cut (predictor variable) ggplot(diamonds, aes(x = cut)) + geom_bar() diamonds %>% count(cut) ``` ## Exercise 2: Visualizations Start by visualizing this relationship of interest, that between `price` and `cut`. a. We just don't see anything clearly on a scatterplot. With the small number of unique values of the predictor variable, all of the points are bunched up on each other. ```{r eval = TRUE} # Try a scatterplot ggplot(diamonds, aes(y = price, x = cut)) + geom_point() ``` b. Separately run each chunk below, with two plots. Comment (`#`) on what changes in the code / output. ```{r eval = TRUE} # Univariate boxplot ggplot(diamonds, aes(y = price)) + geom_boxplot() ``` ```{r eval = TRUE} # Separate boxes by category ggplot(diamonds, aes(y = price, x = cut)) + geom_boxplot() ``` ```{r eval = TRUE} # Univariate density plot ggplot(diamonds, aes(x = price)) + geom_density() ``` ```{r eval = TRUE} # Separate density plots by category ggplot(diamonds, aes(x = price, color = cut)) + geom_density() ``` ```{r eval = TRUE} # Univariate histogram ggplot(diamonds, aes(x = price)) + geom_histogram() ``` ```{r eval = TRUE} # Separate histograms by category ggplot(diamonds, aes(x = price)) + geom_histogram() + facet_wrap(~ cut) ``` c. The relationship between price and cut seems to be opposite what we would expect. The diamonds with the best cut (Ideal) have the lowest average price, and the ones with the worst cut (Fair) are woth the most. Maybe something else is different between the diamonds with the best and worst cuts...size maybe? ## Exercise 3: Numerical summaries Let's follow up our plots with some numerical summaries. a. Mean `price` across all diamonds: ```{r eval = TRUE} diamonds %>% summarize(mean(price)) ``` b. Mean `price` for each type of `cut`: ```{r eval = TRUE} diamonds %>% group_by(cut) %>% summarize(mean(price)) ``` c. Group means should reflect what you see in the plots (easiest to see in the boxplots) d. Create our new `cutFair` variable: ```{r eval = TRUE} diamonds <- diamonds %>% mutate(cutFair=ifelse(cut == "Fair", 1, 0)) ``` e. Calculate the group means based on this new variable ```{r eval = TRUE} diamonds %>% group_by(cutFair) %>% summarize(mean(price)) ``` ## Exercise 4: Modeling trend using a categorical predictor with exactly 2 categories ```{r eval = TRUE} # Construct the model diamond_mod0 <- lm(price ~ cutFair, data = diamonds) # Summarize the model coef(summary(diamond_mod0)) ``` The intercept is the expected value (mean) of the price for all diamonds with a cut quality that isn't Fair (Good, Very Good, Premium, or Ideal, i.e. when `cutFair` = 0)--the same as we saw in exercise 3e. c. When we add the intercept and coefficient for `cutFair`, we get `3919.69 + 439.06 = 4358.75`--this is the mean price for all diamonds with a Fair cut quality that we saw in exercise 3e! Therefore, the coefficient of `cutFair` (439.06) is interpreted as the difference between the mean value of diamonds with a Fair cut quality and the mean value of diamonds with a higher cut quality. ## Exercise 5: Modeling trend using a categorical predictor with >2 categories ```{r eval = TRUE} # Construct the model diamond_mod <- lm(price ~ cut, data = diamonds) # Summarize the model coef(summary(diamond_mod)) ``` a. We are seeing 4 coefficient estimates because each category is being assigned to a separate indicator variable--`cutGood` = 1 when `cut == "Good"` and 0 otherwise, `cutVery Good` = 1 when `cut == "Very Good" and 0 otherwise, and so on. b. E[price | cut] = 4358.7578 - 429.8933 cutGood - 376.9979 cutVery Good + 225.4999 cutPremium - 901.2158 cutIdeal ## Exercise 6: Making sense of the model a. Expected/typical price for diamonds of *Good* cut: E[price | cut] = 4358.7578 - 429.8933 * 1 - 376.9979 * 0 + 225.4999 * 0 - 901.2158 * 0 = 4358.7578 - 429.8933 = \$3928.865 ```{r eval = TRUE} predict(diamond_mod, newdata = data.frame(cut = "Good")) ``` b. Expected/typical price for diamonds of *Fair* cut: E[price | cut] = 4358.7578 - 429.8933 * 0 - 376.9979 * 0 + 225.4999 * 0 - 901.2158 * 0 = \$4358.7578 ```{r eval = TRUE} predict(diamond_mod, newdata = data.frame(cut = "Fair")) ``` c. These come from our group mean calculations in Exercise 3b! The predicted value for diamonds of *Fair* cut is also the same as what we obtained using the SLR model in exercise 4 with only a single `cutFair` indicator variable. ## Exercise 7: Interpreting coefficients Recall that our model formula is *not a formula for a line*. Thus we can't interpret the coefficients as "slopes" as we have before. Taking this into account and reflecting upon your calculations above... a. The average price of a Fair cut diamonds is \$4358.7578. b. - Interpretation of `cutGood` coefficient: On average, Good cut diamonds are worth \$429.89 less than Fair cut diamonds. - Interpretation of `cutVery Good` coefficient: On average, Very Good cut diamonds are worth \$377.00 less than Fair cut diamonds. ## Exercise 8: Modeling choices (CHALLENGE) Why do we fit this model in this way (using 4 indicator variables `cutGood`, `cutVery Good`, `cutPremium`, `cutIdeal`)? Instead, suppose that we created a single variable `cutCat` that gave each category a numerical value: 0 for Fair, 1 for Good, 2 for Very Good, 3 for Premium, and 4 for Ideal. - If we used 0-4 instead of creating indicator variables, we would be constraining the change from 0 to 1, from 1 to 2, etc. to always be of the same magnitude. That is, a 1 unit change in the cut variable would always have the same change in price in our model. - Using separate indicator variables allows the difference between subsequent categories to be different, which allows our model to be a bit more nuanced. It is possible to take nuance too far though. For example, in our previous investigations of bikeshare data, we modeled ridership versus temperature. We treated temperature as a quantitative predictor. Imagine if we had created an indicator variable for each unique temperature in the data---that would be so many variables! Having so many variables creates a very complex model which can be hard to make sense of. (These ideas are addressed further in STAT 253: Statistical Machine Learning!) ## Exercise 9: The least squares criterion ```{r eval = TRUE} # Observed price = 326 diamonds %>% select(price, cut) %>% head(1) # Predicted price = 3457.542 predict(diamond_mod, newdata = data.frame(cut = "Ideal")) # Residual 326 - 3457.542 ``` ## Exercise 10: Diamond color Consider modeling `price` by `color`. - The best color diamonds are J, and worst are D. We would expect D diamonds to have the lowest price and increase steadily as we get to J. This is in fact what we see in the boxplots. ```{r eval = TRUE} ggplot(diamonds, aes(x = color, y = price)) + geom_boxplot() diamonds %>% group_by(color) %>% summarize(mean(price)) ``` - We fit a linear model and obtain the model formula: E[price | color] = 3169.95 - 93.20 colorE + 554.93 colorF + 829.18 colorG + 1316.72 colorH + 1921.92 colorI + 2153.86 colorJ ```{r eval = TRUE} diamond_mod2 <- lm(price ~ color, data = diamonds) coef(summary(diamond_mod2)) ``` - Color D is the reference level because we don't see its indicator variable in the model output. - Interpretation of the intercept: Diamonds with D color cost \$3169.95 on average. - Interpretation of the `colorE` coefficient: Diamonds with E color cost \$93.20 **less** than D color diamonds on average. - Interpretation of the `colorF` coefficient: Diamonds with F color cost \$554.93 **more** than D color diamonds on average. ## Exercise 11: Diamond clarity We see the unexpected result that diamonds of better clarity (VS1 and higher) have lower average prices. In fact the best clarity diamonds (VVS1 and IF) have the *lowest* average prices. What might be going on? What if the most clear diamonds were also quite small... ```{r eval = TRUE} ggplot(diamonds, aes(x = clarity, y = price)) + geom_boxplot() diamonds %>% group_by(clarity) %>% summarize(mean(price)) diamond_mod3 <- lm(price ~ clarity, data = diamonds) coef(summary(diamond_mod3)) ``` ## Exercise 12: Evaluating model strength ```{r eval = TRUE} flipper_model_1 <- lm(flipper_len ~ sex, data = penguins) flipper_model_2 <- lm(flipper_len ~ species, data = penguins) ``` ### Part a `species` appears to be the stronger predictor - the `flipper_len` values are more distinct between species (there's less overlap in the boxes) ```{r eval = TRUE} data(penguins) penguins <- penguins %>% filter(!is.na(sex), !is.na(species)) penguins %>% ggplot(aes(y = flipper_len, x = sex)) + geom_boxplot() penguins %>% ggplot(aes(y = flipper_len, x = species)) + geom_boxplot() ``` ### Part b Correlation only works when both variables in the pair are *quantitative*, but `sex` and `species` are *categorical*. ### Part c Using `sex`: R-squared = 0.06511 Using `species`: R-squared = 0.7782. Thus `species` is the stronger predictor of `flipper_len`. It explains roughly 78% of the variability in flipper length from penguin to penguin. ```{r eval = TRUE} summary(flipper_model_1) summary(flipper_model_2) ``` ## Exercise 13: Evaluating model correctness Since there are only 3 different `species`, there are only 3 possible predictions of `flipper_len` in this model (one per species). Thus each of the 3 groups on the x-axis corresponds to a species. Within each species, the observed `flipper_len` deviates from the species-based prediction, thus so too do the residuals. That's why the points in each group have different y coordinates. Overall, our model doesn't seem too wrong! No matter the x value (species!), the residuals are normally balanced above and below 0 with roughly equal variance in each species. ```{r eval = TRUE} # Residual plot flipper_model_2 %>% ggplot(aes(x = .fitted, y = .resid)) + geom_point() + geom_hline(yintercept = 0) # Residual plot using boxes! flipper_model_2 %>% ggplot(aes(x = .fitted, y = .resid, group = .fitted)) + geom_boxplot() + geom_hline(yintercept = 0) ``` ## Exercise 14: Really understand how the coefficients work ```{r eval = TRUE} penguins %>% group_by(island) %>% summarize(mean(flipper_len)) ``` a. E[`flipper_len` | `island`] = 209.56 - 16.37 islandDream - 18.03 islandTorgersen Why? - The island terms are `islandDream` and `islandTorgersen` since Biscoe is the reference. - The intercept is the average flipper length for the reference, Biscoe. - The `islandDream` coefficient reflects the fact that the average flipper length on Dream is 16.37mm *lower* than that on Biscoe (193.1870 - 209.5583) - The `islandTorgersen` coefficient reflects the fact that the average flipper length on Dream is 18.03mm *lower* than that on Biscoe (191.5319 - 209.5583) b. ```{r eval = TRUE} flipper_model_3 <- lm(flipper_len ~ island, penguins) summary(flipper_model_3) ```