# Chapter 5 Basic Regression

Now that we are equipped with data visualization skills from Chapter 2, an understanding of the "tidy" data format from Chapter 4, and data wrangling skills from Chapter 3, we now proceed with data modeling. The fundamental premise of data modeling is *to make explicit the relationship* between:

- an outcome variable \(y\), also called a dependent variable and
- an explanatory/predictor variable \(x\), also called an independent variable or covariate.

Another way to state this is using mathematical terminology: we will model the outcome variable \(y\) *as a function* of the explanatory/predictor variable \(x\). Why do we have two different labels, explanatory and predictor, for the variable \(x\)? That's because roughly speaking data modeling can be used for two purposes:

**Modeling for prediction**: You want to predict an outcome variable \(y\) based on the information contained in a set of predictor variables. You don't care so much about understanding how all the variables relate and interact, but so long as you can make good predictions about \(y\), you're fine. For example, if we know many individuals' risk factors for lung cancer, such as smoking habits and age, can we predict whether or not they will develop lung cancer? Here we wouldn't care so much about distinguishing the degree to which the different risk factors contribute to lung cancer, but instead only on whether or not they could be put together to make reliable predictions.**Modeling for explanation**: You want to explicitly describe the relationship between an outcome variable \(y\) and a set of explanatory variables, determine the significance of any found relationships, and have measures summarizing these. Continuing our example from above, we would now be interested in describing the individual effects of the different risk factors and quantifying the magnitude of these effects. One reason could be to design an intervention to reduce lung cancer cases in a population, such as targeting smokers of a specific age group with an advertisement for smoking cessation programs. In this book, we'll focus more on this latter purpose.

Data modeling is used in a wide variety of fields, including statistical inference, causal inference, artificial intelligence, and machine learning. There are many techniques for data modeling, such as tree-based models, neural networks and deep learning, and supervised learning. In this chapter, we'll focus on one particular technique: *linear regression*, one of the most commonly-used and easy-to-understand approaches to modeling. Recall our discussion in Subsection 1.4.3 on numerical and categorical variables. Linear regression involves:

- an outcome variable \(y\) that is
*numerical*and - explanatory variables \(x_i\) (e.g. \(x_1, x_2, ...\)) that are either
*numerical*or*categorical*.

With linear regression there is always only one numerical outcome variable \(y\) but we have choices on both the number and the type of explanatory variables to use. We're going to cover the following regression scenarios:

- In this current chapter on basic regression, we'll always have only one explanatory variable.
- In the next chapter, Chapter 6 on
*multiple regression*, we'll have more than one explanatory variable:

We'll study all four of these regression scenarios using real data, all easily accessible via R packages!

### Needed packages

Let's now load all the packages needed for this chapter (this assumes you've already installed them). In this chapter we introduce some new packages:

- The
`tidyverse`

"umbrella" package. Recall from our discussion in Section 4.4.1 that loading the`tidyverse`

package by running`library(tidyverse)`

loads the following commonly used data science packages all at once:`ggplot2`

for data visualization`dplyr`

for data wrangling`tidyr`

for converting data to "tidy" format`readr`

for importing spreadsheet data into R- As well as the more advanced
`purrr`

,`tibble`

,`stringr`

, and`forcats`

packages

- The
`skimr`

(Waring et al. 2021) package, which provides a simple-to-use function to quickly compute a wide array of commonly-used summary statistics. - The
`gapminder`

package, which provides excerpts of data available from Gapminder.org - The
`moderndive`

package, which accompanies this book and includes data sets we will analyze

If needed, read Section 1.3 for information on how to install and load R packages.

```
library(tidyverse)
library(skimr)
library(gapminder)
library(moderndive)
```

## 5.1 One numerical explanatory variable

Why do some professors and instructors at universities and colleges get high teaching evaluations from students while others don't? What factors can explain these differences? Are there biases? These are questions that are of interest to university/college administrators, as teaching evaluations are among the many criteria considered in determining which professors and instructors should get promotions. Researchers at the University of Texas in Austin, Texas (UT Austin) tried to answer this question: what factors can explain differences in instructor's teaching evaluation scores? To this end, they collected information on \(n = 463\) instructors. A full description of the study can be found at openintro.org.

We'll keep things simple for now and try to explain differences in instructor evaluation scores as a function of one numerical variable: their "beauty score." The specifics on how this score was calculated will be described shortly.

Could it be that instructors with higher beauty scores also have higher teaching evaluations? Could it be instead that instructors with higher beauty scores tend to have lower teaching evaluations? Or could it be there is no relationship between beauty score and teaching evaluations?

We'll achieve ways to address these questions by modeling the relationship between these two variables with a particular kind of linear regression called *simple linear regression*. Simple linear regression is the most basic form of linear regression. With it we have

- A numerical outcome variable \(y\). In this case, an instructor's teaching score.
- A single numerical explanatory variable \(x\). In this case, an instructor's beauty score.

### 5.1.1 Exploratory data analysis

A crucial step before doing any kind of modeling or analysis is performing an *exploratory data analysis*, or EDA, of all our data. Exploratory data analysis can give you a sense of the distribution of the data and whether there are outliers and/or missing values. Most importantly, it can inform how to build your model. There are many approaches to exploratory data analysis; here are three:

- Most fundamentally: just looking at the raw values, in a spreadsheet for example. While this may seem trivial, many people ignore this crucial step!
- Computing summary statistics like means, medians, and standard deviations.
- Creating data visualizations.

Let's load the `evals`

data (which is built into the `moderndive`

package), `select`

only a subset of the variables, and look at the raw values. Recall you can look at the raw values by running `View()`

in the console in RStudio to pop-up the spreadsheet viewer with the data frame of interest as the argument to `View()`

. Here, however, we present only a snapshot of five randomly chosen rows:

```
evals_ch5 <- evals %>%
select(score, bty_avg, age)
```

```
evals_ch5 %>%
sample_n(5)
```

score | bty_avg | age |
---|---|---|

3.7 | 3.00 | 62 |

4.7 | 4.33 | 46 |

4.8 | 5.50 | 62 |

2.8 | 2.00 | 62 |

4.0 | 2.33 | 64 |

While a full description of each of these variables can be found at openintro.org, let's summarize what each of these variables represents.

`score`

: Numerical variable of the average teaching score based on students' evaluations between 1 and 5. This is the outcome variable \(y\) of interest.`bty_avg`

: Numerical variable of average "beauty" rating based on a panel of 6 students' scores between 1 and 10. This is the numerical explanatory variable \(x\) of interest. Here 1 corresponds to a low beauty rating and 10 to a high beauty rating.`age`

: A numerical variable of age in years as an integer value.

