# 7 Multiple Regression

In Chapter 6 we introduced ideas related to modeling, in particular that the fundamental premise of modeling is *to make explicit the relationship* between an outcome variable \(y\) and an explanatory/predictor variable \(x\). Recall further the synonyms that we used to also denote \(y\) as the dependent variable and \(x\) as an independent variable or covariate.

There are many modeling approaches one could take, among the most well-known being linear regression, which was the focus of the last chapter. Whereas in the last chapter we only focused on regression scenarios where there is only one explanatory/predictor variable, in this chapter, we now focus on modeling scenarios where there is more than one; this is known as multiple regression. You can imagine when trying to model a particular outcome variable, like teaching evaluation score as in Section 6.1 or life expectancy as in Section 6.2, it would be very useful to incorporate more than one explanatory variable.

Since our regression models will now consider more than one explanatory/predictor variable, the interpretation of the associated effect of any one explanatory/predictor variables must be made in conjunction with the others. For example, say we are modeling individuals’ incomes as a function of their number of years of education and their parents’ wealth. When interpreting the effect of education on income, one has to consider the effect of their parents’ wealth at the same time, as these two variables are almost certainly related. Make note of this throughout this chapter and as you work on interpreting the results of multiple regression models into the future.

### Needed packages

Let’s load all the packages needed for this chapter (this assumes you’ve already installed them). Read Section 2.3 for information on how to install and load R packages.

```
library(ggplot2)
library(dplyr)
library(moderndive)
library(ISLR)
```

## 7.1 Two numerical explanatory variables

Let’s now attempt to identify factors that are associated with how much credit card debt an individual will have. The textbook An Introduction to Statistical Learning with Applications in R by Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani is an intermediate-level textbook on statistical and machine learning freely available here. It has an accompanying R package called `ISLR`

with datasets that the authors use to demonstrate various machine learning methods. One dataset that is frequently used by the authors is the `Credit`

dataset where predictions are made on the credit card balance held by \(n = 400\) credit card holders based on information about them like income, credit limit, and education level.

Since no information was provided as to who these \(n\) = 400 individuals are and how they came to be included in this dataset, it will be hard to make any scientific claims based on this data. Recall our discussion from the previous chapter that correlation does not necessarily imply causation. That being said, we’ll still use `Credit`

to demonstrate multiple regression with:

- A numerical outcome variable \(y\), in this case credit card balance.
- Two explanatory variables:
- A first numerical explanatory variable \(x_1\). In this case, their credit limit.
- A second numerical explanatory variable \(x_2\). In this case, their income (in thousands of dollars).

In the forthcoming Learning Checks, we’ll consider a different scenario:

- The same numerical outcome variable \(y\): credit card balance.
- Two new explanatory variables:
- A first numerical explanatory variable \(x_1\): their credit rating.
- A second numerical explanatory variable \(x_2\): their age.

### 7.1.1 Exploratory data analysis

Let’s load the `Credit`

data and `select()`

only the needed subset of variables.

```
library(ISLR)
Credit <- Credit %>%
select(Balance, Limit, Income, Rating, Age)
```

Let’s look at the raw data values both by bringing up RStudio’s spreadsheet viewer and the `glimpse()`

function. Although in Table 7.1 we only show 5 randomly selected credit card holders out of 400:

`View(Credit)`

Balance | Limit | Income | Rating | Age | |
---|---|---|---|---|---|

141 | 1425 | 6045 | 39.8 | 459 | 32 |

9 | 279 | 3300 | 15.1 | 266 | 66 |

13 | 204 | 5308 | 80.6 | 394 | 57 |

262 | 1050 | 9310 | 180.4 | 665 | 67 |

267 | 15 | 4952 | 88.8 | 360 | 86 |

`glimpse(Credit)`

