15  Simple Regression

We begin with a simple bivariate (two variable) model with one outcome and one predictor. The reason we start with this basic model is to help new users understand the output and details without adding too many other details into the mix. So, we start with the easiest model - the bivariate regression model.

15.1 Learning Objectives

  1. Understand the purpose of the model.
  2. Learn how to interpret the output of a bivariate regression model.
  3. Understand how to use the model to make predictions.
  4. Learn how to evaluate the model’s performance.

15.2 A Little History (You Can Skip It, But Don’t)

Regression is old. Adrien-Marie Legendre published the method of least squares in 1805, and Carl Friedrich Gauss claimed he had been using it since 1795 (priority disputes like this are common in the history of mathematics). But the word regression comes from Sir Francis Galton, who in 1885 noticed that the tall sons of tall fathers tended to be a little shorter than their fathers, and the short sons of short fathers a little taller, both regressing toward the mean. Yule (1897) and Pearson (1903) took Galton’s idea and gave us the machinery we still use today.

Why do we tell you this? Because regression was invented to answer an honest, human question: if I know something about you (X), what is my best guess about something else (Y)? That is the whole game. Everything below is just bookkeeping in service of that one question.

15.3 The Model

A regression model is a claim that one variable is a function of another:

\[Y = f(X)\]

For the simple linear case, that function is a straight line:

\[y = bx + a\]

where \(b\) is the slope (how much \(y\) changes for a one-unit change in \(x\)) and \(a\) is the intercept (the value of \(y\) when \(x = 0\)). You may have seen this in algebra class as \(y = mx + b\); statisticians use different letters for the same idea.

Let’s make some data and look at it. We will keep it tiny - five points - so you can check every number by hand.

x <- c(1, 2, 3, 4, 5)
y <- c(4, 8, 10, 7, 8)

plot(x, y, xlim = c(0, 6), pch = 19,
     main = "Five points and the line that fits them best")

Now we fit the model. In R, the workhorse is lm() (for linear model). The formula y ~ x reads “y is predicted by x.”

lm1 <- lm(y ~ x)
plot(x, y, xlim = c(0, 6), pch = 19)
abline(lm1, col = "red", lwd = 2)

summary(lm1)

Call:
lm(formula = y ~ x)

Residuals:
   1    2    3    4    5 
-2.0  1.3  2.6 -1.1 -0.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.3000     2.2898   2.315    0.104
x             0.7000     0.6904   1.014    0.385

Residual standard error: 2.183 on 3 degrees of freedom
Multiple R-squared:  0.2552,    Adjusted R-squared:  0.006944 
F-statistic: 1.028 on 1 and 3 DF,  p-value: 0.3853
REGRESSION
  /STATISTICS COEFF R ANOVA
  /DEPENDENT y
  /METHOD=ENTER x.
using GLM
m1 = lm(@formula(y ~ x), df)
coeftable(m1)

R just handed us a slope and an intercept. Here is what matters most in this book: do not trust the machine until you can reproduce what it did. The rest of this chapter rebuilds every number in that output by hand. If you can do that, you understand regression, rather than just knowing how to type summary().

15.4 What “Best” Means: Residuals

The red line is “best” in a very specific sense. For each point, the residual is the vertical distance between what we observed (\(y\)) and what the line predicts (\(\hat{y}\)):

\[e = y - \hat{y}\]

Let’s draw those distances.

plot(x, y, xlim = c(0, 6), pch = 19,
     main = "Residuals: the vertical gaps the line leaves behind")
abline(lm1, col = "red", lwd = 2)
for (i in 1:5) {
  segments(x[i], y[i], x[i], coef(lm1)[[2]] * x[i] + coef(lm1)[[1]],
           col = "blue", lwd = 2)
}

Here are the residuals themselves:

resid(lm1)
   1    2    3    4    5 
-2.0  1.3  2.6 -1.1 -0.8

Two facts worth burning into memory:

sum(resid(lm1))        # the residuals sum to (essentially) zero
[1] 7.771561e-16
sum(resid(lm1)^2)      # the sum of SQUARED residuals - this is what we minimize
[1] 14.3
ImportantLeast Squares, Defined

The “best” line is the one line (out of infinitely many possible lines) that makes sum(resid^2) as small as it can possibly be. That is all “least squares” means: least (smallest) squares (of the residuals). The residuals sum to zero for any line through the means, so we cannot minimize their plain sum - it is always zero. Squaring first solves that, and it punishes big misses more than small ones. We shall return to that squared quantity again and again; it is central to everything that follows.

15.5 Solving for the Slope and Intercept by Hand

R did not use magic. It used two formulas. Here is the slope:

\[b = \frac{N \sum xy - \sum x \sum y}{N \sum x^2 - \left(\sum x\right)^2}\]

N <- length(x)
b <- (N * sum(x * y) - sum(x) * sum(y)) / (N * sum(x^2) - sum(x)^2)
b
[1] 0.7

Compare that b to the slope in the summary(lm1) output above. Same number. Good.