An alternative way to look at the raw data values is by choosing a random sample of the rows in `evals_ch5`

by piping it into the `sample_n()`

function from the `dplyr`

package. Here we set the `size`

argument to be `5`

, indicating that we want a random sample of 5 rows. We display the results in Table 5.2. Note that due to the random nature of the sampling, you will likely end up with a different subset of 5 rows.

```
evals_ch5 %>%
sample_n(size = 5)
```

score | bty_avg | age |
---|---|---|

4.4 | 5.67 | 57 |

3.1 | 7.00 | 33 |

4.1 | 4.17 | 45 |

5.0 | 4.33 | 46 |

4.8 | 4.83 | 52 |

Now that we've looked at the raw values in our `evals_ch5`

data frame and got a preliminary sense of the data, let's move on to the next common step in an exploratory data analysis: computing summary statistics. Let's start by computing the mean and median of our numerical outcome variable `score`

and our numerical explanatory variable "beauty" score denoted as `bty_avg`

. We'll do this by using the `summarize()`

function from `dplyr`

along with the `mean()`

and `median()`

summary functions we saw in Section 3.3.

```
evals_ch5 %>%
summarize(mean_bty_avg = mean(bty_avg), mean_score = mean(score),
median_bty_avg = median(bty_avg), median_score = median(score))
```

```
# A tibble: 1 × 4
mean_bty_avg mean_score median_bty_avg median_score
<dbl> <dbl> <dbl> <dbl>
1 4.42 4.17 4.33 4.3
```

However, what if we want other summary statistics as well, such as the standard deviation (a measure of spread), the minimum and maximum values, and various percentiles? Typing out all these summary statistic functions in `summarize()`

would be long and tedious. Instead, let's use the convenient `skim()`

function from the `skimr`

package. This function takes in a data frame, "skims" it, and returns commonly used summary statistics. Let's take our `evals_ch5`

data frame, `select()`

only the outcome and explanatory variables teaching `score`

and `bty_avg`

, and pipe them into the `skim()`

function:

```
evals_ch5 %>%
select(score, bty_avg) %>%
skim()
```

```
Skim summary statistics
n obs: 463
n variables: 2
── Variable type:numeric ───────────────────────────────────────────────────────
variable missing complete n mean sd p0 p25 p50 p75 p100
bty_avg 0 463 463 4.42 1.53 1.67 3.17 4.33 5.5 8.17
score 0 463 463 4.17 0.54 2.3 3.8 4.3 4.6 5
```

(Note that for formatting purposes, the inline histogram that is usually printed with `skim()`

has been removed. This can be done by running `skim_with(numeric = list(hist = NULL))`

prior to using the `skim()`

function as well.)

For our two numerical variables teaching `score`

and "beauty" score `bty_avg`

it returns:

`missing`

: the number of missing values`complete`

: the number of non-missing or complete values`n`

: the total number of values`mean`

: the average`sd`

: the standard deviation`p0`

: the 0^{th}percentile: the value at which 0% of observations are smaller than it (the*minimum*value)`p25`

: the 25^{th}percentile: the value at which 25% of observations are smaller than it (the*1*)^{st}quartile`p50`

: the 50^{th}percentile: the value at which 50% of observations are smaller than it (the*2*quartile and more commonly called the^{nd}*median*)`p75`

: the 75^{th}percentile: the value at which 75% of observations are smaller than it (the*3*)^{rd}quartile`p100`

: the 100^{th}percentile: the value at which 100% of observations are smaller than it (the*maximum*value)

Looking at this output, we get an idea of how the values of both variables distribute. For example, the mean teaching score was 4.17 out of 5 whereas the mean "beauty" score was 4.42 out of 10. Furthermore, the middle 50% of teaching scores were between 3.80 and 4.6 (the first and third quartiles) whereas the middle 50% of "beauty" scores were between 3.17 and 5.5 out of 10.

The `skim()`

function only returns what are known as *univariate* summary statistics: functions that take a single variable and return some numerical summary of that variable. However, there also exist *bivariate* summary statistics: functions that take in two variables and return some summary of those two variables. In particular, when the two variables are numerical, we can compute the *correlation coefficient*. Generally speaking, *coefficients* are quantitative expressions of a specific phenomenon. A *correlation coefficient* is a quantitative expression of the *strength of the linear relationship between two numerical variables*. Its value ranges between -1 and 1 where:

- -1 indicates a perfect
*negative relationship*: As the value of one variable goes up, the value of the other variable tends to go down following along a straight line. - 0 indicates no relationship: The values of both variables go up/down independently of each other.
- +1 indicates a perfect
*positive relationship*: As the value of one variable goes up, the value of the other variable tends to go up as well in a linear fashion.

Figure 5.1 gives examples of different correlation coefficient values for hypothetical numerical variables \(x\) and \(y\). We see that while for a correlation coefficient of -0.75 there is still a negative relationship between \(x\) and \(y\), it is not as strong as the negative relationship between \(x\) and \(y\) when the correlation coefficient is -1.

The correlation coefficient is computed using the `cor()`

function, where the inputs to the function are the two numerical variables for which we want to quantify the strength of the linear relationship.

```
evals_ch5 %>%
summarise(correlation = cor(score, bty_avg))
```

```
# A tibble: 1 × 1
correlation
<dbl>
1 0.187
```

You can also use the `cor()`

function directly instead of using it inside `summarise`

, but you will need to use the `$`

syntax to access the specific variables within a data frame (See Subsection 1.4.3):

`cor(x = evals_ch5$bty_avg, y = evals_ch5$score)`

`[1] 0.187`

In our case, the correlation coefficient of 0.187 indicates that the relationship between teaching evaluation score and beauty average is "weakly positive." There is a certain amount of subjectivity in interpreting correlation coefficients, especially those that aren't close to -1, 0, and 1. For help developing such intuition and more discussion on the correlation coefficient see Subsection 5.3.1 below.

Let's now proceed by visualizing this data. Since both the `score`

and `bty_avg`

variables are numerical, a scatterplot is an appropriate graph to visualize this data. Let's do this using `geom_point()`

and set informative axes labels and title and display the result in Figure 5.2.

```
ggplot(evals_ch5, aes(x = bty_avg, y = score)) +
geom_point() +
labs(x = "Beauty Score", y = "Teaching Score",
title = "Relationship of teaching and beauty scores")
```

Observe the following:

- Most "beauty" scores lie between 2 and 8.
- Most teaching scores lie between 3 and 5.
- Recall our earlier computation of the correlation coefficient, which describes the strength of the linear relationship between two numerical variables. Looking at Figure 5.3, it is not immediately apparent that these two variables are positively related. This is to be expected given the positive, but rather weak (close to 0), correlation coefficient of 0.187.

