2. Simple linear regression

The Toluca Company builds refrigeration equipment and the replacement parts that keep it running. One part had long been produced in batches, or lots, of whatever size seemed convenient at the time. When the company started a cost-cutting review, its engineers wanted to settle a basic question: what lot size is cheapest to produce? Every production run carries a fixed setup cost that does not care about lot size, plus machining and assembly work that grows with the number of parts. To find the sweet spot, the engineers first needed one ingredient: how do the labor hours a run consumes depend on the size of the lot?

They pulled records from 25 recent runs made under stable conditions over a six-month period. Each record is a pair: the lot size (units produced) and the work hours the run took. The lot sizes are all multiples of ten, a scheduling convenience, and they range from 20 to 120 units. The scatter of these 25 points is the starting picture of the whole chapter.

Two things stand out in Figure 1. Bigger lots take more hours, and the cloud of points is close enough to a straight line that a line is a reasonable summary. But which line? An engineer could lay a ruler on the plot and draw one by eye, and a second engineer would draw a slightly different one. This chapter replaces the ruler with a rule: a precise, defensible recipe for the single best line, an honest account of how much to trust it, and the assumptions that make the whole enterprise valid. By the end you will be able to fit that line, read its slope in the units of the problem (hours per unit of lot size), say how uncertain the slope is, and prove that under stated conditions no competing method does better.

The modeling workflow¶

Every chapter of this book is one pass through the same five-stage loop. Naming the stages now gives you a map to carry through the whole course, and every later chapter opens by saying which stage it serves.

• ASK. Fix the question and the data that can answer it. Here: how do work hours depend on lot size, and what can 25 production runs tell us?

• EXPLORE. Plot and summarize before fitting anything. The scatter in Figure 1 is this stage; it is what tells us a straight line is even worth trying.

• FIT. Estimate the model. Most of this chapter lives here, turning the ruler-and-eye line into the single least-squares line.

• CHECK. Ask whether the fitted model can be trusted, using the residuals and the assumptions behind them. We begin this in Section 2.3 and devote all of Chapter 9 to it.

• USE. Interpret the estimates, predict, and decide. Reading the slope in hours per unit of lot size is the first taste; inference and prediction fill Chapter 3.

As Figure 2 shows, the stages form a loop, not a one-way street: a failed CHECK sends you back to EXPLORE or FIT with a better model. Simple linear regression is the smallest complete trip around this loop, which is why the course starts here.

2.1 The simple linear regression model¶

Intuition¶

Regression starts from a plain idea: one thing tends to move with another, and we want to pin that relationship down in numbers. Here we describe how a response $Y$ (Toluca work hours) depends on a single predictor $X$ (lot size). No sane model claims the hours are an exact function of lot size: two runs of the same size took different amounts of time, because of tool wear, the crew on shift, the weather, and a hundred small things nobody recorded. So we split each observation into two pieces: a systematic part that lot size explains, and a random part that it does not.

The systematic part is a straight line. On average, a run of size $X$ takes some baseline amount of time plus a fixed number of extra hours per additional unit. The random part is a disturbance that pushes the actual hours above or below that average. We assume the disturbance has no built-in tendency to be positive or negative, has the same typical size regardless of lot size, and does not coordinate from one run to the next.

Formula¶

The simple linear regression model (Definition 2.1) for observations $i = 1, \dots, n$ is

• $Y_i$ is the observed response in the $i$-th case (work hours for run $i$).

• $X_i$ is the predictor in the $i$-th case (lot size for run $i$), treated as a known, fixed constant. Some books call this the independent variable or the explanatory variable; this book says predictor throughout.

• $\beta_0$ is the intercept parameter: the mean response when $X = 0$.

• $\beta_1$ is the slope parameter: the change in the mean response for a one-unit increase in $X$.

• $\varepsilon_i$ is the random error in the $i$-th case, the amount by which $Y_i$ departs from its mean.

The parameters $\beta_0$ and $\beta_1$ are fixed numbers we do not know and want to estimate. The model puts four conditions on the errors:

In words: each error averages to zero (the line is right on average), every error has the same variance $\sigma^2$ (constant spread, called homoscedasticity), and distinct errors are uncorrelated (one run’s luck says nothing about another’s). Notice what is not yet assumed: no particular shape for the error distribution. Normality enters only in 2.6 Maximum likelihood under normal errors, and only when we want to build intervals and tests.

Because $X_i$ is a fixed constant and $\beta_0 + \beta_1 X_i$ is therefore also constant, the error conditions transfer directly to the responses:

The quantity $E\{Y_i\}$ is the mean response at level $X_i$. The function $E\{Y\} = \beta_0 + \beta_1 X$ that gives the mean response at any $X$ is the regression function (Definition 2.2). It is a line, and estimating it is the whole game.

Figure 3 is the mental picture to carry through the chapter: a straight mean line, with a distribution of possible $Y$ values stacked at each $X$, every stack the same width.