```
Observations: 400
Variables: 5
$ Balance <int> 333, 903, 580, 964, 331, 1151, 203, 872, 279, 1350, 1407, 0...
$ Limit <int> 3606, 6645, 7075, 9504, 4897, 8047, 3388, 7114, 3300, 6819,...
$ Income <dbl> 14.9, 106.0, 104.6, 148.9, 55.9, 80.2, 21.0, 71.4, 15.1, 71...
$ Rating <int> 283, 483, 514, 681, 357, 569, 259, 512, 266, 491, 589, 138,...
$ Age <int> 34, 82, 71, 36, 68, 77, 37, 87, 66, 41, 30, 64, 57, 49, 75,...
```

Let’s look at some summary statistics:

```
Credit %>%
select(Balance, Limit, Income) %>%
summary()
```

```
Balance Limit Income
Min. : 0 Min. : 855 Min. : 10.4
1st Qu.: 69 1st Qu.: 3088 1st Qu.: 21.0
Median : 460 Median : 4622 Median : 33.1
Mean : 520 Mean : 4736 Mean : 45.2
3rd Qu.: 863 3rd Qu.: 5873 3rd Qu.: 57.5
Max. :1999 Max. :13913 Max. :186.6
```

We observe, for example,

- The mean and median credit card balance is around $500. 25% of card holders had debts of 69 dollars or less.
- The mean and median credit card limit is just under $5000.
- 75% of these card holders had incomes of $57,500 or less.

Since our outcome variable `Balance`

and the explanatory variables `Limit`

and `Rating`

are numerical, we can compute the correlation coefficient between pairs of these variables. There are two ways of doing this. First, we could run the `cor()`

command as seen in Subsection 6.1.1 twice, once for each explanatory variable:

```
cor(Credit$Balance, Credit$Limit)
cor(Credit$Balance, Credit$Income)
```

Or we can simultaneously compute them by returning a *correlation matrix* in Table 7.2. We can read off the correlation coefficient for any pair of variables by looking them up in the appropriate row/column combination.

```
Credit %>%
select(Balance, Limit, Income) %>%
cor()
```

Balance | Limit | Income | |
---|---|---|---|

Balance | 1.000 | 0.862 | 0.464 |

Limit | 0.862 | 1.000 | 0.792 |

Income | 0.464 | 0.792 | 1.000 |

For example, the correlation coefficient of:

`Balance`

with itself is 1 as we would expect based on the definition of the correlation coefficient.`Balance`

with`Limit`

is 0.862. This indicates a strong positive linear relationship, which makes sense as only individuals with large credit limits can accrue large credit card balances.`Balance`

with`Income`

is 0.464. This is suggestive of another positive linear relationship, although not as strong as the relationship between`Balance`

and`Limit`

.- As an added bonus, we can read off the correlation coefficient of the two explanatory variables,
`Limit`

and`Income`

of 0.792. In this case, we say there is a high degree of*collinearity*between these two explanatory variables.

Collinearity (or multicollinearity) is a phenomenon in which one explanatory variable in a multiple regression model can be linearly predicted from the others with a substantial degree of accuracy. So in this case, if we knew someone’s credit card `Limit`

and since `Limit`

and `Income`

are highly correlated, we could make a fairly accurate guess as to that person’s `Income`

. Or put loosely, these two variables provided redundant information. For now let’s ignore any issues related to collinearity and press on.

Let’s visualize the relationship of the outcome variable with each of the two explanatory variables in two separate plots:

```
ggplot(Credit, aes(x = Limit, y = Balance)) +
geom_point() +
labs(x = "Credit limit (in $)", y = "Credit card balance (in $)",
title = "Relationship between balance and credit limit") +
geom_smooth(method = "lm", se = FALSE)
ggplot(Credit, aes(x = Income, y = Balance)) +
geom_point() +
labs(x = "Income (in $1000)", y = "Credit card balance (in $)",
title = "Relationship between balance and income") +
geom_smooth(method = "lm", se = FALSE)
```

First, there is a positive relationship between credit limit and balance, since as credit limit increases so also does credit card balance; this is to be expected given the strongly positive correlation coefficient of 0.862. In the case of income, the positive relationship doesn’t appear as strong, given the weakly positive correlation coefficient of 0.464. However the two plots in Figure 7.1 only focus on the relationship of the outcome variable with each of the explanatory variables independently. To get a sense of the *joint* relationship of all three variables simultaneously through a visualization, let’s display the data in a 3-dimensional (3D) scatterplot, where