Once we have the slope, the intercept falls right out. The least-squares line always passes through the point \((\bar{x}, \bar{y})\) - the “center of mass” of the data - so we can substitute the means and solve:

\[a = \bar{y} - b\bar{x}\]

a <- mean(y) - b * mean(x)
a
[1] 5.3

There is a more cryptic-looking version that is algebraically identical:

a <- (sum(y) - b * sum(x)) / N
a
[1] 5.3
NoteWhy are those two the same?

Because \(\bar{y} = \frac{\sum y}{N}\) and \(\bar{x} = \frac{\sum x}{N}\). Substitute those into \(a = \bar{y} - b\bar{x}\), put everything over the common denominator \(N\), and you land on the second formula. Try it with a pencil. Do it once and these “scary” formulas lose their power to intimidate, because you will see they are usually just a tidy version of something simple.

15.6 Quantifying Our Uncertainty: Standard Errors

We have a slope and an intercept. But this is a statistics book, and the entire point of statistics is to quantify our uncertainty. How wrong might these estimates be? For that we need standard errors, and they all start from one quantity: the standard error of estimate (also called \(S_{xy}\) or \(S_{ee}\)), which is basically the typical size of a residual:

\[S_{xy} = \sqrt{\frac{\sum (Y_{obs} - Y_{exp})^2}{N - 2}}\]

We divide by \(N - 2\) (not \(N\)) because we spent two degrees of freedom estimating two things: the slope and the intercept.

Sxy <- sqrt(sum((y - predict(lm1))^2) / (N - 2))
Sxy
[1] 2.18327

That is the Residual standard error line in the summary() output. From it we get the standard error of the slope and of the intercept:

# standard error of b (the slope)
Sb <- Sxy / sqrt(sum((x - mean(x))^2))
Sb
[1] 0.6904105
# standard error of a (the intercept)
Sa <- Sxy * sqrt(sum(x^2) / (N * sum((x - mean(x))^2)))
Sa
[1] 2.289833

Now the payoff. A \(t\) statistic is just an estimate divided by its own standard error, and we compare it to a \(t\) distribution with \(N-2\) degrees of freedom:

t_b <- b / Sb
p_b <- 2 * pt(-abs(t_b), df = N - 2)
c(slope = b, SE = Sb, t = t_b, p = p_b)
    slope        SE         t         p 
0.7000000 0.6904105 1.0138896 0.3852987

Look back at the x row of summary(lm1). The estimate, the standard error, the \(t\) value, the \(p\) value: we just rebuilt the entire regression table from scratch. That is the main lesson of the chapter: R is a fast calculator, not an oracle, and you now understand every number it produced.

NoteFrom the Mailbag (a real question from PSYC 754)

“When calculating the standard error of estimate and the standard error of \(b\), why are the calculations for the degrees of freedom different? (I.e., standard error of estimate \(N-2\) but standard error of \(b\) \(N-1\).)”

Sharp catch, and the honest answer is that they should not be different. There is only one rule here: the degrees of freedom equal \(N\) minus the number of parameters you estimated. In simple regression you estimate two (a slope and an intercept), so every standard error in this chapter - the estimate, the slope, the intercept - rests on the same \(N-2\). If a formula sheet ever shows \(N-1\) for the SE of \(b\), that is a slip, not a second rule. (When we add predictors in the next chapter, that same one rule becomes \(N-k-1\).)

NoteAnother route to the slope’s uncertainty: bootstrap it

The standard-error formula above assumes normal, equal-variance errors. When you doubt that, you can skip the formula entirely and bootstrap the slope: resample your rows with replacement a few thousand times, refit lm() each time, and read the SE and a 95% interval straight off the spread of the resampled slopes - no distributional assumption required. Same slope, an honest interval. See Beyond the Normal Curve.

15.7 The Meaning of the Slope (Mind Your Units)

Students memorize “b is the change in y per unit change in x” and then promptly forget to interpret it. Units are your friend here. Suppose \(x\) is height in inches and \(y\) is weight in pounds:

  • \(b\) is measured in lbs per inch (rise over run: \(\frac{\text{lbs}}{\text{inches}}\))
  • \(a\) is measured in lbs
  • so the equation reads: \(y\ (\text{lbs}) = \frac{\text{lbs}}{\text{inches}} \times \text{inches} + \text{lbs}\), i.e. \(\text{lbs} = \text{lbs} + \text{lbs}\) ✓

The units have to balance, and checking them is a free error-detector. If your slope came out in the wrong units, you did something wrong. This is a habit engineers use constantly, and it is worth adopting.

15.8 Standardized Coefficients

Sometimes we want a slope that does not depend on the units at all - so we can compare predictors measured on wildly different scales. We get it by standardizing both variables (converting each to a z-score with scale()) before fitting. When both variables are standardized, the intercept is zero by construction, so we drop it with - 1:

yz <- scale(y)
xz <- scale(x)
lm2 <- lm(yz ~ xz - 1)
summary(lm2)

Call:
lm(formula = yz ~ xz - 1)

Residuals:
      1       2       3       4       5 