2.2 Least squares from first principles¶

Intuition¶

We have a cloud of points and want the one line that fits best. “Best” needs a definition. Least squares defines it by penalizing misses: for a candidate line, measure the vertical gap from each point to the line, square the gaps so that overshooting and undershooting both count as errors and large misses hurt disproportionately, and add them up. The best line is the one that makes this total, the sum of squared errors, as small as possible. Figure 4 shows the idea with three candidate lines through the Toluca data.

Formula¶

Write a candidate line with intercept $b_0$ and slope $b_1$. Its sum of squared errors (Definition 2.3) is

• $Y_i - b_0 - b_1 X_i$ is the vertical gap from the point $(X_i, Y_i)$ to the candidate line.

• $Q$ is a function of the two numbers $b_0, b_1$ we get to choose; the data are fixed.

The least-squares estimates $b_0$ and $b_1$ are the values that minimize $Q$. To write them compactly, define the three basic sums of squares and cross-products around the means:

where $\bar{X}$ and $\bar{Y}$ are the sample means of the predictor and response. The next two derivations establish the minimizers, stated here first.

Derivation by calculus¶

Proof (by calculus). $Q(b_0, b_1)$ is a sum of squares, so it is smooth and bounded below by zero. At a minimum, both partial derivatives vanish. Differentiate:

Setting each to zero and dividing by -2 gives the normal equations:

Expand the sums. The first equation becomes $\sum Y_i = n b_0 + b_1 \sum X_i$, and dividing by $n$ gives $\bar{Y} = b_0 + b_1 \bar{X}$, so

Substitute this into the second normal equation $\sum X_i Y_i = b_0 \sum X_i + b_1 \sum X_i^2$. Replacing $b_0$ and collecting the $b_1$ terms,

The left side equals $\sum X_i Y_i - n \bar{X}\bar{Y} = \sum (X_i - \bar X)(Y_i - \bar Y) = S_{xy}$, and the right side’s bracket equals $\sum X_i^2 - n\bar{X}^2 = \sum (X_i - \bar X)^2 = S_{xx}$. Therefore

This is the only stationary point, and it is a minimum by the second-derivatives test for a function of two variables (the same test you would use in calculus for one variable, extended to two). The three second partial derivatives of $Q$ are constants: $Q_{b_0 b_0} = 2n$, $Q_{b_1 b_1} = 2\sum X_i^2$, and $Q_{b_0 b_1} = 2\sum X_i$. Form the test quantity $D = Q_{b_0 b_0}\,Q_{b_1 b_1} - (Q_{b_0 b_1})^2 = 4\big(n\sum X_i^2 - (\sum X_i)^2\big) = 4n\,S_{xx}$. Since $D > 0$ whenever the $X_i$ are not all equal, and $Q_{b_0 b_0} = 2n > 0$, the stationary point is a local minimum. Because $Q$ is a sum of squares (so it is bounded below by zero) and this is its only stationary point, that local minimum is also the global minimum. (Chapter 7 revisits this same fact in matrix form, once you have the tools of Chapter 6.) $\blacksquare$

Derivation by an algebraic identity¶

The calculus argument finds the minimizer but leaves the geometry a little hidden. Here is a second proof that uses no derivatives and shows directly that no line beats the least-squares line. It rests on two facts we will prove in 2.3 Fitted values and the properties of residuals: writing the fitted values as $\hat{Y}_i = b_0 + b_1 X_i$ and the residuals as $e_i = Y_i - \hat{Y}_i$, the residuals satisfy $\sum e_i = 0$ and $\sum X_i e_i = 0$.

Proof (by an algebraic identity). Take any competing line with intercept $a$ and slope $c$. Split the gap from each point to that line by routing it through the least-squares fitted value:

Square and sum. The sum of squared errors for the competing line is

The cross term vanishes:

So for every competing line,

with equality only when $\sum g_i^2 = 0$, that is when $g_i = 0$ for all $i$, which (since the $X_i$ are not all equal) forces $a = b_0$ and $c = b_1$. The least-squares line is the unique minimizer. $\blacksquare$

The two normal equations have a clean picture. Each one is a line in the $(b_0, b_1)$ plane, and the least-squares estimates sit where they cross, at the bottom of the SSE bowl. Figure 5 shows the bowl and Figure 6 shows it from above.

R¶

Fitting the model in R uses lm, which reads a formula response ~ predictor. First read the data and look at it.

toluca <- read.csv("data/toluca.csv")
dim(toluca)
head(toluca, 3)

[1] 25  2
  lotsize hours
1      80   399
2      30   121
3      50   221

Being able to reproduce lm from the raw sums $S_{xx}$ and $S_{xy}$ matters: it means you can compute a regression anywhere, check software output, and, later, extend the idea when there is no lm button for the model you need.