- The numerical outcome variable \(y\)
`Balance`

is on the z-axis (vertical axis) - The two numerical explanatory variables form the “floor” axes. In this case
- The first numerical explanatory variable \(x_1\)
`Income`

is on of the floor axes. - The second numerical explanatory variable \(x_2\)
`Limit`

is on the other floor axis.

- The first numerical explanatory variable \(x_1\)

Click on the following image to open an interactive 3D scatterplot in your browser:

Previously in Figure 6.5, we plotted a “best-fitting” regression line through a set of points where the numerical outcome variable \(y\) was teaching `score`

and a single numerical explanatory variable \(x\) `bty_avg`

. What is the analogous concept when we have *two* numerical predictor variables? Instead of a best-fitting line, we now have a best-fitting *plane*, which is 3D generalization of lines which exist in 2D. Click on the following image to open an interactive plot of the regression plane in your browser. Move the image around, zoom in, and think about how this plane generalizes the concept of a linear regression line to three dimensions.

*Learning check*

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

but with `Rating`

and `Age`

as the new explanatory variables \(x_1\) and \(x_2\). 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 a credit card holder’s balance and their credit rating and age?

### 7.1.2 Multiple regression

Just as we did when we had a single numerical explanatory variable \(x\) in Subsection 6.1.2 and when we had a single categorical explanatory variable \(x\) in Subsection 6.2.2, we fit a regression model and get the regression table in our two numerical explanatory variable scenario. To fit a regression model and get a table using `get_regression_table()`

, we now use a `+`

to consider multiple explanatory variables. In this case since we want to preform a regression of `Limit`

and `Income`

simultaneously, we input `Balance ~ Limit + Income`

.

```
Balance_model <- lm(Balance ~ Limit + Income, data = Credit)
get_regression_table(Balance_model)
```

term | estimate | std_error | statistic | p_value | lower_ci | upper_ci |
---|---|---|---|---|---|---|

intercept | -385.179 | 19.465 | -19.8 | 0 | -423.446 | -346.912 |

Limit | 0.264 | 0.006 | 45.0 | 0 | 0.253 | 0.276 |

Income | -7.663 | 0.385 | -19.9 | 0 | -8.420 | -6.906 |

How do we interpret these three values that define the regression plane?

- Intercept: -$385.18 (rounded to two decimal points to represent cents). The intercept in our case represents the credit card balance for an individual who has both a credit
`Limit`

of $0 and`Income`

of $0. In our data however, the intercept has limited practical interpretation as no individuals had`Limit`

or`Income`

values of $0 and furthermore the smallest credit card balance was $0. Rather, it is used to situate the regression plane in 3D space. - Limit: $0.26. Now that we have multiple variables to consider, we have to add a caveat to our interpretation:
*all other things being equal, for every increase of one unit in credit*. Note:`Limit`

(dollars), there is an associated increase of on average $0.26 in credit card balance- Just as we did in Subsection 6.1.2, we are not making any causal statements, only statements relating to the association between credit limit and balance
- The
*all other things being equal*is making a statement about all other explanatory variables, in this case only one:`Income`

. This is equivalent to saying “holding`Income`

constant, we observed an associated increase of $0.26 in credit card balance for every dollar increase in credit limit.”

- Income: -$7.66. Similarly,
*all other things being equal, for every increase of one unit in*.`Income`

(in other words, $1000 in income), there is an associated decrease of on average $7.66 in credit card balance

However, recall in Figure 7.1 that when considered separately, both `Limit`

and `Income`

had positive relationships with the outcome variable `Balance`

: as card holders’ credit limits increased their credit card balances tended to increase as well, and a similar relationship held for incomes and balances. In the above multiple regression, however, the slope for `Income`

is now -7.66, suggesting a *negative relationship* between income and credit card balance. What explains these contradictory results?

This is known as Simpson’s Paradox, a phenomenon in which a trend appears in several different groups of data but disappears or reverses when these groups are combined. We expand on this in Subsection 7.3.2 where we’ll look at the relationship between credit `Limit`