Before we continue, we bring to light an important fact about this dataset: it suffers from *overplotting*. Recall from the data visualization Subsection 2.3.2 that overplotting occurs when several points are stacked directly on top of each other thereby obscuring the number of points. For example, let's focus on the 6 points in the top-right of the plot with a beauty score of around 8 out of 10: are there truly only 6 points, or are there many more just stacked on top of each other? You can think of these as *ties*. Let's break up these ties with a little random "jitter" added to the points in Figure 5.3.

Jittering adds a little random bump to each of the points to break up these ties: just enough so you can distinguish them, but not so much that the plot is overly altered. Furthermore, jittering is strictly a visualization tool; it does not alter the original values in the dataset.

Let's compare side-by-side the regular scatterplot in Figure 5.2 with the jittered scatterplot in Figure 5.3 in Figure 5.4.

We make several further observations:

- Focusing our attention on the top-right of the plot again, as noted earlier where there seemed to only be 6 points in the regular scatterplot, we see there were in fact really 9 as seen in the jittered scatterplot.
- A further interesting trend is that the jittering revealed a large number of instructors with beauty scores of between 3 and 4.5, towards the lower end of the beauty scale.

To keep things simple in this chapter, we'll present regular scatterplots rather than the jittered scatterplots, though we'll keep the overplotting in mind whenever looking at such plots. Going back to scatterplot in Figure 5.2, let's improve on it by adding a "regression line" in Figure 5.5. This is easily done by adding a new layer to the `ggplot`

code that created Figure 5.3: `+ geom_smooth(method = "lm")`

. A regression line is a "best fitting" line in that of all possible lines you could draw on this plot, it is "best" in terms of some mathematical criteria. We discuss the criteria for "best" in Subsection 5.3.3 below, but we suggest you read this only after covering the concept of a *residual* coming up in Subsection 5.1.3.

```
ggplot(evals_ch5, aes(x = bty_avg, y = score)) +
geom_point() +
labs(x = "Beauty Score", y = "Teaching Score",
title = "Relationship of teaching and beauty scores") +
geom_smooth(method = "lm")
```

``geom_smooth()` using formula 'y ~ x'`

When viewed on this plot, the regression line is a visual summary of the relationship between two numerical variables, in our case the outcome variable `score`

and the explanatory variable `bty_avg`

. The positive slope of the blue line is consistent with our observed correlation coefficient of 0.187 suggesting that there is a positive relationship between `score`

and `bty_avg`

. We'll see later however that while the correlation coefficient is not equal to the slope of this line, they always have the same sign: positive or negative.

What are the grey bands surrounding the blue line? These are *standard error* bands, which can be thought of as error/uncertainty bands. Let's skip this idea for now and suppress these grey bars by adding the argument `se = FALSE`

to `geom_smooth(method = "lm")`

. We'll introduce standard errors in Chapter 9 on sampling distributions, use them for constructing *confidence intervals* and conducting *hypothesis tests* in Chapters 10 and 12, and consider them when we revisit regression in Chapter 12.

```
ggplot(evals_ch5, aes(x = bty_avg, y = score)) +
geom_point() +
labs(x = "Beauty Score", y = "Teaching Score",
title = "Relationship of teaching and beauty scores") +
geom_smooth(method = "lm", se = FALSE)
```

``geom_smooth()` using formula 'y ~ x'`

*Learning check*

**(LC5.1)** Conduct a new exploratory data analysis with the same outcome variable \(y\) being `score`

but with `age`

as the new explanatory variable \(x\). Remember, this involves three things:

- Looking at the raw values.
- Computing summary statistics of the variables of interest.
- Creating informative visualizations.

What can you say about the relationship between age and teaching scores based on this exploration?

### 5.1.2 Simple linear regression

You may recall from secondary school / high school algebra, in general, the equation of a line is \(y = a + bx\), which is defined by two coefficients. Recall we defined this earlier as "quantitative expressions of a specific property of a phenomenon." These two coefficients are:

- the intercept coefficient \(a\), or the value of \(y\) when \(x = 0\), and
- the slope coefficient \(b\), or the increase in \(y\) for every increase of one in \(x\).

However, when defining a line specifically for regression, like the blue regression line in Figure 5.6, we use slightly different notation: the equation of the regression line is \(\widehat{y} = b_0 + b_1 \cdot x\) where

- the intercept coefficient is \(b_0\), or the value of \(\widehat{y}\) when \(x=0\), and
- the slope coefficient \(b_1\), or the increase in \(\widehat{y}\) for every increase of one in \(x\).

Why do we put a "hat" on top of the \(y\)? It's a form of notation commonly used in regression, which we'll introduce in the next Subsection 5.1.3 when we discuss *fitted values*. For now, let's ignore the hat and treat the equation of the line as you would from secondary school / high school algebra recognizing the slope and the intercept. We know looking at Figure 5.6 that the slope coefficient corresponding to `bty_avg`

should be positive. Why? Because as `bty_avg`

increases, professors tend to roughly have higher teaching evaluation `scores`

. However, what are the specific values of the intercept and slope coefficients? Let's not worry about computing these by hand, but instead let the computer do the work for us. Specifically, let's use R!

Let's get the value of the intercept and slope coefficients by outputting something called the *linear regression table*. We will fit the linear regression model to the `data`

using the `lm()`

function and save this to `score_model`

. `lm`

stands for "linear model." When we say "fit", we are saying find the best fitting line to this data.

The `lm()`

function that "fits" the linear regression model is typically used as `lm(y ~ x, data = data_frame_name)`

where:

`y`

is the outcome variable, followed by a tilde (`~`

). This is likely the key to the left of "1" on your keyboard. In our case,`y`

is set to`score`

.`x`

is the explanatory variable. In our case,`x`

is set to`bty_avg`

. We call the combination`y ~ x`

a*model formula*.`data_frame_name`

is the name of the data frame that contains the variables`y`

and`x`

. In our case,`data_frame_name`

is the`evals_ch5`

data frame.

```
score_model <- lm(score ~ bty_avg, data = evals_ch5)
score_model
```

```
Call:
lm(formula = score ~ bty_avg, data = evals_ch5)
Coefficients:
(Intercept) bty_avg
3.8803 0.0666
```

This output is telling us that the `Intercept`

coefficient \(b_0\) of the regression line is 3.8803, and the slope coefficient for `by_avg`

is 0.0666. Therefore the blue regression line in Figure 5.6 is

\[\widehat{\text{score}} = b_0 + b_{\text{bty avg}} \cdot\text{bty avg} = 3.8803 + 0.0666\cdot\text{ bty avg}\]

where