Small datasets are worth computing by hand once, so the formulas stop being black boxes. The file toluca_mini.csv is a six-row slice of the Toluca data, small enough to do on paper. Figure 7 shows its points and their deviations from the means, the raw material of $S_{xx}$ and $S_{xy}$.

Before you take the formulas on trust, spend a minute trying to beat them by hand.

What to notice. No pair of slider settings pushes $Q$ below $54{,}825.5$, and the settings that come closest read 62.5 and 3.57, the least-squares estimates. Try this. Hold the slope at 3.57 and sweep the intercept alone: the sum of residuals crosses zero at exactly the setting where $Q$ bottoms out, which is the first normal equation happening in front of you. Back to Section 2.2.

2.3 Fitted values and the properties of residuals¶

Intuition¶

Once the line is fixed, every observed point splits into two parts. The fitted value $\hat{Y}_i$ (Definition 2.5) is where the line says the point should be; the residual $e_i = Y_i - \hat{Y}_i$ is how far the actual point lands from it. Residuals are the model’s leftovers, the part of $Y$ that lot size did not explain. Figure 9 draws them as the vertical sticks connecting each point to the line. They are not just scenery: least squares forces the residuals to satisfy exact algebraic identities, and those identities are both a proof tool and a practical way to catch a mistake in a fit.

Formula¶

The fitted values and residuals are

• $\hat{Y}_i$ is the estimated mean response at $X_i$, the height of the fitted line there.

• $e_i$ is the residual, the vertical gap from the observed $Y_i$ down to the line.

Do not confuse $e_i$ with the model error $\varepsilon_i = Y_i - (\beta_0 + \beta_1 X_i)$. The error is the gap to the true, unknown line; the residual is the gap to our estimated line. We never see $\varepsilon_i$, but we can always compute $e_i$, and it is our best stand-in. Figure 10 puts the two gaps side by side so the difference is impossible to miss.

Derivation (properties of residuals)¶

Proof. All of these follow from the two normal equations of 2.2 Least squares from first principles.

1. The residuals sum to zero. The first normal equation is exactly $\sum (Y_i - b_0 - b_1 X_i) = \sum e_i = 0$.

2. The predictor is orthogonal to the residuals. The second normal equation is exactly $\sum X_i (Y_i - b_0 - b_1 X_i) = \sum X_i e_i = 0$.

3. The fitted values are orthogonal to the residuals. Using properties 1 and 2, $\sum \hat{Y}_i e_i = \sum (b_0 + b_1 X_i) e_i = b_0 \sum e_i + b_1 \sum X_i e_i = 0$.

4. The fitted values reproduce the total of $Y$. From property 1, $\sum e_i = \sum (Y_i - \hat{Y}_i) = 0$, so $\sum \hat{Y}_i = \sum Y_i$, and dividing by $n$, the mean fitted value equals $\bar{Y}$.

5. The line passes through the center of the data. Setting $X = \bar{X}$ in the fitted line and using $b_0 = \bar{Y} - b_1 \bar{X}$ gives $\hat{Y} = b_0 + b_1 \bar{X} = \bar{Y}$. So $(\bar{X}, \bar{Y})$ is always on the least-squares line. $\blacksquare$

Property 1 is why the vertical sticks in Figure 9 balance out. Properties 1 and 2 were the ingredients that made the cross term vanish in the algebraic proof of 2.2 Least squares from first principles, so the identities repay the loan we took there.

R and Python¶

A plot of the residuals against the predictor is the first picture to draw after any fit. If the model’s constant-variance and straight-line assumptions hold, the residuals should hover around zero with no trend and no fanning, as in Figure 11. A curve would signal that the relationship is not straight; a widening funnel would signal nonconstant variance. Chapter 9 builds a full diagnostic toolkit on this one habit.

The identities in Theorem 2.6 are exact, so the fastest way to believe them is to try to break one.

What to notice. However badly you mangle the data, $\sum e_i$ and $\sum X_i e_i$ stay at zero and the mean of $\hat{Y}$ stays equal to $\bar{Y}$. Try this. Drag one run to the very top of the plot: the residuals become huge, the fit becomes useless, and the identities do not move, which is why they check your arithmetic and never your model. Back to Section 2.3.

2.4 Estimating the error variance¶

Intuition¶

The slope and intercept describe the line, but the model has a third unknown: $\sigma^2$, the variance of the errors, which controls how tightly the points hug the line. A small $\sigma^2$ means predictions can be sharp; a large one means even a perfect line leaves big uncertainty. We cannot see the errors $\varepsilon_i$, but the residuals $e_i$ are their visible stand-ins, so we estimate $\sigma^2$ from the spread of the residuals. Figure 13 shows what that spread looks like as a band around the fitted line.

Formula¶

The error sum of squares (Definition 2.7) and the estimator of $\sigma^2$ are

• $\mathrm{SSE}$ is the total squared residual, the leftover variation the line did not explain.

• $\mathrm{MSE}$ (mean square error), also written $s^2$, is $\mathrm{SSE}$ divided by its degrees of freedom $n-2$; it estimates $\sigma^2$.