-0.9129  0.5934  1.1867 -0.5021 -0.3651 
attr(,"scaled:center")
[1] 7.4
attr(,"scaled:scale")
[1] 2.191

Coefficients:
   Estimate Std. Error t value Pr(>|t|)
xz   0.5052     0.4315   1.171    0.307

Residual standard error: 0.863 on 4 degrees of freedom
Multiple R-squared:  0.2552,    Adjusted R-squared:  0.06901 
F-statistic: 1.371 on 1 and 4 DF,  p-value: 0.3067
* SPSS prints the standardized Beta in the coefficients table by default,
* so the ordinary regression yields the standardized slope with no extra steps.
REGRESSION
  /STATISTICS COEFF R ANOVA
  /DEPENDENT y
  /METHOD=ENTER x.
using GLM, Statistics
z(v) = (v .- mean(v)) ./ std(v)
df2 = DataFrame(y = z(df.y), x = z(df.x))
lm2 = lm(@formula(y ~ x - 1), df2)
coeftable(lm2)

This standardized slope is the beta weight (\(\beta\)). In the bivariate case it equals the correlation between \(x\) and \(y\) - a fact worth remembering, and one we will lean on heavily when we get to multiple predictors (Aiken 2002).

15.9 A Real Example

Five points are great for hand calculation but a little silly. Let’s use a real dataset - one of ours, from years of teaching this material: the high-school and college GPAs (plus SAT scores and letter-of-recommendation quality) of 100 students. Does high-school GPA predict college GPA?

gpa <- read.csv("data/gpa.csv")
str(gpa)
'data.frame':   100 obs. of  4 variables:
 $ CollegeGPA: num  2.04 2.56 3.75 1.1 3 0.05 1.38 1.5 1.38 4.01 ...
 $ HSGPA     : num  2.01 3.4 3.68 1.54 3.32 0.33 0.36 1.97 2.03 2.05 ...
 $ SATtotal  : int  1070 1254 1466 706 1160 756 1058 1008 1104 1200 ...
 $ QLOR      : int  5 6 6 4 5 3 2 7 4 7 ...
lm3 <- lm(CollegeGPA ~ HSGPA, data = gpa)

plot(CollegeGPA ~ HSGPA, data = gpa, pch = 19, col = "grey40",
     main = "Predicting college GPA from high-school GPA")
abline(lm3, col = "red", lwd = 2)

summary(lm3)

Call:
lm(formula = CollegeGPA ~ HSGPA, data = gpa)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.48231 -0.41083 -0.00236  0.37232  2.02871 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.82204    0.19069   4.311 3.88e-05 ***
HSGPA        0.56549    0.08783   6.438 4.49e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6313 on 98 degrees of freedom
Multiple R-squared:  0.2972,    Adjusted R-squared:  0.2901 
F-statistic: 41.45 on 1 and 98 DF,  p-value: 4.492e-09
REGRESSION
  /STATISTICS COEFF R ANOVA
  /DEPENDENT CollegeGPA
  /METHOD=ENTER HSGPA.
using GLM
lm3 = lm(@formula(CollegeGPA ~ HSGPA), df)
coeftable(lm3)

Interpret it out loud: for every one-point increase in high-school GPA, college GPA is predicted to rise by the slope you see above (in GPA-points per GPA-point). Real data are messier than the toy five points - the scatter is wider, the fit imperfect - but the machinery is identical, and you now know how to read every number in that summary, and how uncertain each one is. (Hold onto this dataset; we bring back its other columns - SAT and recommendation quality - when we add predictors in the next chapter.)

15.10 Challenge

TipDo One Yourself

Take the five-point toy data (x and y) from the top of this chapter and, without calling summary(), calculate by hand (in R, using the formulas above):

  1. the slope \(b\) and intercept \(a\),
  2. the standard error of each,
  3. the \(t\) statistic for the slope, and
  4. the \(p\) value associated with that \(t\).

Then run summary(lm(y ~ x)) and confirm every number matches. Be sure you can interpret each value in a sentence a non-statistician would understand. If you can do that, you have mastered simple regression and you are ready for more than one predictor.

15.11 Where We Go Next

Everything here generalizes. Add a second predictor and \(y = bx + a\) becomes \(y = b_1x_1 + b_2x_2 + a\) - the plane instead of the line. The logic (minimize the squared residuals, then quantify our uncertainty about each coefficient) does not change one bit. That is the subject of the next chapter, Multiple Regression. For a deeper, book-length treatment of everything we introduce here and there, the classic reference is Cohen, Cohen, West, and Aiken (Aiken 2002); for a gentler modern companion, see Lewis-Beck and Lewis-Beck (Lewis-Beck and Lewis-Beck 2024).

Aiken, Stephen G. West, Patricia Cohen. 2002. Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences. 3rd ed. Routledge. https://doi.org/10.4324/9780203774441.
Lewis-Beck, Colin, and Michael Lewis-Beck. 2024. Applied Regression: An Introduction. Second Edition. Thousand Oaks, California. https://doi.org/10.4135/9781483396774.