- The intercept coefficient \(b_0 = 3.8803\) means for instructors that had a hypothetical beauty score of 0, we would expect them to have on average a teaching score of 3.8803. In this case however, while the intercept has a mathematical interpretation when defining the regression line, there is no
*practical*interpretation since`score`

is an average of a panel of 6 students' ratings from 1 to 10, a`bty_avg`

of 0 would be impossible. Furthermore, no instructors had a beauty score anywhere near 0 in this data. Of more interest is the slope coefficient associated with

`bty_avg`

: \(b_{\text{bty avg}} = +0.0666\). This is a numerical quantity that summarizes the relationship between the outcome and explanatory variables. Note that the sign is positive, suggesting a positive relationship between beauty scores and teaching scores, meaning as beauty scores go up, so also do teaching scores go up. The slope's precise interpretation is:For every increase of 1 unit in

`bty_avg`

, there is an*associated*increase of,*on average*, 0.0666 units of`score`

.

Such interpretations need be carefully worded:

- We only stated that there is an
*associated*increase, and not necessarily a*causal*increase. For example, perhaps it's not that beauty directly affects teaching scores, but instead younger instructors tend to be percieved as better teachers (perhaps because they are more energetic or use more modern techniques), and younger instructors are also perceived to be more beautiful. Avoiding such reasoning can be summarized by the adage "correlation is not necessarily causation." In other words, just because two variables are correlated, it doesn't mean one directly causes the other. We discuss these ideas more in Subsection 5.3.2 and in Chapter 7.

- We say that this associated increase is
*on average*0.0666 units of teaching`score`

and not that the associated increase is*exactly*0.0666 units of`score`

across all values of`bty_avg`

. This is because the slope is the average increase across all points as shown by the regression line in Figure 5.6.

Now that we've learned how to compute the equation for the blue regression line in Figure 5.6 and interpreted all its terms, let's take our modeling one step further. This time after fitting the model using the `lm()`

, let's get something called the *regression table* using the `summary()`

function:

```
# Fit regression model:
score_model <- lm(score ~ bty_avg, data = evals_ch5)
# Get regression results:
summary(score_model)
```

```
Call:
lm(formula = score ~ bty_avg, data = evals_ch5)
Residuals:
Min 1Q Median 3Q Max
-1.925 -0.369 0.142 0.398 0.931
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.8803 0.0761 50.96 < 0.0000000000000002 ***
bty_avg 0.0666 0.0163 4.09 0.000051 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.535 on 461 degrees of freedom
Multiple R-squared: 0.035, Adjusted R-squared: 0.0329
F-statistic: 16.7 on 1 and 461 DF, p-value: 0.0000508
```

Note how we took the output of the model fit saved in `score_model`

and used it as an input to the subsequent `summary()`

function. The raw output of the `summary()`

function above gives lots of information about the regression model that we won't cover in this introductory course (e.g., Multiple R-squared, F-statistic, etc.). We will only consider the "Coefficients" section of the output. We can print these relevant results only by accessing the `coefficients`

object stored in the summary results.

`summary(score_model)$coefficients`

Estimate | Std. Error | t value | Pr(>|t|) | |
---|---|---|---|---|

(Intercept) | 3.880 | 0.076 | 50.96 | 0 |

bty_avg | 0.067 | 0.016 | 4.09 | 0 |

For now since we are only using regression as an exploratory data analysis tool, we will only focus on the "Estimate" column that contains the estimates of the intercept and slope for the best fit line for our data. The remaining three columns refer to statistical concepts known as standard errors, t-statistics, and p-values, which we will cover in Parts III and IV of the book when we talk about statistical theory and inference.

*Learning check*

**(LC5.2)** Fit a new simple linear regression using `lm(score ~ age, data = evals_ch5)`

where `age`

is the new explanatory variable \(x\). Get information about the "best-fitting" line from the regression table by applying the `summary()`

function. How do the regression results match up with the results from your exploratory data analysis above?

### 5.1.3 Observed/fitted values and residuals

We just saw how to get the value of the intercept and the slope of the regression line from the regression table generated by `summary()`

. Now instead, say we want information on individual points. In this case, we focus on one of the \(n = 463\) instructors in this dataset, corresponding to a single row of `evals_ch5`

.

For example, say we are interested in the 21st instructor in this dataset:

score | bty_avg | age |
---|---|---|

4.9 | 7.33 | 31 |

What is the value on the blue line corresponding to this instructor's `bty_avg`

of 7.333? In Figure 5.7 we mark three values in particular corresponding to this instructor.

- Red circle: This is the
*observed value*\(y\) = 4.9 and corresponds to this instructor's actual teaching score. - Red square: This is the
*fitted value*\(\widehat{y}\) and corresponds to the value on the regression line for \(x\) = 7.333. This value is computed using the intercept and slope in the regression table above: \[\widehat{y} = b_0 + b_1 \cdot x = 3.88 + 0.067 * 7.333 = 4.369\] - Blue arrow: The length of this arrow is the
*residual*and is computed by subtracting the fitted value \(\widehat{y}\) from the observed value \(y\). The residual can be thought of as the error or "lack of fit" of the regression line. In the case of this instructor, it is \(y - \widehat{y}\) = 4.9 - 4.369 = 0.531. In other words, the model was off by 0.531 teaching score units for this instructor.

``geom_smooth()` using formula 'y ~ x'`

What if we want both

- the fitted value \(\widehat{y} = b_0 + b_1 \cdot x\) and
- the residual \(y - \widehat{y}\)

for not only the 21st instructor but for all 463 instructors in the study? Recall that each instructor corresponds to one of the 463 rows in the `evals_ch5`

data frame and also one of the 463 points in the regression plot in Figure 5.6.

We could repeat the above calculations by hand 463 times, but that would be tedious and time consuming. Instead, let's use our data wrangling tools from Chapter 3 and the functions `fitted()`

and `residuals()`

to create a data frame with these values. Similar to the `summary()`

function, the `fitted()`

and `residuals()`

functions also take a model object as their input. `fitted()`

calculates all the *fitted* \(\hat{y}\) values by plugging in each observed \(x\) value in the dataset into the regression equation, and `residuals()`

similarly calculates all the residuals (\(y - \hat{y}\)) for each observation in the dataset. The following code will store these new \(\hat{y}\) and residual variables, which we will name `score_hat`

and `residual`

respectively, along with their corresponding `score`

and `bty_avg`

values from the original data `evals_ch5`

.

```
score_model_data <- evals_ch5 %>%
select(score, bty_avg) %>%
mutate(score_hat = fitted(score_model),
residual = residuals(score_model)) %>%
rownames_to_column("ID")
```

Note that we also used the `rownames_to_column()`