• $s$ is the estimated standard deviation of the errors, in the units of $Y$ (hours).

The divisor is $n-2$, not $n$ or $n-1$. Two degrees of freedom are spent estimating the two coefficients $b_0$ and $b_1$ that define the line the residuals are measured from. The next derivation shows that this exact divisor, and no other, makes $\mathrm{MSE}$ correct on average.

Derivation (why the divisor is $n-2$)¶

Proof. We show $E\{\mathrm{SSE}\} = (n-2)\sigma^2$, so that dividing by $n-2$ gives an estimator that is right on average. Start from a tidy identity for SSE. Since $\hat{Y}_i = \bar{Y} + b_1(X_i - \bar{X})$ (using $b_0 = \bar{Y} - b_1 \bar X$),

so, expanding the square and using $S_{xy} = b_1 S_{xx}$,

Take expectations of the two pieces. For any random variable, $E\{W^2\} = \operatorname{Var}\{W\} + (E\{W\})^2$.

The slope piece. We show in 2.5 Sampling behavior and the Gauss-Markov theorem that $E\{b_1\} = \beta_1$ and $\operatorname{Var}\{b_1\} = \sigma^2 / S_{xx}$. Hence

The total piece. Write $S_{yy} = \sum Y_i^2 - n\bar{Y}^2$. Each $Y_i$ has mean $\mu_i = \beta_0 + \beta_1 X_i$ and variance $\sigma^2$, so $E\{Y_i^2\} = \sigma^2 + \mu_i^2$ and $E\{\sum Y_i^2\} = n\sigma^2 + \sum \mu_i^2$. The sample mean $\bar{Y}$ has mean $\bar{\mu} = \beta_0 + \beta_1 \bar{X}$ and variance $\sigma^2/n$, so $E\{n\bar{Y}^2\} = n(\sigma^2/n + \bar{\mu}^2) = \sigma^2 + n\bar{\mu}^2$. Subtracting,

The final bracket is $\sum (\mu_i - \bar{\mu})^2 = \sum \big(\beta_1 (X_i - \bar{X})\big)^2 = \beta_1^2 S_{xx}$. So $E\{S_{yy}\} = (n-1)\sigma^2 + \beta_1^2 S_{xx}$.

Combine. Subtracting the slope piece from the total piece,

Therefore $E\{\mathrm{MSE}\} = E\{\mathrm{SSE}\}/(n-2) = \sigma^2$: the divisor $n-2$ is exactly what makes $\mathrm{MSE}$ an unbiased estimator of $\sigma^2$. $\blacksquare$

R and Python¶

Theorem 2.8 says the divisor $n-2$ is the one that is right on average, and repeated sampling is the way to watch that claim settle.

What to notice. Individual values of $\mathrm{MSE}$ scatter widely, most of them landing anywhere from about 1200 to about 4000, yet their average walks steadily toward 2383.4. Try this. Press Draw 1000 twice and compare the two means: the $n-2$ divisor lands on $\sigma^2$ while the $n$ divisor stays about 8 percent low no matter how many samples you add, which is what bias means. Back to Section 2.4.

2.5 Sampling behavior and the Gauss-Markov theorem¶

Intuition¶

The estimates $b_0$ and $b_1$ are computed from a particular sample of 25 runs. A different 25 runs, made under the same conditions, would give slightly different estimates. So $b_0$ and $b_1$ are themselves random: they have distributions, called sampling distributions, across the repeated samples we could have drawn. Figure 15 makes this concrete by simulating 60 fresh samples at the Toluca lot sizes and fitting each: the fitted lines scatter around the true mean line, and they fan out most where $X$ is far from $\bar{X}$.

Two questions organize this section. Are $b_0$ and $b_1$ correct on average (unbiased)? And how much do they vary from sample to sample (their variances)? The Gauss-Markov theorem then delivers a strong payoff: among all the sensible estimators one might invent, least squares has the smallest variance.

Formula¶

Both estimators are unbiased, and their variances are

• $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$ shrinks when the errors are small ($\sigma^2$ low) or the predictor is spread wide ($S_{xx}$ large). Spreading your $X$ values out buys precision.

• $\operatorname{Var}\{b_0\}$ carries an extra $\bar{X}^2/S_{xx}$ term because the intercept lives far off at $X = 0$, and tilting the line pivots it there.

Figure 16 makes the first fact visible. When the lot sizes are bunched together, a little noise tips the slope wildly; when they are spread wide, the same noise barely moves it.

Replacing the unknown $\sigma^2$ by its estimate $\mathrm{MSE}$ gives the estimated variances $s^2\{b_1\} = \mathrm{MSE}/S_{xx}$ and $s^2\{b_0\} = \mathrm{MSE}\big(1/n + \bar{X}^2/S_{xx}\big)$, whose square roots are the standard errors $s\{b_1\}$ and $s\{b_0\}$ that software prints.

