Linear Regression

Linear regression fits a straight line (or, with several inputs, a flat surface) through data so you can describe how one variable changes with others and make predictions. How many miles per gallon do you lose for every extra 1,000 lbs of car weight? How much does ozone rise per degree of temperature? Regression answers these with numbers and tells you how sure you can be.

In R, one function does the work: lm(), for linear model. In this lesson you will fit a model, read its summary, predict new values, add more variables and check whether the model is trustworthy.

Fitting a model with lm()

lm() takes a formula, response ~ predictor, read as "mpg explained by weight", and the data frame:

model <- lm(mpg ~ wt, data = mtcars)
model
coef(model)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Coefficients:
(Intercept)           wt  
     37.285       -5.344  

(Intercept)          wt 
  37.285126   -5.344472 

The fitted line is mpg = 37.29 − 5.34 × wt. The slope says each extra 1,000 lbs of weight goes with about 5.3 fewer miles per gallon. The intercept (37.29) is the predicted mpg at weight zero; it is needed to place the line but has no real meaning here, because no car weighs nothing. coef() returns the two numbers as a named vector. This is the same red line you drew with abline() in Base Graphics.

Reading summary()

summary() gives the full report:

model <- lm(mpg ~ wt, data = mtcars)
summary(model)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,	Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

The important parts:

• Residuals: the differences between the actual and predicted mpg. They range from about −4.5 to 6.9 and their median is near 0, which is good.
• Coefficients table: each estimate with its standard error (its uncertainty), a t value and a p-value (Pr(>|t|)). The p-value for wt tests H0 "the true slope is 0". At 1.29e-10 (0.000000000129), weight clearly matters. The stars are a quick visual guide to the p-values.
• Residual standard error: predictions are typically about 3 mpg off.
• R-squared: 0.75, so weight explains about 75% of the variation in mpg between cars. Adjusted R-squared penalises extra predictors and is better for comparing models of different sizes.
• F-statistic: tests whether the model as a whole is better than just using the average.

Making predictions

predict() applies the model to new data. The new data frame must have a column with the same name as the predictor:

model <- lm(mpg ~ wt, data = mtcars)
new_cars <- data.frame(wt = c(2.5, 3.5))
predict(model, new_cars)
predict(model, new_cars, interval = "confidence")
round(summary(model)$r.squared, 3)

       1        2 
23.92395 18.57948 
       fit      lwr      upr
1 23.92395 22.55284 25.29506
2 18.57948 17.43342 19.72553
[1] 0.753

A car weighing 2,500 lbs is predicted to do about 23.9 mpg. interval = "confidence" adds a 95% range for the average mpg of all cars of that weight. Use interval = "prediction" for the range of a single new car, which is wider because individual cars vary. Only predict inside the range of your data (here about 1,500 to 5,400 lbs); a straight line can give nonsense far outside it.

More predictors and factors

Add predictors with +. Wrap a categorical variable in factor() so R treats it as groups, not numbers:

model1 <- lm(mpg ~ wt, data = mtcars)
model2 <- lm(mpg ~ wt + hp + factor(am), data = mtcars)
round(coef(summary(model2)), 4)
round(summary(model2)$adj.r.squared, 3)
anova(model1, model2)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  34.0029     2.6427 12.8669   0.0000
wt           -2.8786     0.9050 -3.1808   0.0036
hp           -0.0375     0.0096 -3.9018   0.0005
factor(am)1   2.0837     1.3764  1.5139   0.1413
[1] 0.823
Analysis of Variance Table

Model 1: mpg ~ wt
Model 2: mpg ~ wt + hp + factor(am)
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1     30 278.32                                
2     28 180.29  2    98.031 7.6123 0.002291 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• In a multiple regression each coefficient is the effect of that variable holding the others constant. With horsepower in the model, the weight effect shrinks to −2.9 mpg per 1,000 lbs, because heavy cars also tend to be powerful.
• factor(am)1 is the difference between manual (1) and the reference group, automatic (0): about 2.1 mpg more, but with p = 0.14 that difference is not significant once weight and power are known. The big gap you found with a t-test in Hypothesis Tests is mostly explained by weight and power.
• Adjusted R-squared rises from 0.745 to 0.823.
• anova(model1, model2) tests whether the extra predictors improve the fit; p = 0.002 says they do.

Checking the model

A regression is only trustworthy if its residuals behave: no clear pattern, roughly constant spread. Plot residuals against fitted values:

model <- lm(mpg ~ wt, data = mtcars)
head(round(residuals(model), 2))
plot(model, which = 1)

        Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
            -2.28             -0.92             -2.09              1.30 
Hornet Sportabout           Valiant 
            -0.20             -0.69 

The chart shows residuals on the vertical axis against fitted mpg, with a red smooth line through them. Here the red line bends into a gentle U shape, positive at both ends and negative in the middle, and three cars are labelled as the largest residuals. The U hints that the true relationship curves slightly, which a model with I(wt^2) or log(wt) could capture. plot(model) without which shows four diagnostic plots in turn.

Missing values are dropped automatically. In airquality, 37 days have no ozone reading, so the model uses 116 days, which nobs() reports:

model <- lm(Ozone ~ Temp + Wind, data = airquality)
nobs(model)
round(coef(model), 2)

[1] 116
(Intercept)        Temp        Wind 
     -71.03        1.84       -3.06 

Each extra degree Fahrenheit goes with 1.84 more units of ozone and each extra mph of wind with 3.06 fewer, holding the other constant. Wind disperses pollution; heat helps form it.

Common mistakes

• Forgetting factor() on coded categories such as cyl or am, so R fits a slope through the codes.
• Reading the intercept literally when zero is outside the data.
• Predicting far outside the data range.
• Judging a model by R-squared alone. Check residual plots and whether the coefficients make sense.
• Wrong column names in predict(). The new data needs the exact predictor names.

Key points

• lm(y ~ x, data = df) fits a linear model; coef() gives the intercept and slopes.
• summary() shows estimates, standard errors, p-values, R-squared and the residual standard error.
• predict(model, newdata, interval = ...) gives predictions with confidence or prediction intervals.
• Add predictors with +; wrap categories in factor(); compare nested models with anova().
• Check residual plots with plot(model); coefficients show association, not causation.

m1 <- lm(Petal.Width ~ Petal.Length, data = iris)
round(coef(m1), 3)
round(summary(m1)$r.squared, 3)
predict(m1, data.frame(Petal.Length = c(1.5, 4, 6)))
m2 <- lm(Petal.Width ~ Petal.Length + Species, data = iris)
anova(m1, m2)$`Pr(>F)`[2]

 (Intercept) Petal.Length 
      -0.363        0.416 
[1] 0.927
        1         2         3 
0.2605576 1.2999461 2.1314570 
[1] 5.48225e-10