function from the `tibble`

package (part of the `tidyverse`

suite) to create a convenient ID column from the rownames. In the table below we only present the results for the 21st through the 24th instructors.

ID | score | bty_avg | score_hat | residual |
---|---|---|---|---|

21 | 4.9 | 7.33 | 4.37 | 0.531 |

22 | 4.6 | 7.33 | 4.37 | 0.231 |

23 | 4.5 | 7.33 | 4.37 | 0.131 |

24 | 4.4 | 5.50 | 4.25 | 0.153 |

Let's recap what we're seeing in Table 5.5 in terms of the linear regression model (`score_model`

) that we fit to the `evals_ch5`

data:

- The
`ID`

column represents the row number where this observation appears in the original`evals_ch5`

data - The
`score`

column represents the observed value of the outcome variable \(y\). - The
`bty_avg`

column represents the values of the explanatory variable \(x\). - The
`score_hat`

column represents the fitted values \(\widehat{y}\). - The
`residual`

column represents the residuals \(y - \widehat{y}\).

Just as we did for the 21st instructor in the `evals_ch5`

dataset (in the first row of the table above), let's repeat the above calculations for the 24th instructor in the `evals_ch5`

dataset (in the fourth row of the table above):

`score`

= 4.4 is the observed value \(y\) for this instructor.`bty_avg`

= 5.50 is the value of the explanatory variable \(x\) for this instructor.`score_hat`

= 4.25 = 3.88 + 0.067 * \(x\) = 3.88 + 0.067 * 5.50 is the fitted value \(\widehat{y}\) for this instructor.`residual`

= 0.153 = 4.4 - 4.25 is the value of the residual for this instructor. In other words, the model was off by 0.153 teaching score units for this instructor.

More development of this idea appears in Section 5.3.3 and we encourage you to read that section after you investigate residuals.

## 5.2 One categorical explanatory variable

It's an unfortunate truth that life expectancy is not the same across various countries in the world; there are a multitude of factors that are associated with how long people live. International development agencies are very interested in studying these differences in the hope of understanding where governments should allocate resources to address this problem. In this section, we'll explore differences in life expectancy in two ways:

- Differences between continents: Are there significant differences in life expectancy, on average, between the five continents of the world: Africa, the Americas, Asia, Europe, and Oceania?
- Differences within continents: How does life expectancy vary within the world's five continents? For example, is the spread of life expectancy among the countries of Africa larger than the spread of life expectancy among the countries of Asia?

To answer such questions, we'll study the `gapminder`

dataset in the `gapminder`

package. Recall we mentioned this dataset in Subsection 2.1.2 when we first studied the "Grammar of Graphics" introduced in Figure 2.1. This dataset has international development statistics such as life expectancy, GDP per capita, and population by country (\(n\) = 142) for 5-year intervals between 1952 and 2007.

We'll use this data for linear regression again, but note that our explanatory variable \(x\) is now categorical, and not numerical like when we covered simple linear regression in Section 5.1. More precisely, we have:

- A numerical outcome variable \(y\). In this case, life expectancy.
- A single categorical explanatory variable \(x\). In this case, the continent the country is part of.

When the explanatory variable \(x\) is categorical, the concept of a "best-fitting" line is a little different than the one we saw previously in Section 5.1 where the explanatory variable \(x\) was numerical. We'll study these differences shortly in Subsection 5.2.2, but first we conduct our exploratory data analysis.

### 5.2.1 Exploratory data analysis

The data on the 142 countries can be found in the `gapminder`

data frame included in the `gapminder`

package. However, to keep things simple, let's `filter()`

for only those observations/rows corresponding to the year 2007. Additionally, let's `select()`

only the subset of the variables we'll consider in this chapter. We'll save this data in a new data frame called `gapminder2007`

:

```
library(gapminder)
gapminder2007 <- gapminder %>%
filter(year == 2007) %>%
select(country, lifeExp, continent, gdpPercap)
```

Recall from Section 5.1.1 that there are three common steps in an exploratory data analysis:

- Most crucially: Looking at the raw data values.
- Computing summary statistics, like means, medians, and interquartile ranges.
- Creating data visualizations.

Let's perform the first common step in an exploratory data analysis: looking at the raw data values. You can do this by using RStudio's spreadsheet viewer or by using the `glimpse()`

command as introduced in Section 1.4.3 on exploring data frames:

`glimpse(gapminder2007)`

```
Rows: 142
Columns: 4
$ country <fct> "Afghanistan", "Albania", "Algeria", "Angola", "Argentina", …
$ lifeExp <dbl> 43.8, 76.4, 72.3, 42.7, 75.3, 81.2, 79.8, 75.6, 64.1, 79.4, …
$ continent <fct> Asia, Europe, Africa, Africa, Americas, Oceania, Europe, Asi…
$ gdpPercap <dbl> 975, 5937, 6223, 4797, 12779, 34435, 36126, 29796, 1391, 336…
```

Observe that `Observations: 142`

indicates that there are 142 rows/observations in `gapminder2007`

, where each row corresponds to one country. In other words, the *observational unit* is an individual country. Furthermore, observe that the variable `continent`

is of type `<fct>`

, which stands for *factor*, which is R's way of encoding categorical variables.

A full description of all the variables included in `gapminder`

can be found by reading the associated help file (run `?gapminder`

in the console). However, let's fully describe the 4 variables we selected in `gapminder2007`

:

`country`

: An identification text variable used to distinguish the 142 countries in the dataset.`lifeExp`

: A numerical variable of that country's life expectancy at birth. This is the outcome variable \(y\) of interest.`continent`

: A categorical variable with five levels. Here "levels" corresponds to the possible categories: Africa, Asia, Americas, Europe, and Oceania. This is the explanatory variable \(x\) of interest.`gdpPercap`

: A numerical variable of that country's GDP per capita in US inflation-adjusted dollars that we'll use as another outcome variable \(y\) in the Learning Check at the end of this subsection.

Furthermore, let's look at a random sample of five out of the 142 countries in Table 5.6. Again, note due to the random nature of the sampling, you will likely end up with a different subset of 5 rows.

```
gapminder2007 %>%
sample_n(size = 5)
```

country | lifeExp | continent | gdpPercap |
---|---|---|---|

Togo | 58.4 | Africa | 883 |

Sao Tome and Principe | 65.5 | Africa | 1598 |

Congo, Dem. Rep. | 46.5 | Africa | 278 |

Lesotho | 42.6 | Africa | 1569 |

Bulgaria | 73.0 | Europe | 10681 |

Now that we've looked at the raw values in our `gapminder2007`

data frame and got a sense of the data, let's move on to computing summary statistics. Let's once again apply the `skim()`