Derivation (unbiasedness and variances)¶

Proof (mean and variance of $b_1$). The key is to write $b_1$ as a weighted sum of the responses. Since $\sum (X_i - \bar{X}) \bar{Y} = 0$,

The weights $k_i$ depend only on the fixed $X$ values and satisfy three facts, each a one-line computation:

Unbiasedness. Using $E\{Y_i\} = \beta_0 + \beta_1 X_i$,

Variance. The $Y_i$ are uncorrelated with common variance $\sigma^2$, so variances of a weighted sum add with squared weights:

Proof (mean and variance of $b_0$). Write $b_0 = \bar{Y} - b_1 \bar{X} = \sum m_i Y_i$ with $m_i = \tfrac{1}{n} - \bar{X} k_i$. Then $\sum m_i = 1$ and $\sum m_i X_i = \bar X - \bar X = 0$, so $E\{b_0\} = \beta_0 \sum m_i + \beta_1 \sum m_i X_i = \beta_0$. For the variance, using $\sum k_i = 0$ and $\sum k_i^2 = 1/S_{xx}$,

Derivation (Gauss-Markov: $b_1$ is BLUE)¶

Proof. Consider any linear estimator $\hat\beta_1^{\ast} = \sum c_i Y_i$ with fixed weights $c_i$. Its expectation is $\sum c_i(\beta_0 + \beta_1 X_i) = \beta_0 \sum c_i + \beta_1 \sum c_i X_i$. For this to equal $\beta_1$ for every possible $\beta_0$ and $\beta_1$, the weights must satisfy

Write $c_i = k_i + d_i$, splitting each weight into the least-squares weight $k_i$ and a deviation $d_i$. Because both $c_i$ and $k_i$ satisfy the two unbiasedness constraints, the deviations satisfy $\sum d_i = 0$ and $\sum d_i X_i = 0$. Now the cross term drops out:

Therefore

Since $\sigma^2 \sum d_i^2 \ge 0$, every linear unbiased competitor has variance at least that of $b_1$, with equality only when all $d_i = 0$, that is when the competitor is $b_1$. The least-squares slope is BLUE. The same argument with weights $m_i$ shows $b_0$ is BLUE. $\blacksquare$

R and Python¶

The variance formulas are not just algebra; they describe the fanning bundle of lines in Figure 15. We can confirm both the unbiasedness and the slope-variance formula by simulation, drawing 5000 samples from the estimated model and refitting each. Figure 17 shows the histogram of the 5000 slope estimates.

set.seed(4210)
beta0 <- 62.366; beta1 <- 3.5702; sigma <- 48.82
xfix <- toluca$lotsize
Sxx_fix <- sum((xfix - mean(xfix))^2)
b1_sim <- replicate(5000, {
  ysim <- beta0 + beta1 * xfix + rnorm(length(xfix), 0, sigma)
  coef(lm(ysim ~ xfix))[["xfix"]]
})
round(c(mean_of_b1_hats = mean(b1_sim),
        true_beta1 = beta1,
        sd_of_b1_hats = sd(b1_sim),
        formula_se = sqrt(sigma^2 / Sxx_fix)), 4)

mean_of_b1_hats      true_beta1   sd_of_b1_hats      formula_se
         3.5699          3.5702          0.3479          0.3469

rng = np.random.default_rng(4210)
beta0, beta1, sigma = 62.366, 3.5702, 48.82
xfix = toluca["lotsize"].to_numpy()
Sxx_fix = np.sum((xfix - xfix.mean()) ** 2)
b1_hats = np.empty(5000)
for i in range(5000):
    ysim = beta0 + beta1 * xfix + rng.normal(0, sigma, size=xfix.size)
    b1_hats[i] = np.polyfit(xfix, ysim, 1)[0]
print(round(b1_hats.mean(), 4), beta1,
      round(b1_hats.std(ddof=1), 4),
      round(np.sqrt(sigma ** 2 / Sxx_fix), 4))

3.5675 3.5702 0.352 0.3469

The average of the 5000 slope estimates sits right on $\beta_1 = 3.5702$ (unbiasedness), and their simulated standard deviation, 0.348 in R and 0.352 in Python, matches the formula value $\sigma/\sqrt{S_{xx}} = 0.347$. The two languages differ only in the last digit because they use different random number generators; the agreement with theory is the point.

The bundle of fitted lines in Figure 15 is easier to believe when you build it one sample at a time.

What to notice. Draw samples one at a time and no single $b_1$ equals 3.5702, yet the pile they build is centered there. Try this. Press Draw 1000 and compare the SD of $b_1$ readout with $\sigma/\sqrt{S_{xx}} = 0.3469$ sitting beside it: the variance formula of Theorem 2.9 is not an approximation, it is the width of this pile. Back to Section 2.5.

2.6 Maximum likelihood under normal errors¶

Intuition¶