and credit card balance but split by different income bracket groups.

*Learning check*

**(LC7.2)** Fit a new simple linear regression using `lm(Balance ~ Rating + Age, data = Credit)`

where `Rating`

and `Age`

are the new numerical explanatory variables \(x_1\) and \(x_2\). Get information about the “best-fitting” line from the regression table by applying the `get_regression_table()`

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

### 7.1.3 Observed/fitted values and residuals

As we did previously, in Table 7.4 let’s unpack the output of the `get_regression_points()`

function for our model for credit card balance for all 400 card holders in the dataset. Recall that each card holder corresponds to one of the 400 rows in the `Credit`

data frame and also for one of the 400 points in the 3D scatterplots in Subsection 7.1.1.

```
regression_points <- get_regression_points(Balance_model)
regression_points
```

ID | Balance | Limit | Income | Balance_hat | residual |
---|---|---|---|---|---|

1 | 333 | 3606 | 14.9 | 454 | -120.8 |

2 | 903 | 6645 | 106.0 | 559 | 344.3 |

3 | 580 | 7075 | 104.6 | 683 | -103.4 |

4 | 964 | 9504 | 148.9 | 986 | -21.7 |

5 | 331 | 4897 | 55.9 | 481 | -150.0 |

Recall the format of the output:

`Balance`

corresponds to \(y\) (the observed value)`Balance_hat`

corresponds to \(\widehat{y}\) (the fitted value)`residual`

corresponds to \(y - \widehat{y}\) (the residual)

### 7.1.4 Residual analysis

Recall in Section 6.1.4, our first residual analysis plot investigated the presence of any systematic pattern in the residuals when we had a single numerical predictor: `bty_age`

. For the `Credit`

card dataset, since we have two numerical predictors, `Limit`

and `Income`

, we must perform this twice:

```
ggplot(regression_points, aes(x = Limit, y = residual)) +
geom_point() +
labs(x = "Credit limit (in $)", y = "Residual", title = "Residuals vs credit limit")
ggplot(regression_points, aes(x = Income, y = residual)) +
geom_point() +
labs(x = "Income (in $1000)", y = "Residual", title = "Residuals vs income")
```

In this case, there **does** appear to be a systematic pattern to the residuals. As the scatter of the residuals around the line \(y=0\) is definitely not consistent. This behavior of the residuals is further evidenced by the histogram of residuals in Figure 7.3. We observe that the residuals have a slight right-skew (recall we say that data is right-skewed, or positively-skewed, if there is a tail to the right). Ideally, these residuals should be bell-shaped around a residual value of 0.

```
ggplot(regression_points, aes(x = residual)) +
geom_histogram(color = "white") +
labs(x = "Residual")
```

Another way to interpret this histogram is that since the residual is computed as \(y - \widehat{y}\) = `balance`

- `balance_hat`

, we have some values where the fitted value \(\widehat{y}\) is very much lower than the observed value \(y\). In other words, we are underestimating certain credit card holders’ balances by a very large amount.

*Learning check*

**(LC7.3)** Continuing with our regression using `Rating`

and `Age`

as the explanatory variables and credit card `Balance`

as the outcome variable, use the `get_regression_points()`

function to get the observed values, fitted values, and residuals for all 400 credit card holders. Perform a residual analysis and look for any systematic patterns in the residuals.

## 7.2 One numerical & one categorical explanatory variable

Let’s revisit the instructor evaluation data introduced in Section 6.1, where we studied the relationship between instructor evaluation scores and their beauty scores. This analysis suggested that there is a positive relationship between `bty_avg`

and `score`

, in other words as instructors had higher beauty scores, they also tended to have higher teaching evaluation scores. Now let’s say instead of `bty_avg`

we are interested in the numerical explanatory variable \(x_1\) `age`

and furthermore we want to use a second explanatory variable \(x_2\), the (binary) categorical variable `gender`

. Our modeling scenario now becomes