function from the `skimr`

package. Recall from our previous EDA that this function takes in a data frame, "skims" it, and returns commonly used summary statistics. Let's take our `gapminder2007`

data frame, `select()`

only the outcome and explanatory variables `lifeExp`

and `continent`

, and pipe them into the `skim()`

function:

```
gapminder2007 %>%
select(lifeExp, continent) %>%
skim()
```

```
Skim summary statistics
n obs: 142
n variables: 2
── Variable type:factor ────────────────────────────────────────────────────────
variable missing complete n n_unique top_counts ordered
continent 0 142 142 5 Afr: 52, Asi: 33, Eur: 30, Ame: 25 FALSE
── Variable type:numeric ───────────────────────────────────────────────────────
variable missing complete n mean sd p0 p25 p50 p75 p100
lifeExp 0 142 142 67.01 12.07 39.61 57.16 71.94 76.41 82.6
```

The `skim()`

output now reports summaries for categorical variables (`Variable type:factor`

) separately from the numerical variables (`Variable type:numeric`

). For the categorical variable `continent`

, it reports:

`missing`

,`complete`

, and`n`

, which are the number of missing, complete, and total number of values as before, respectively.`n_unique`

: The number of unique levels to this variable, corresponding to Africa, Asia, Americas, Europe, and Oceania. This refers to the how many countries are in the data for each continent.`top_counts`

: In this case the top four counts:`Africa`

has 52 countries,`Asia`

has 33,`Europe`

has 30, and`Americas`

has 25. Not displayed is`Oceania`

with 2 countries.`ordered`

: This tells us whether the categorical variable is "ordinal": whether there is an encoded hierarchy (like low, medium, high). In this case,`continent`

is not ordered.

Turning our attention to the summary statistics of the numerical variable `lifeExp`

, we observe that the global median life expectancy in 2007 was 71.94. Thus, half of the world's countries (71 countries) had a life expectancy less than 71.94. The mean life expectancy of 67.01 is lower however. Why is the mean life expectancy lower than the median?

We can answer this question by performing the last of the three common steps in an exploratory data analysis: creating data visualizations. Let's visualize the distribution of our outcome variable \(y\) = `lifeExp`

in Figure 5.8.

```
ggplot(gapminder2007, aes(x = lifeExp)) +
geom_histogram(binwidth = 5, color = "white") +
labs(x = "Life expectancy", y = "Number of countries",
title = "Histogram of distribution of worldwide life expectancies")
```

We see that this data is *left-skewed*, also known as *negatively* skewed: there are a few countries with low life expectancy that are bringing down the mean life expectancy. However, the median is less sensitive to the effects of such outliers; hence, the median is greater than the mean in this case.

Remember, however, that we want to compare life expectancies both between continents and within continents. In other words, our visualizations need to incorporate some notion of the variable `continent`

. We can do this easily with a faceted histogram. Recall from Section 2.6 that facets allow us to split a visualization by the different values of another variable. We display the resulting visualization in Figure 5.9 by adding a `facet_wrap(~ continent, nrow = 2)`

layer.

```
ggplot(gapminder2007, aes(x = lifeExp)) +
geom_histogram(binwidth = 5, color = "white") +
labs(x = "Life expectancy", y = "Number of countries",
title = "Histogram of distribution of worldwide life expectancies") +
facet_wrap(~ continent, nrow = 2)
```

Observe that unfortunately the distribution of African life expectancies is much lower than the other continents, while in Europe life expectancies tend to be higher and furthermore do not vary as much. On the other hand, both Asia and Africa have the most variation in life expectancies. There is the least variation in Oceania, but keep in mind that there are only two countries in Oceania: Australia and New Zealand.

Recall that an alternative method to visualize the distribution of a numerical variable split by a categorical variable is by using a side-by-side boxplot. We map the categorical variable `continent`

to the \(x\)-axis and the different life expectancies within each continent on the \(y\)-axis in Figure 5.10.

```
ggplot(gapminder2007, aes(x = continent, y = lifeExp)) +
geom_boxplot() +
labs(x = "Continent", y = "Life expectancy (years)",
title = "Life expectancy by continent")
```

Some people prefer comparing the distributions of a numerical variable between different levels of a categorical variable using a boxplot instead of a faceted histogram. This is because we can make quick comparisons between the categorical variable's levels with imaginary horizontal lines. For example, observe in Figure 5.10 that we can quickly convince ourselves that Oceania has the highest median life expectancies by drawing an imaginary horizontal line at \(y\) = 80. Furthermore, as we observed in the faceted histogram in Figure 5.9, Africa and Asia have the largest variation in life expectancy as evidenced by their large interquartile ranges (the heights of the boxes).

It’s important to remember however that the solid lines in the middle of the boxes correspond to the medians (the middle value) rather than the mean (the average). So for example, if you look at Asia, the solid line denotes the median life expectancy of around 72 years. This tells us that half of all countries in Asia have a life expectancy below 72 years whereas half have a life expectancy above 72 years.

Let's compute the median and mean life expectancy for each continent with a little more data wrangling and display the results in Table 5.7.

```
lifeExp_by_continent <- gapminder2007 %>%
group_by(continent) %>%
summarize(median = median(lifeExp), mean = mean(lifeExp))
```

continent | median | mean |
---|---|---|

Africa | 52.9 | 54.8 |

Americas | 72.9 | 73.6 |

Asia | 72.4 | 70.7 |

Europe | 78.6 | 77.6 |

Oceania | 80.7 | 80.7 |

Observe the order of the second column `median`

life expectancy: Africa is lowest, the Americas and Asia are next with similar medians, then Europe, then Oceania. This ordering corresponds to the ordering of the solid black lines inside the boxes in our side-by-side boxplot in Figure 5.10.

Let's now turn our attention to the values in the third column `mean`

. Using Africa's mean life expectancy of 54.8 as a *baseline for comparison*, let's start making relative comparisons to the life expectancies of the other four continents:

- The mean life expectancy of the Americas is 73.6 - 54.8 = 18.8 years higher.
- The mean life expectancy of Asia is 70.7 - 54.8 = 15.9 years higher.
- The mean life expectancy of Europe is 77.6 - 54.8 = 22.8 years higher.
- The mean life expectancy of Oceania is 80.7 - 54.8 = 25.9 years higher.

Let's put these values in Table 5.8, which we'll revisit later on in this section.

continent | mean | Difference versus Africa |
---|---|---|

Africa | 54.8 | 0.0 |

Americas | 73.6 | 18.8 |

Asia | 70.7 | 15.9 |

Europe | 77.6 | 22.8 |

Oceania | 80.7 | 25.9 |

*Learning check*