So far we chose least squares because squared error is a reasonable penalty, and Gauss-Markov rewarded that choice with the best-linear-unbiased property, all without assuming any shape for the error distribution. Now suppose we add one assumption: the errors are normally distributed. A different principle, maximum likelihood, asks which parameter values make the data we actually observed most probable. Remarkably, under normal errors that principle points at the very same line. Figure 19 shows the two criteria, sum of squared errors and negative log-likelihood, bottoming out at exactly the same slope.

Formula¶

The normal error regression model (Definition 2.11) keeps everything from 2.1 The simple linear regression model and adds a shape:

Under normality the uncorrelated errors are now fully independent. The likelihood is the joint density of the observed $Y_i$ read as a function of the parameters:

• $L$ is large for parameter values that make the observed data look typical, small for values that make it look freakish.

• Maximum likelihood estimates are the $\beta_0, \beta_1, \sigma^2$ that maximize $L$.

Derivation (maximum likelihood equals least squares)¶

Proof. Maximizing $L$ is easier through its logarithm, which is maximized at the same place because $\log$ is increasing:

Coefficients. Fix any $\sigma^2 > 0$. The first two terms do not involve $\beta_0$ or $\beta_1$, and the last term has $-1/(2\sigma^2) < 0$ multiplying $\sum (Y_i - \beta_0 - \beta_1 X_i)^2$. Maximizing $\ln L$ over $\beta_0, \beta_1$ therefore means minimizing that sum of squares, which is exactly the least-squares criterion $Q$ from 2.2 Least squares from first principles. So the maximum-likelihood estimates of the coefficients are the least-squares estimates:

Variance. Substitute $b_0, b_1$ so the sum of squares becomes $\mathrm{SSE}$, then differentiate $\ln L$ with respect to $\sigma^2$ and set it to zero:

The maximum-likelihood variance estimate divides by $n$, not $n-2$. From 2.4 Estimating the error variance we know the $n$ divisor is biased low, understating $\sigma^2$, because it ignores the two degrees of freedom spent on $b_0$ and $b_1$. That is why the book, and all the inference in Chapter 3, uses the unbiased $\mathrm{MSE} = \mathrm{SSE}/(n-2)$ rather than the maximum-likelihood $\hat\sigma^2_{\text{MLE}}$. $\blacksquare$

R and Python¶

2.7 Chapter summary¶

This chapter built a simple linear regression from the ground up and proved the properties that make it trustworthy. Starting from the model $Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$ and its three error assumptions, we defined the least-squares criterion and derived the estimators two independent ways, by calculus through the normal equations and by an algebraic identity that shows no other line has a smaller SSE. The residuals obey exact identities that double as a correctness check; the error variance is estimated by $\mathrm{MSE} = \mathrm{SSE}/(n-2)$, whose $n-2$ divisor is exactly what makes it unbiased. Treating $b_0$ and $b_1$ as random over repeated samples gave their means and variances, the Gauss-Markov theorem showed they are best among linear unbiased estimators with no normality needed, and maximum likelihood under normal errors reproduced the very same line.

Key results at a glance

For the Toluca data these give $b_0 = 62.37$, $b_1 = 3.5702$, $\mathrm{SSE} = 54{,}825$, $\mathrm{MSE} = 2384$, $s = 48.82$, $s\{b_1\} = 0.347$, $s\{b_0\} = 26.18$, and $\hat\sigma^2_{\text{MLE}} = 2193$.

Key terms

Simple linear regression model, regression function, mean response, random error, least-squares criterion, least-squares estimates, normal equations, fitted value, residual, error sum of squares (SSE), mean square error (MSE), homoscedasticity, unbiased estimator, Gauss-Markov theorem, best linear unbiased estimator (BLUE), normal error regression model, maximum likelihood.

You should now be able to

• State the simple linear regression model and explain the job each error assumption does.

• Derive the least-squares estimators two ways, by calculus and by the algebraic identity.

• Compute $b_0$, $b_1$, fitted values, and residuals by hand and with software.

• Prove the residual identities and use them as a built-in check on a fit.

• Explain the $n-2$ divisor by deriving $E\{\mathrm{MSE}\} = \sigma^2$, and compute $\mathrm{MSE}$ and $s$.

• Derive the means and variances of $b_0$ and $b_1$ and prove that $b_1$ is the best linear unbiased estimator.

• Show that maximum likelihood under normal errors reproduces the least-squares estimates.