- A numerical outcome variable \(y\). As before, instructor evaluation score.
- Two explanatory variables:
- A numerical explanatory variable \(x_1\): in this case, their age.
- A categorical explanatory variable \(x_2\): in this case, their binary gender.

### 7.2.1 Exploratory data analysis

Let’s reload the `evals`

data and `select()`

only the needed subset of variables.

```
load(url("http://www.openintro.org/stat/data/evals.RData"))
evals <- evals %>%
select(score, age, gender)
```

Let’s look at the raw data values both by bringing up RStudio’s spreadsheet viewer and the `glimpse()`

function, although in Table 7.5 we only show 5 randomly selected instructors out of 463:

`View(evals)`

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

290 | 3.6 | 6.67 | 34 | male |

341 | 4.9 | 3.50 | 43 | male |

199 | 3.3 | 2.33 | 47 | male |

47 | 4.4 | 4.67 | 33 | female |

215 | 4.7 | 3.67 | 60 | male |

Let’s look at some summary statistics:

`summary(evals)`

```
score bty_avg age gender
Min. :2.30 Min. :1.67 Min. :29.0 female:195
1st Qu.:3.80 1st Qu.:3.17 1st Qu.:42.0 male :268
Median :4.30 Median :4.33 Median :48.0
Mean :4.17 Mean :4.42 Mean :48.4
3rd Qu.:4.60 3rd Qu.:5.50 3rd Qu.:57.0
Max. :5.00 Max. :8.17 Max. :73.0
```

In Figure 7.4, we plot a scatterplot of `score`

over `age`

, but given that `gender`

is a binary categorical variable

- We can assign a color to points from each of the two levels of
`gender`

: female and male. - Furthermore, the
`geom_smooth(method = "lm", se = FALSE)`

layer automatically fits a different regression line for each.

```
ggplot(evals, aes(x = age, y = score, col = gender)) +
geom_jitter() +
labs(x = "Age", y = "Teaching Score", color = "Gender") +
geom_smooth(method = "lm", se = FALSE)
```

We notice some interesting trends:

- There are almost no women faculty over the age of 60.
- Fitting separate regression lines for men and women, we see they have different slopes. We see that the associated effect of increasing age seems to be much harsher for women than men. In other words, as women age, the drop in their teaching score appears to be more faster.

### 7.2.2 Multiple regression

Much like we started to consider multiple explanatory variables using the `+`

sign in Subsection 7.1.2, let’s fit a regression model and get the regression table, this time saving our regression model fit in `score_model_2`

so as to not overwrite the model `score_model`

from Section 6.1.2.

```
score_model_2 <- lm(score ~ age + gender, data = evals)
get_regression_table(score_model_2)
```

term | estimate | std_error | statistic | p_value | lower_ci | upper_ci |
---|---|---|---|---|---|---|

intercept | 4.484 | 0.125 | 35.79 | 0.000 | 4.238 | 4.730 |

age | -0.009 | 0.003 | -3.28 | 0.001 | -0.014 | -0.003 |

gendermale | 0.191 | 0.052 | 3.63 | 0.000 | 0.087 | 0.294 |

The modeling equation for this scenario is:

\[
\begin{align}
\widehat{y} &= b_0 + b_1 * x_1 + b_2 * x_2 \\
\widehat{\mbox{score}} &= b_0 + b_{\mbox{age}} * \mbox{age} + b_{\mbox{male}} * \mathbb{1}[\mbox{is male}] \\
\end{align}
\] where \(\mathbb{1}[\mbox{is male}]\) is an *indicator function* for `sex == male`

. In other words, \(\mathbb{1}[\mbox{is male}]\) equals one if the current observation corresponds to a male professor, and 0 if the current observation corresponds to a female professor. This model can be visualized in Figure 7.5.

We see that:

- Females are treated as the baseline for comparison for no other reason than “female” is alphabetically earlier than “male.” The \(b_{male} = 0.1906\) is the vertical “bump” that men get in their teaching evaluation scores. Or more precisely, it is the average difference in teaching score that men get
*relative to the baseline of women*. - Accordingly, the intercepts are (which in this case make no sense since no instructor can have age 0):
- for women: \(b_0\) = 4.484
- for men: \(b_0 + b_{male}\) = 4.484 + 0.191 = 4.675