**(LC5.3)** Conduct a new exploratory data analysis with the same explanatory variable \(x\) being `continent`

but with `gdpPercap`

as the new outcome variable \(y\). Remember, this involves three things:

- Most crucially: Looking at the raw data values.
- Computing summary statistics, such as means, medians, and interquartile ranges.
- Creating data visualizations.

What can you say about the differences in GDP per capita between continents based on this exploration?

### 5.2.2 Linear regression

In Subsection 5.1.2 we introduced *simple* linear regression, which involves modeling the relationship between a numerical outcome variable \(y\) as a function of a numerical explanatory variable \(x\), in our life expectancy example, we now have a categorical explanatory variable \(x\) `continent`

. While we still can fit a regression model, given our categorical explanatory variable we no longer have a concept of a "best-fitting" line, but rather "differences relative to a baseline for comparison."

Before we fit our regression model, let's create a table similar to Table 5.7, but

- Report the mean life expectancy for each continent.
- Report the difference in mean life expectancy
*relative*to Africa's mean life expectancy of 54.806 in the column "mean vs Africa"; this column is simply the "mean" column minus 54.806.

Think back to your observations from the eyeball test of Figure 5.10 at the end of the last subsection. The column "mean vs Africa" is the same idea of comparing a summary statistic to a baseline for comparison, in this case the countries of Africa, but using means instead of medians.

continent | mean | mean vs Africa |
---|---|---|

Africa | 54.8 | 0.0 |

Americas | 73.6 | 18.8 |

Asia | 70.7 | 15.9 |

Europe | 77.6 | 22.8 |

Oceania | 80.7 | 25.9 |

Now, let's use the `summary()`

function again to get the *regression table* for the `gapminder2007`

analysis:

```
lifeExp_model <- lm(lifeExp ~ continent, data = gapminder2007)
summary(lifeExp_model)$coefficients
```

Estimate | Std. Error | t value | Pr(>|t|) | |
---|---|---|---|---|

(Intercept) | 54.8 | 1.02 | 53.45 | 0 |

continentAmericas | 18.8 | 1.80 | 10.45 | 0 |

continentAsia | 15.9 | 1.65 | 9.68 | 0 |

continentEurope | 22.8 | 1.70 | 13.47 | 0 |

continentOceania | 25.9 | 5.33 | 4.86 | 0 |

Just as before, we're only interested in the "Estimate" column, but unlike before, we now have 5 rows corresponding to 5 outputs in our table: an intercept like before, but also `continentAmericas`

, `continentAsia`

, `continentEurope`

, and `continentOceania`

. What are these values? First, we must describe the equation for fitted value \(\widehat{y}\), which is a little more complicated when the \(x\) explanatory variable is categorical:

\[ \begin{align} \widehat{\text{life exp}} &= b_0 + b_{\text{Amer}}\cdot\mathbb{1}_{\mbox{Amer}}(x) + b_{\text{Asia}}\cdot\mathbb{1}_{\mbox{Asia}}(x) + b_{\text{Euro}}\cdot\mathbb{1}_{\mbox{Euro}}(x) + b_{\text{Ocean}}\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot\mathbb{1}_{\mbox{Amer}}(x) + 15.9\cdot\mathbb{1}_{\mbox{Asia}}(x) + 22.8\cdot\mathbb{1}_{\mbox{Euro}}(x) + 25.9\cdot\mathbb{1}_{\mbox{Ocean}}(x) \end{align} \]

Let's break this down. First, \(\mathbb{1}_{A}(x)\) is what's known in mathematics as an "indicator function" that takes one of two possible values:

\[ \mathbb{1}_{A}(x) = \left\{ \begin{array}{ll} 1 & \text{if } x \text{ is in } A \\ 0 & \text{otherwise} \end{array} \right. \]

In a statistical modeling context this is also known as a "dummy variable". In our case, let's consider the first such indicator variable:

\[ \mathbb{1}_{\mbox{Amer}}(x) = \left\{ \begin{array}{ll} 1 & \text{if } \text{country } x \text{ is in the Americas} \\ 0 & \text{otherwise}\end{array} \right. \]

Now let's interpret the terms in the estimate column of the regression table. First \(b_0 =\) `intercept = 54.8`

corresponds to the mean life expectancy for countries in Africa, since for country \(x\) in Africa we have the following equation:

\[ \begin{align} \widehat{\text{life exp}} &= b_0 + b_{\text{Amer}}\cdot\mathbb{1}_{\mbox{Amer}}(x) + b_{\text{Asia}}\cdot\mathbb{1}_{\mbox{Asia}}(x) + b_{\text{Euro}}\cdot\mathbb{1}_{\mbox{Euro}}(x) + b_{\text{Ocean}}\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot\mathbb{1}_{\mbox{Amer}}(x) + 15.9\cdot\mathbb{1}_{\mbox{Asia}}(x) + 22.8\cdot\mathbb{1}_{\mbox{Euro}}(x) + 25.9\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot 0 + 15.9\cdot 0 + 22.8\cdot 0 + 25.9\cdot 0\\ &= 54.8 \end{align} \]

i.e. All four of the indicator variables are equal to 0. Recall we stated earlier that we would treat Africa as the baseline for comparison group. Furthermore, this value corresponds to the group mean life expectancy for all African countries in Table 5.8.

Next, \(b_{\text{Amer}}\) = `continentAmericas = 18.8`

is the difference in mean life expectancy of countries in the Americas relative to Africa, or in other words, on average countries in the Americas had life expectancy 18.8 years greater. The fitted value yielded by this equation is:

\[ \begin{align} \widehat{\text{life exp}} &= b_0 + b_{\text{Amer}}\cdot\mathbb{1}_{\mbox{Amer}}(x) + b_{\text{Asia}}\cdot\mathbb{1}_{\mbox{Asia}}(x) + b_{\text{Euro}}\cdot\mathbb{1}_{\mbox{Euro}}(x) + b_{\text{Ocean}}\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot\mathbb{1}_{\mbox{Amer}}(x) + 15.9\cdot\mathbb{1}_{\mbox{Asia}}(x) + 22.8\cdot\mathbb{1}_{\mbox{Euro}}(x) + 25.9\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot 1 + 15.9\cdot 0 + 22.8\cdot 0 + 25.9\cdot 0\\ &= 54.8 + 18.8\\ &= 73.6 \end{align} \]

i.e. in this case, only the indicator function \(\mathbb{1}_{\mbox{Amer}}(x)\) is equal to 1, but all others are 0. Recall that 73.6 corresponds to the group mean life expectancy for all countries in the Americas in Table 5.8.

Similarly, \(b_{\text{Asia}}\) = `continentAsia = 15.9`

is the difference in mean life expectancy of Asian countries relative to Africa countries, or in other words, on average countries in the Asia had life expectancy 15.9 years greater than Africa. The fitted value yielded by this equation is:

\[ \begin{align} \widehat{\text{life exp}} &= b_0 + b_{\text{Amer}}\cdot\mathbb{1}_{\mbox{Amer}}(x) + b_{\text{Asia}}\cdot\mathbb{1}_{\mbox{Asia}}(x) + b_{\text{Euro}}\cdot\mathbb{1}_{\mbox{Euro}}(x) + b_{\text{Ocean}}\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot\mathbb{1}_{\mbox{Amer}}(x) + 15.9\cdot\mathbb{1}_{\mbox{Asia}}(x) + 22.8\cdot\mathbb{1}_{\mbox{Euro}}(x) + 25.9\cdot\mathbb{1}_{\mbox{Ocean}}(x)\\ &= 54.8 + 18.8\cdot 0 + 15.9\cdot 1 + 22.8\cdot 0 + 25.9\cdot 0\\ &= 54.8 + 15.9\\ &= 70.7 \end{align} \]

i.e. in this case, only the indicator function \(\mathbb{1}_{\mbox{Asia}}(x)\) is equal to 1, but all others are 0. Recall that 70.7 corresponds to the group mean life expectancy for all countries in Asia in Table 5.8. The same logic applies to \(b_{\text{Euro}} = 22.8\) and \(b_{\text{Ocean}} = 25.9\); they correspond to the "offset" in mean life expectancy for countries in Europe and Oceania, relative to the mean life expectancy of the baseline group for comparison of African countries.

Let's generalize this idea a bit. If we fit a linear regression model using a categorical explanatory variable \(x\) that has \(k\) levels, a regression model will return an intercept and \(k - 1\) "slope" coefficients. When \(x\) is a numerical explanatory variable the interpretation is of a "slope" coefficient, but when \(x\) is categorical the meaning is a little trickier. They are *offsets* relative to the baseline.

In our case, since there are \(k = 5\) continents, the regression model returns an intercept corresponding to the baseline for comparison Africa and \(k - 1 = 4\) slope coefficients corresponding to the Americas, Asia, Europe, and Oceania. Africa was chosen as the baseline by R for no other reason than it is first alphabetically of the 5 continents. You can manually specify which continent to use as baseline instead of the default choice of whichever comes first alphabetically, but we leave that to a more advanced course. (The `forcats`

package is particularly nice for doing this and we encourage you to explore using it.)

*Learning check*

**(LC5.4)** Fit a new linear regression using `lm(gdpPercap ~ continent, data = gapminder2007)`

where `gdpPercap`

is the new outcome variable \(y\). Get information about the "best-fitting" line from the regression table by applying the `summary()`

function. How do the regression results match up with the results from your exploratory data analysis above?

### 5.2.3 Observed/fitted values and residuals

Recall in Subsection 5.1.3 we defined the following three concepts:

- Observed values \(y\), or the observed value of the outcome variable
- Fitted values \(\widehat{y}\), or the value on the regression line for a given \(x\) value
- Residuals \(y - \widehat{y}\), or the error between the observed value and the fitted value

Let's use similar code (from Subsection 5.1.3) to obtain these values for the `lifeExp_model`

fit to the `gapminder2007`

dataset:

```
lifeExp_model_data <- gapminder2007 %>%
select(country, lifeExp, continent) %>%
mutate(lifeExp_hat = fitted(lifeExp_model),
residual = residuals(lifeExp_model))
```

country | lifeExp | continent | lifeExp_hat | residual |
---|---|---|---|---|

Afghanistan | 43.8 | Asia | 70.7 | -26.900 |

Albania | 76.4 | Europe | 77.6 | -1.226 |

Algeria | 72.3 | Africa | 54.8 | 17.495 |

Angola | 42.7 | Africa | 54.8 | -12.075 |

Argentina | 75.3 | Americas | 73.6 | 1.712 |

Australia | 81.2 | Oceania | 80.7 | 0.516 |

Austria | 79.8 | Europe | 77.6 | 2.180 |

Bahrain | 75.6 | Asia | 70.7 | 4.907 |

Bangladesh | 64.1 | Asia | 70.7 | -6.666 |

Belgium | 79.4 | Europe | 77.6 | 1.792 |

Observe in Table 5.11 that `lifeExp_hat`

contains the fitted values \(\hat{y} = \widehat{lifeExp}\). If you look closely, there are only 5 possible values for `lifeExp_hat`

. These correspond to the five mean life expectancies for the 5 continents that we displayed in Table 5.8 and computed using the "Estimate" column of the regression results in Table 5.10.

The `residual`

column is simply \(y - \widehat{y}\) = `lifeexp - lifeexp_hat`

. These values can be interpreted as that particular country's deviation from the mean life expectancy of the respective continent's mean. For example, the first row of Table 5.11 corresponds to Afghanistan, and the residual of \(y - \widehat{y} = 43.8 - 70.7 = -26.9\) is telling us that Afghanistan's life expectancy is a whopping 26.9 years lower than the mean life expectancy of all Asian countries. This can in part be explained by the many years of war that country has suffered.

*Learning check*

**(LC5.5)** Using either the sorting functionality of RStudio's spreadsheet viewer or using the data wrangling tools you learned in Chapter 3, identify the five countries with the five smallest (most negative) residuals? What do these negative residuals say about their life expectancy relative to their continents?

**(LC5.6)** Repeat this process, but identify the five countries with the five largest (most positive) residuals. What do these positive residuals say about their life expectancy relative to their continents?

## 5.4 Conclusion

### 5.4.1 Additional resources

As we suggested in Subsection 5.1.1, interpreting coefficients that are not close to the extreme values of -1, 0, and 1 can be somewhat subjective. To help develop your sense of correlation coefficients, we suggest you play the following 80's-style video game called "Guess the correlation" at http://guessthecorrelation.com/.

### 5.4.2 What's to come?

In this chapter, you've studied the term *basic regression*, where you fit models that only have one explanatory variable. In Chapter 6, we'll study *multiple regression*, where our regression models can now have more than one explanatory variable! In particular, we'll consider two scenarios: regression models with one numerical and one categorical explanatory variable and regression models with two numerical explanatory variables. This will allow you to construct more sophisticated and more powerful models, all in the hopes of better explaining your outcome variable \(y\).

### References

Waring, Elin, Michael Quinn, Amelia McNamara, Eduardo Arino de la Rubia, Hao Zhu, and Shannon Ellis. 2021. *Skimr: Compact and Flexible Summaries of Data*. https://CRAN.R-project.org/package=skimr.