Where this fits. This chapter is the first full trip around the workflow of The modeling workflow. It carried the Toluca question through ASK and EXPLORE, settled into FIT, and took first steps into CHECK (the residual identities and the residual plot) and USE (reading the slope in hours per unit). Chapter 1 argued that regression is a way of thinking about how one quantity moves with another; this chapter turned that idea into a specific, computable line and proved the properties that make it trustworthy. Everything here is point estimation: single best guesses for $\beta_0$, $\beta_1$, and $\sigma^2$. What we have not yet done is quantify uncertainty with intervals and tests. That is Chapter 3, which takes the sampling distributions sketched here, adds the normal error assumption in full, and builds $t$ intervals for the slope and intercept, the ANOVA decomposition (3.6 The analysis of variance approach), the $F$ test (3.7 The F test and its equivalence to the t test), and the difference between a confidence interval for a mean response and a prediction interval for a new run (3.5 Prediction interval for a new observation). The exact Toluca fit set here ($b_1 = 3.5702$, $\mathrm{MSE} = 2384$) carries forward unchanged, now with inference attached. It travels further still: Chapter 5 reruns the same slope through a permutation test (5.2 The permutation test for the slope) and checks its p-value against Chapter 3’s, and in Chapter 7 the same estimators reappear in matrix form (7.1 The model and least squares in matrix form), where the two-line normal equations become a single equation $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{Y}$ and the Gauss-Markov proof (7.6 The Gauss-Markov theorem) generalizes to any number of predictors at once.

2.8 Frequently asked questions¶

Q1. Why square the residuals instead of taking absolute values? Both are reasonable losses. Squaring makes the objective smooth and differentiable, so calculus gives a clean closed-form solution (the normal equations), and it is exactly the loss for which the Gauss-Markov and normal-likelihood results hold. Absolute-value loss gives “least absolute deviations” regression, which is more resistant to outliers but has no simple formula and needs iterative solving. The course starts with squared error because it is where the theory is cleanest.

Q2. Is the intercept $b_0 = 62.37$ really “hours to produce zero units”? No, not meaningfully. A lot size of zero is far outside the data (the smallest run was 20 units), so the intercept is an extrapolation. Read it as the number that anchors the line’s height, computed so the line passes through $(\bar{X}, \bar{Y})$, not as a physical setup time. Interpreting an intercept literally only makes sense when $X = 0$ is inside or near the observed range.

Q3. What is the difference between an error $\varepsilon_i$ and a residual $e_i$? The error $\varepsilon_i = Y_i - (\beta_0 + \beta_1 X_i)$ is the gap to the true, unknown line, and we can never observe it. The residual $e_i = Y_i - \hat{Y}_i$ is the gap to our estimated line, and we compute it directly. Residuals are our visible estimates of the invisible errors, which is why we estimate $\sigma^2$ from them.

Q4. Why divide SSE by $n-2$ and not $n-1$ like an ordinary sample variance? A sample variance divides by $n-1$ because it estimates one quantity, the mean, from the data. Regression estimates two, the intercept and slope, before measuring residuals, so it spends two degrees of freedom, leaving $n-2$. The derivation in 2.4 Estimating the error variance shows $E\{\mathrm{SSE}\} = (n-2)\sigma^2$ exactly, so $n-2$ is the divisor that makes $\mathrm{MSE}$ unbiased.

Q5. Do I need the errors to be normal for least squares to work? No. The least-squares estimates, their unbiasedness, their variance formulas, and the Gauss-Markov optimality all hold with no assumption about the error distribution’s shape, only mean zero, constant variance, and uncorrelated errors. Normality is an extra assumption, added in 2.6 Maximum likelihood under normal errors and used in Chapter 3 to build exact $t$ and $F$ intervals and tests.

Q6. Why does the bundle of fitted lines in Figure 15 pinch in the middle? Every least-squares line passes through $(\bar{X}, \bar{Y})$ (residual property 5), so near $\bar{X}$ the lines can only wobble a little. Away from $\bar{X}$, a small change in the slope pivots the line a large vertical distance, so the lines fan out. The variance of a fitted value grows with distance from $\bar{X}$, which Chapter 3 turns into the curved shape of a confidence band.

Q7. R prints 4.45e-10 and Python prints 4.45e-10 for the slope’s p-value, but we did not assume normality. Where did that come from? Those p-values and the standard errors’ t-values do use the normal error model, and we will justify them properly in Chapter 3. This chapter shows the standard errors themselves, which come from the variance formulas and need no normality; the significance stars and p-values are a preview of the inference to come.

2.9 Practice problems¶

1. (A) State the simple linear regression model and list its three assumptions on the errors, saying in one phrase what each assumption guarantees.

2. (A) Explain the difference between the regression function $E\{Y\} = \beta_0 + \beta_1 X$ and a fitted line $\hat{Y} = b_0 + b_1 X$. Which involves parameters and which involves estimates?

3. (A) In the Toluca fit, interpret the slope 3.5702 in the units of the problem. Why is it unsafe to use the line to predict hours for a lot of 500 units?

4. (A) Give a one-sentence interpretation of $s = 48.82$ hours for a manager who has never taken statistics.

5. (A) Why is the least-squares intercept for the Toluca data not a meaningful “setup time”? What would have to be true of the data for an intercept to be interpretable?

6. (A) Without computing, state the values of $\sum e_i$ and $\sum X_i e_i$ for any least-squares fit with an intercept, and say which normal equation gives each.

