In multilinear regression, we assume that the average y-values are a
(affine) linear function of \(p\)
explanatory variables \(x_1, \ldots,
x_p\): \[\mu_y = \beta_0 + \beta_1 x_1
+ \ldots + \beta_p x_p.\] We can use the lm()
command in R to estimate the coefficients in the model.
The model we get will have the form \[\hat{y} = b_0 + b_1 x_1 + \ldots + b_p
x_p.\]
Regression (hyper)planes
If there are two explanatory variables \(x_1\) and \(x_2\), then you get a regression plane
instead of a line.
## place bib gender age state time minutes
## 1 1 66 M 35 VA 1:10:28 70.47
## 2 2 87 M 29 VA 1:18:08 78.13
## 3 3 112 F 32 VA 1:25:47 85.78
## 4 4 116 M 32 VA 1:27:02 87.03
## 5 5 32 M 38 VA 1:27:14 87.23
## 6 6 115 F 31 VA 1:28:15 88.25
myLM <-lm(minutes ~ age + gender, data = results)myLM
##
## Call:
## lm(formula = minutes ~ age + gender, data = results)
##
## Coefficients:
## (Intercept) age genderM
## 105.2479 0.9703 -21.0004
What does plot(linear_model) do?
You can plot a linear model in R. You get several images that help
you check the assumptions for linear regression inference (normal
residuals and constant variance). The par() function lets
you organize the plots into rows & columns.
par(mfrow=c(2,2))plot(myLM)
Example: Course evaluations
What factors explain differences in teaching evaluation scores? Here
is data from 463 courses at the University of Texas - Austin that can be
found in the R library moderndive (from the online textbook
ModernDive)
In a simple linear model, the coefficient on the age
variable does not depend on gender, so the lines for
predicting teaching scores are parallel.
simple_lm <-lm(score ~ age + gender, data = evals_data)simple_lm
##
## Call:
## lm(formula = score ~ age + gender, data = evals_data)
##
## Coefficients:
## (Intercept) age gendermale
## 4.484116 -0.008678 0.190571
Linear models with interaction
A better linear model might have different slopes for men versus
women. You can get models like that by adding an iteraction
term: \[\mu_y = \beta_0 + \beta_1
x_1 + \beta_2 x_2 + \underbrace{\beta_{12} x_1 x_2}_{\text{interaction
term}}.\]
The two trendlines are not parallel when you include interaction
effects. Women appear to have a steeper age penalty on course
evaluations than male instructors do.
Interaction model
Here’s how you make a linear model with interactions in R. Notice
that you use * instead of + in the
lm() function.
interaction_lm <-lm(score ~ age * gender, data = evals_data)interaction_lm
##
## Call:
## lm(formula = score ~ age * gender, data = evals_data)
##
## Coefficients:
## (Intercept) age gendermale age:gendermale
## 4.88299 -0.01752 -0.44604 0.01353
Based on the output above, what is the slope for male instructors?
What is the slope for female instructors?
The babies data frame contains data on baby birth
weights. The variables are:
bwt - birth weight (in ounces)
gestation - length of the pregnancy (in days)
parity - 1 if baby was first born, 0 otherwise
age - mother’s age (in years)
height - mother’s height (in inches)
weight - mother’s weight (in lbs.)
smoke - 1 if the mother is a smoker, 0 otherwise
Cleaning up the data
There are a lot of cells with NA (not available) entries, and these
could mess up our analysis below. The na.omit() command is
a fast way to remove these.
babies <-na.omit(babies)dim(babies)
## [1] 1174 8
Whenever you omit data, you should make sure that you aren’t omitting
a large percentage of the sample, and you might also want to check the
way that the data was collected to make sure that individuals with
missing data are not systematically different from other individuals in
the sample. In this example, we are omitting 62 rows of data (out of
1236). That’s only about 5% of the data, so we probably aren’t affecting
our results too much.
The linear model
babies_lm <-lm(bwt ~ gestation + parity + age + height + weight + smoke, data = babies)babies_lm
We need to check that there is a roughly linear relationship between
each of the explanatory variables and the response variable. This also
lets you see how the residuals depend on each explanatory variable.
par(mfrow=c(2, 3))plot(babies$gestation, babies$bwt, xlab='Gestation time (in days)', ylab='Birth weight (in ounces)')plot(babies$age, babies$bwt, xlab="Mother's age", ylab='Birth weight (in ounces)')plot(babies$height, babies$bwt, xlab="Mother's height (inches)", ylab='Birth weight (in ounces)')plot(babies$weight, babies$bwt, xlab="Mother's weight (lbs.)", ylab='Birth weight (in ounces)')plot(babies$parity, babies$bwt, xlab="Parity", ylab='Birth weight (in ounces)')plot(babies$smoke, babies$bwt, xlab="Smoker", ylab='Birth weight (in ounces)')
What does plot(data_frame) do?
If you plot a data frame, then you get a matrix of pairwise plots
between every pair of variables in the data frame. That can help check
that there is a linear relationship between the variables.
plot(babies)
Checking the residuals
We also need to check the residuals to see that they are roughly
normally distributed with constant variance. This is harder to do with
so many variables. Here are two of the most important things to
check:
Plot the residuals versus predicted \(\hat{y}\) values (residuals versus fitted
data)
A normal quantile plot of residuals (to check for normality)
par(mfrow=c(1,2))qqnorm(resid(babies_lm))qqline(resid(babies_lm))plot(fitted(babies_lm),resid(babies_lm),xlab='Predicted values',ylab='Residuals',main='Residuals vs. Predicted Values')abline(0,0)
You could also get these plots automatically (with some extra junk)
if you plot() the linear model.
Confidence & prediction intervals
These work exactly the same as the single variable case.
With more than two variables, you can include several interactions.
For example,
lm(y ~ x1*x2 + x1*x3 + x2*x3, data = df)
would include both the main effects and also the interaction terms
for all three pairs of variables in the model: \[\mu_y = \beta_0 + \underbrace{\beta_{1} x_1 +
\beta_2 x_2 + \beta_3 x_3}_\text{main effects} + \underbrace{\beta_{12}
x_1 x_2 + \beta_{13} x_1 x_3 + \beta_{23} x_2 x_3}_\text{interaction
terms}.\]
Including interaction terms between two quantitative variables makes
the model nonlinear in the explanatory variables. That isn’t necessarily
a bad things, but it adds complexity.
Practice
A group of students wanted to investigate which factors influence the
price of a house (Koester, Davis, and Ross, 2003). They used http://www.househunt.com, limiting their search to
single family homes in California. They collected a stratified sample by
stratifying on three different regions in CA (northern, southern, and
central), and then randomly selecting a sample from within each strata.
They decided to focus on homes that were less than 5000 square feet and
sold for less than $1.5 million.