- Both men and women have the same slope. In other words,
*in this model*the associated effect of age is the same for men and women. All other things being equal, for every increase in 1 in age, there is on average an associated decrease of \(b_{age}\) = -0.0086 in teaching score.

But wait, why is Figure 7.5 different than Figure 7.4! What is going on? What we have in the original plot is know as an *interaction effect* between age and gender. Focusing on fitting a model for each of men and women we see that the resulting regression lines are different. Thus, `gender`

appears to interact in different ways for men and women with the different values of `age`

.

### 7.2.3 Multiple regression with interaction effects

We say a model has an *interaction effect* if the associated effect of one variable *depends on the value of another variable*. These types of models usually prove to be tricky to view on first glance because of their complexity. In this case, the effect of `age`

will depend on the value of `gender`

. Put differently, the effect of age on teaching scores will differ for men and for women, as was suggested by the different slopes for men and women in our visual exploratory data analysis in Figure 7.4.

Let’s fit a regression with an interaction term. Instead of using the `+`

sign in the enumeration of explanatory variables, we use the `*`

sign. Let’s fit this regression and save it in `score_model_3`

, then we get the regression table using the `get_regression_table()`

function as before.

```
score_model_interaction <- lm(score ~ age * gender, data = evals)
get_regression_table(score_model_interaction)
```

term | estimate | std_error | statistic | p_value | lower_ci | upper_ci |
---|---|---|---|---|---|---|

intercept | 4.883 | 0.205 | 23.80 | 0.000 | 4.480 | 5.286 |

age | -0.018 | 0.004 | -3.92 | 0.000 | -0.026 | -0.009 |

gendermale | -0.446 | 0.265 | -1.68 | 0.094 | -0.968 | 0.076 |

age:gendermale | 0.014 | 0.006 | 2.45 | 0.015 | 0.003 | 0.024 |

The modeling equation for this scenario is:

\[ \begin{align} \widehat{y} &= b_0 + b_1*x_1 + b_2*x_2 + b_3*x_1*x_2\\ \widehat{\mbox{score}} &= b_0 + b_{\mbox{age}}*\mbox{age} + b_{\mbox{male}}*\mathbb{1}[\mbox{is male}] + b_{\mbox{age,male}}*\mbox{age}*\mathbb{1}[\mbox{is male}] \\ \end{align} \]

Oof, that’s a lot of rows in the regression table output and a lot of terms in the model equation. The fourth term being added on the right hand side of the equation corresponds to the *interaction term*. Let’s simplify things by considering men and women separately. First, recall that \(\mathbb{1}[\mbox{is male}]\) equals 1 if a particular observation (or row in `evals`

) corresponds to a male instructor. In this case, using the values from the regression table the fitted value of \(\widehat{\mbox{score}}\) is:

\[ \begin{align} \widehat{\mbox{score}} &= b_0 + b_{\mbox{age}}*\mbox{age} + b_{\mbox{male}}*\mathbb{1}[\mbox{is male}] + b_{\mbox{age,male}}*\mbox{age}*\mathbb{1}[\mbox{is male}] \\ &= b_0 + b_{\mbox{age}}*\mbox{age} + b_{\mbox{male}}*1 + b_{\mbox{age,male}}\mbox{age}*1 \\ &= \left(b_0 + b_{\mbox{male}}\right) + \left(b_{\mbox{age}} + b_{\mbox{age,male}} \right)*\mbox{age} \\ &= \left(4.883 + -0.446\right) + \left(-0.018 + 0.014 \right)*\mbox{age} \\ &= 4.437 -0.004*\mbox{age} \end{align} \]

Second, recall that \(\mathbb{1}[\mbox{is male}]\) equals 0 if a particular observation corresponds to a female instructors. Again, using the values from the regression table the fitted value of \(\widehat{\mbox{score}}\) is:

\[ \begin{align} \widehat{\mbox{score}} &= b_0 + b_{\mbox{age}}*\mbox{age} + b_{\mbox{male}}*\mathbb{1}[\mbox{is male}] + b_{\mbox{age,male}}*\mbox{age}*\mathbb{1}[\mbox{is male}] \\ &= b_0 + b_{\mbox{age}}*\mbox{age} + b_{\mbox{male}}*0 + b_{\mbox{age,male}}\mbox{age}*0 \\ &= b_0 + b_{\mbox{age}}*\mbox{age}\\ &= 4.883 -0.018*\mbox{age} \end{align} \]

Let’s summarize these values in a table:

Gender | Intercept | Slope for age |
---|---|---|

Male instructors | 4.44 | -0.004 |

Female instructors | 4.88 | -0.018 |

We see that while male instructors have a lower intercept, as they age, they have a less steep associated average decrease in teaching scores: 0.004 teaching score units per year as opposed to -0.018 for women. This is consistent with the different slopes and intercepts of the red and blue regression lines fit in Figure 7.4. Recall our definition of a model having an interaction effect: when the associated effect of one variable, in this case `age`

, depends on the value of another variable, in this case `gender`

.

But how do we know when it’s appropriate to include an interaction effect? For example, which is the more appropriate model? The regular multiple regression model without an interaction term we saw in Section 7.2.2 or the multiple regression model with the interaction term we just saw? We’ll revisit this question in Chapter 11 on “inference for regression.”

### 7.2.4 Observed/fitted values and residuals

Now say we want to apply the above calculations for male and female instructors for all 463 instructors in the `evals`

dataset. As our multiple regression models get more and more complex, computing such values by hand gets more and more tedious. The `get_regression_points()`

function spares us this tedium and returns all fitted values and all residuals. For simplicity, let’s focus only on the fitted interaction model, which is saved in `score_model_interaction`

.

```
regression_points <- get_regression_points(score_model_interaction)
regression_points
```

ID | score | age | gender | score_hat | residual |
---|---|---|---|---|---|

1 | 4.7 | 36 | female | 4.25 | 0.448 |

2 | 4.1 | 36 | female | 4.25 | -0.152 |

3 | 3.9 | 36 | female | 4.25 | -0.352 |

4 | 4.8 | 36 | female | 4.25 | 0.548 |

5 | 4.6 | 59 | male | 4.20 | 0.399 |

Recall the format of the output:

`score`

corresponds to \(y\) the observed value`score_hat`

corresponds to \(\widehat{y} = \widehat{\mbox{score}}\) the fitted value`residual`

corresponds to the residual \(y - \widehat{y}\)

### 7.2.5 Residual analysis

As always, let’s perform a residual analysis first with a histogram:

```
ggplot(regression_points, aes(x = residual)) +
geom_histogram(binwidth = 0.25, color = "white") +
labs(x = "Residual")
```

Second, the residuals as compared to the predictor variables:

- \(x_1\): numerical explanatory/predictor variable of
`age`

- \(x_2\): categorical explanatory/predictor variable of
`gender`

```
ggplot(regression_points, aes(x = age, y = residual)) +
geom_point() +
labs(x = "age", y = "Residual") +
geom_hline(yintercept = 0, col = "blue", size = 1) +
facet_wrap(~gender)
```

## 7.4 Conclusion

### 7.4.1 What’s to come?

Congratulations! We’re ready to proceed to the third portion of this book: “statistical inference” using a new package called `infer`

. Once we’ve covered Chapters 8 on sampling, 9 on confidence intervals, and 10 on hypothesis testing, we’ll come back to the models we’ve seen in “data modeling” in Chapter 11 on inference for regression. As we said at the end of Chapter 6, we’ll see why we’ve been conducting the residual analyses from Subsections {#model3residuals} and {#model4residuals}; we are actually verifying some very important assumptions that must be met for the `std_error`

(standard error), `p_value`

, `conf_low`

and `conf_high`

(the end-points of the confidence intervals) columns in our regression tables to have valid interpretation.

### 7.4.2 Script of R code

An R script file of all R code used in this chapter is available here.