7. (A) A fitted line has $\sum e_i = 0$ but $\sum X_i e_i = 120$. Can it be a least-squares fit? Explain.

8. (A) Explain why spreading the predictor values farther apart (increasing $S_{xx}$) makes the slope estimate more precise, referring to its variance formula.

9. (B) Starting from $Q(b_0, b_1) = \sum (Y_i - b_0 - b_1 X_i)^2$, derive both normal equations by differentiation, and solve them for $b_0$ and $b_1$ (Theorem 2.4).

10. (B) Prove the algebraic identity $\sum (Y_i - a - cX_i)^2 = \mathrm{SSE} + \sum g_i^2$ for an arbitrary line $(a, c)$, identifying $g_i$ and showing the cross term is zero. Conclude that least squares is the unique minimizer.

11. (B) Show that the least-squares line passes through $(\bar{X}, \bar{Y})$, and that $\sum \hat{Y}_i = \sum Y_i$ (Theorem 2.6).

12. (B) Writing $b_1 = \sum k_i Y_i$ with $k_i = (X_i - \bar{X})/S_{xx}$, prove $\sum k_i = 0$, $\sum k_i X_i = 1$, and $\sum k_i^2 = 1/S_{xx}$, then use them to show $E\{b_1\} = \beta_1$ and $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$.

13. (B) Derive $\operatorname{Var}\{b_0\} = \sigma^2(1/n + \bar{X}^2/S_{xx})$ by writing $b_0 = \sum m_i Y_i$ with $m_i = 1/n - \bar{X}k_i$.

14. (B) Prove the Gauss-Markov theorem for the slope (Theorem 2.10): any linear unbiased estimator $\sum c_i Y_i$ of $\beta_1$ has variance at least $\sigma^2/S_{xx}$, with equality only for $b_1$.

15. (B) Derive $E\{\mathrm{SSE}\} = (n-2)\sigma^2$ (Theorem 2.8) using $\mathrm{SSE} = S_{yy} - b_1^2 S_{xx}$ and the facts $E\{b_1\} = \beta_1$, $\operatorname{Var}\{b_1\} = \sigma^2/S_{xx}$.

16. (B) Write down the log-likelihood for the normal error model and show that maximizing it over $\beta_0, \beta_1$ is equivalent to least squares (Theorem 2.12). Then derive $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$.

17. (B) Show that $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$ is biased for $\sigma^2$, and compute the exact bias in terms of $\sigma^2$ and $n$.

18. (B) Prove the decomposition $\mathrm{SSE} = S_{yy} - b_1^2 S_{xx}$ starting from $e_i = (Y_i - \bar{Y}) - b_1(X_i - \bar{X})$.

19. (B) A student proposes estimating the slope by the ratio of means $\bar{Y}/\bar{X}$ (a line through the origin). Show this is generally biased for $\beta_1$, and explain which model assumption it implicitly, and wrongly, imposes.

20. (C) Read toluca.csv in R or Python, fit the model, and reproduce $b_0$, $b_1$, and $s$. Confirm $b_1 = S_{xy}/S_{xx}$ from the raw sums.

21. (C) Compute the residuals in software and verify numerically that $\sum e_i = 0$, $\sum X_i e_i = 0$, and $\sum \hat{Y}_i = \sum Y_i$. Report the three sums.

22. (C) Using toluca_mini.csv, compute $b_0$ and $b_1$ by hand (show the deviation table), then reproduce them in software. Explain why they differ from the full-data estimates.

23. (C) Predict the mean work hours for lot sizes 40, 80, and 120 using the fitted line, and mark which of these are interpolation and which (if any) are extrapolation.

24. (C) Fit the model, extract $\mathrm{MSE}$ two ways (from summary / fit.summary() and by hand as $\mathrm{SSE}/(n-2)$), and confirm they match. Then compute $\hat\sigma^2_{\text{MLE}} = \mathrm{SSE}/n$ and report the percentage difference.

25. (C) Write a short simulation (seed 4210) that draws 2000 samples from the fitted model at the Toluca lot sizes, refits each, and reports the mean and standard deviation of the 2000 intercept estimates $b_0$. Compare the standard deviation to the formula $s\{b_0\}$.

26. (C) Make the residuals-versus-lot-size plot in R or Python. Describe what a violation of the constant-variance assumption would look like in this plot, and say whether the Toluca residuals show it.

27. (C) Reverse the roles: regress lotsize on hours. Show that the new slope is not the reciprocal of the original slope, and explain why regressing $Y$ on $X$ and $X$ on $Y$ give different lines.

28. (B) An analyst has only the summary numbers $\bar{X} = 70$, $\bar{Y} = 312.28$, $S_{xx} = 19{,}800$, and $S_{xy} = 70{,}690$, with the raw data lost. Can they recover $b_0$ and $b_1$? Can they recover $\mathrm{MSE}$? Explain what is and is not reconstructable from these four numbers.