4.1 Linear regression

A regression line is the straight line closest to the points. Its slope says how much the outcome changes per unit of the predictor, and its uncertainty says how much another sample could change that.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from statsmodels.regression.linear_model import RegressionResults

from statly import style

style.use()
rng = np.random.default_rng(1)

engel = (
    sm.datasets.engel.load_pandas().data
)  # income and food spending, Belgian households (1857)
states = sm.datasets.statecrime.load_pandas().data.drop("District of Columbia")  # 50 US states
income, food = engel.income.to_numpy(), engel.foodexp.to_numpy()


def fit_food_on_income(data: pd.DataFrame) -> RegressionResults:
    """Fit a straight line: food spending against household income."""
    return sm.OLS(data.foodexp, sm.add_constant(data.income)).fit()


fit = fit_food_on_income(engel)

4.1.1 The line of least squares

Chapter 1.3 wrote a model as $Y = f(X) + \varepsilon$. Linear regression takes the simplest recipe for $f$, a straight line:

$$Y = a + b \cdot X + \varepsilon$$

with an intercept $a$ (the height of the line at $X = 0$) and a slope $b$ (how much $Y$ changes when $X$ goes up by one). Which line? Every line leaves gaps between itself and the points, the residuals. Least squares picks the line that makes the sum of the squared gaps as small as possible.

def plot_residuals(ax: Axes, shown: int = 40) -> None:
    """The data, the fitted line, and the gaps from some of the points to it."""
    some = rng.choice(len(income), shown, replace=False)
    fitted = fit.fittedvalues.to_numpy()
    ax.vlines(income[some], fitted[some], food[some], color=style.GREY, lw=0.9)
    ax.scatter(income, food, s=9, color=style.BLUE, alpha=0.5)
    grid = np.array([0, 2600])
    ax.plot(grid, fit.params.const + fit.params.income * grid, color=style.CORAL)
    ax.set(xlabel="household income", ylabel="food spending", xlim=(0, 2600), ylim=(0, 1600))
    ax.set_title("The line and its residuals")


def squared_error(slope: float) -> float:
    """Sum of squared residuals for a line with this slope and its best intercept."""
    intercept = food.mean() - slope * income.mean()
    return float(((food - intercept - slope * income) ** 2).sum())


def plot_valley(ax: Axes) -> None:
    """The sum of squared residuals across slopes."""
    slopes = np.linspace(0.3, 0.75, 200)
    ax.plot(slopes, [squared_error(s) / 1e6 for s in slopes], color=style.INK, lw=1.5)
    best = fit.params.income
    ax.plot(best, squared_error(best) / 1e6, "o", ms=7, color=style.CORAL)
    ax.annotate(
        f"slope {best:.2f}",
        (best, squared_error(best) / 1e6),
        xytext=(8, 10),
        textcoords="offset points",
        color=style.CORAL,
        fontsize=9.5,
    )
    ax.set(
        title="Least squares",
        xlabel="slope of the line",
        ylabel="sum of squared residuals (millions)",
    )


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_residuals(axes[0])
plot_valley(axes[1])
fig.tight_layout()
plt.show()

Figure 4.1.1: Left: food spending against income for the 235 households, with the least squares line and the residuals of a few households (grey). Right: the sum of squared residuals for lines of every slope (each with its best intercept); the fitted slope sits at the bottom of the valley.

Why squares? They punish large gaps more than small ones, they give one clear answer that a formula can find, and when the noise is normal, least squares is the same as maximum likelihood (Chapter 2.3). The price is that one far-away point can pull the line a lot; Chapter 4.3 shows how to spot that.

The bottom of the valley has a formula. The best slope is the covariance of x and y divided by the variance of x, and the best line always passes through the point of the two means:

$$\hat b = \frac{\sum (x_i - \bar x)(y_i - \bar y)}{\sum (x_i - \bar x)^2} = r \, \frac{s_y}{s_x}, \qquad \hat a = \bar y - \hat b \, \bar x$$

The second form shows how slope and correlation relate: the slope is the correlation, rescaled from standard deviations into the units of x and y. The hats mark estimates from the data, as opposed to the true a and b we never see.

Squares or absolute distances?

Chapter 1.2 showed that the mean makes the total squared distance smallest and the median the total absolute distance. The same holds for lines. Least squares draws a mean as a line; least absolute deviations (also called median regression) draws a median as a line, by making $\sum_i |y_i - a - b x_i|$ smallest. Which one estimates the true slope better? A simulation answers: draw many samples from a known line, fit both, and see how the slopes spread.

TRUE_SLOPE, N_POINTS, N_SAMPLES = 0.5, 30, 2000


def simulated_slopes(noise: str) -> tuple[np.ndarray, np.ndarray]:
    """Least squares and least absolute deviation slopes for many simulated samples."""
    squares, absolute = [], []
    for _ in range(N_SAMPLES):
        x = rng.uniform(0, 10, N_POINTS)
        e = rng.normal(0, 1, N_POINTS) if noise == "normal" else rng.standard_t(2, N_POINTS)
        X = sm.add_constant(x)
        y = 1 + TRUE_SLOPE * x + e
        squares.append(sm.OLS(y, X).fit().params[1])
        absolute.append(sm.QuantReg(y, X).fit(q=0.5).params[1])
    return np.array(squares), np.array(absolute)


def plot_slopes(ax: Axes, noise: str, title: str) -> None:
    """Histograms of the two kinds of slope, with their spreads."""
    squares, absolute = simulated_slopes(noise)
    bins = np.linspace(0.1, 0.9, 50)
    style.hist(ax, squares, bins=bins, color=style.BLUE, label="squares")
    style.hist(ax, absolute, bins=bins, color=style.CORAL, label="absolute")
    ax.axvline(TRUE_SLOPE, color=style.INK, ls="--", lw=1)
    style.note(
        ax,
        f"spread (sd) of slopes\nsquares {squares.std():.3f}\nabsolute {absolute.std():.3f}",
        x=0.03,
        y=0.95,
        color=style.INK,
    )
    ax.legend(fontsize=8, loc="upper right")
    ax.set(title=title, xlabel="fitted slope", ylabel="samples")


fig, axes = style.subplots(2, height=3.2, width=10)
plot_slopes(axes[0], "normal", "Normal noise")
plot_slopes(axes[1], "heavy", "Heavy-tailed noise")
fig.tight_layout()
plt.show()

Figure 4.1.2: Slopes from 2,000 simulated samples of 30 points on a line with true slope 0.5 (dashed), fitted by least squares (blue) and by least absolute deviations (red). Left: normal noise; least squares slopes cluster more tightly around the truth. Right: noise with heavy tails (a t distribution with 2 degrees of freedom, which throws occasional wild points); now the absolute-distance line is the steadier one.

With well-behaved noise, squares win, which is one reason nearly every regression is least squares. With wild points, absolute distances (or the robust regression of Chapter 4.3) are safer. Quantile regression extends the idea: instead of the median, it can draw a line through the 10th or 90th percentile of y at each x.

4.1.2 Reading the regression table

Every statistics program prints the same table for a fitted line:

table = fit.summary2().tables[1]
table.index = ["intercept", "income (slope)"]
table.columns = ["estimate", "standard error", "t", "p-value", "95% from", "95% to"]
table = table.round(3).astype(object)
table["p-value"] = ["< 0.001" if p < 0.001 else f"{p:.3f}" for p in fit.pvalues]
print(table.to_markdown(colalign=["left"] + ["right"] * 6))

Table 4.1.1: Food spending regressed on income: the regression table.

	estimate	standard error	t	p-value	95% from	95% to
intercept	147.475	15.957	9.242	< 0.001	116.037	178.914
income (slope)	0.485	0.014	33.772	< 0.001	0.457	0.513

• estimate: the intercept and the slope. Each extra unit of income goes with 0.49 more spent on food: of every extra 100, about 49 went to food. Ernst Engel noticed in 1857 that this share falls as income rises, now called Engel’s law.
• standard error: how much the estimate would vary from sample to sample (Chapter 3.2).
• t and p-value: the test of “the true value is 0”, the estimate divided by its standard error (Chapter 3.4). For the slope, t = 34: income and food spending are clearly related.
• 95% interval: the slope lies between 0.46 and 0.51, with the usual 95% confidence.

The intercept, 147, is the predicted spending at an income of 0. No household has an income near 0, so the number is an extrapolation and means little here. That is common: the intercept mainly puts the line at the right height.

4.1.3 Where the p-value of a slope comes from

The slope’s p-value answers: if income had nothing to do with food spending (the true slope is 0), how often would a sample give a slope at least this far from 0? The steps are the same as for every test in Part 3.

1. Estimate the noise. The residuals $e_i = y_i - \hat y_i$ have spread $\hat\sigma = \sqrt{\sum e_i^2 / (n - 2)}$; two degrees of freedom were used up by the intercept and the slope.

2. Turn it into a standard error for the slope:

   $$\text{SE}(\hat b) = \frac{\hat\sigma}{\sqrt{\sum (x_i - \bar x)^2}}$$

3. Divide: $t = \hat b / \text{SE}(\hat b)$, the slope in standard errors from 0.

4. Read the p-value from the t distribution with n − 2 degrees of freedom.

For the food data t = 34, so far out that the picture is not interesting. A weaker relationship shows the steps better: violent crime against poverty in the 50 US states.

x_state, y_state = states.poverty.to_numpy(), states.violent.to_numpy()
crime_fit = sm.OLS(y_state, sm.add_constant(x_state)).fit()
slope, slope_se = crime_fit.params[1], crime_fit.bse[1]
shuffled = np.array([np.polyfit(x_state, rng.permutation(y_state), 1)[0] for _ in range(5_000)])


def plot_crime(ax: Axes) -> None:
    """States with the fitted line."""
    ax.scatter(x_state, y_state, s=14, color=style.BLUE, alpha=0.7)
    grid = np.array([8, 19])
    ax.plot(grid, crime_fit.params[0] + slope * grid, color=style.CORAL)
    style.note(ax, f"slope {slope:.1f}\nSE {slope_se:.1f}", x=0.03, y=0.98, color=style.INK)
    ax.set(title="The data", xlabel="poverty (%)", ylabel="violent crimes per 100,000")


def plot_shuffles(ax: Axes) -> None:
    """Slopes from shuffled data, with the normal curve of the standard error."""
    width = style.hist(ax, shuffled, bins=50, color=style.GREY)
    grid = np.linspace(-45, 45, 300)
    ax.plot(grid, stats.norm.pdf(grid, 0, slope_se) * len(shuffled) * width, color=style.CORAL)
    ax.axvline(slope, color=style.INK, ls="--", lw=1.4)
    share = (abs(shuffled) >= abs(slope)).mean()
    style.note(ax, f"shuffles at least as far: {share:.3f}", x=0.03, y=0.98, color=style.INK)
    ax.set(title="If poverty did not matter", xlabel="slope", ylabel="shuffles")


def plot_t_scale(ax: Axes) -> None:
    """The t distribution with n − 2 degrees of freedom and the observed t."""
    t = slope / slope_se
    x = np.linspace(-4.5, 4.5, 400)
    curve = stats.t.pdf(x, len(x_state) - 2)
    ax.plot(x, curve, color=style.BLUE)
    for tail in (x >= abs(t), x <= -abs(t)):
        ax.fill_between(x[tail], curve[tail], color=style.CORAL, alpha=0.5, lw=0)
    ax.axvline(t, color=style.INK, ls="--", lw=1.4)
    style.note(
        ax, f"t = {t:.2f}\np = {crime_fit.pvalues[1]:.3f}", x=0.03, y=0.98, color=style.INK
    )
    ax.set(title="t with 48 df", xlabel="t = slope / SE", yticks=[])


fig, axes = style.subplots(3, height=3.1, width=10.5)
plot_crime(axes[0])
plot_shuffles(axes[1])
plot_t_scale(axes[2])
fig.tight_layout()
plt.show()

Figure 4.1.3: Testing the slope of violent crime on poverty in the 50 states. Left: the data and the fitted line. Middle: the slopes from 5,000 shuffles of the crime numbers among the states (bars), which breaks any relation, with the normal curve of the slope’s standard error (red); the real slope (dashed) is out in the tail. Right: the same on the t scale: t = slope / standard error, with the two tails beyond the observed t shaded. Both ways give p ≈ 0.03.

The formula for the standard error says what makes a slope precise: little noise ($\hat\sigma$ small), many points, and a wide spread of x ($\sum (x_i - \bar x)^2$ large). The last one is easy to forget:

def plot_spread_of_x(ax: Axes, low: float, high: float, title: str) -> None:
    """Fitted lines from many samples with x in [low, high]."""
    grid = np.array([0, 10])
    slopes = []
    for _ in range(50):
        x = rng.uniform(low, high, 20)
        y = 2 + 0.5 * x + rng.normal(0, 1, 20)
        b, a = np.polyfit(x, y, 1)
        slopes.append(b)
        ax.plot(grid, a + b * grid, color=style.BLUE, lw=0.6, alpha=0.4)
    ax.plot(grid, 2 + 0.5 * grid, color=style.INK, ls="--", lw=1.5)
    ax.axvspan(low, high, color=style.GOLD, alpha=0.12, lw=0)
    style.note(ax, f"slopes vary by ±{np.std(slopes):.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set(title=title, xlabel="x", ylim=(-2, 10))


fig, axes = style.subplots(2, height=3, width=9.5, sharey=True)
plot_spread_of_x(axes[0], 4, 6, "x between 4 and 6")
plot_spread_of_x(axes[1], 0, 10, "x between 0 and 10")
axes[0].set_ylabel("y")
fig.tight_layout()
plt.show()

Figure 4.1.4: Lines fitted to 50 simulated samples of 20 points each, with the same true line (dashed) and the same noise. Left: x spread over a narrow range: the slopes swing wildly. Right: x spread over a wide range: the slopes agree. To learn a slope well, measure over a wide range of x.

4.1.4 How well does the line fit?

Without the line, the best guess for any household’s food spending is the average. With the line, it is the point on the line. The R² compares the two: the share of the variation around the average that the line explains.

mean_line = np.full_like(food, food.mean())
fitted = fit.fittedvalues.to_numpy()


def plot_gaps(ax: Axes, to: np.ndarray, title: str, color: str) -> None:
    """Gaps from every household to a reference: the average or the line."""
    ax.vlines(income, to, food, color=color, lw=0.6, alpha=0.6)
    ax.scatter(income, food, s=7, color=style.BLUE, alpha=0.5)
    order = np.argsort(income)
    ax.plot(
        income[order], to[order], color=style.INK, lw=1.4, ls="--" if np.ptp(to) == 0 else "-"
    )
    squares = ((food - to) ** 2).sum() / 1e6
    style.note(ax, f"sum of squares\n{squares:.1f} million", x=0.03, color=style.INK)
    ax.set(
        title=title,
        xlabel="household income",
        ylabel="food spending",
        xlim=(0, 2600),
        ylim=(0, 1600),
    )


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_gaps(axes[0], mean_line, "Around the average", style.GREY)
plot_gaps(axes[1], fitted, f"Around the line: R² = {fit.rsquared:.2f}", style.CORAL)
fig.tight_layout()
plt.show()

Figure 4.1.5: Left: the gaps from each household to the average food spending (dashed): the total variation. Right: the gaps to the line: what is left. R² is one minus the ratio of the squared gaps.

In symbols, with the total sum of squares $\text{TSS} = \sum (y_i - \bar y)^2$ and the residual sum of squares $\text{RSS} = \sum (y_i - \hat y_i)^2$:

$$R^2 = 1 - \frac{\text{RSS}}{\text{TSS}} = \frac{\text{TSS} - \text{RSS}}{\text{TSS}}$$

Here R² = 0.83: income accounts for 83% of the variation in food spending. With a single predictor, R² is simply the square of the correlation, r² (Chapter 3.6). The rest is the noise $\varepsilon$: the gaps that remain. Their typical size is the residual standard error, $\hat\sigma = \sqrt{\text{RSS}/(n-2)}$, 114 here: a household is usually within about ±2 × 114 of the line.

R² depends as much on the noise as on the slope. The same slope can come with any R²:

def plot_noise_level(ax: Axes, noise: float) -> None:
    """Data around the same line with a given amount of noise, with R²."""
    x = rng.uniform(0, 10, 150)
    y = 2 + 0.5 * x + rng.normal(0, noise, 150)
    model = sm.OLS(y, sm.add_constant(x)).fit()
    ax.scatter(x, y, s=8, color=style.BLUE, alpha=0.5)
    ax.plot([0, 10], model.params[0] + model.params[1] * np.array([0, 10]), color=style.CORAL)
    style.note(
        ax,
        f"R² = {model.rsquared:.2f}\nslope {model.params[1]:.2f} ± {model.bse[1]:.2f}",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(title=f"noise sd {noise}", xlabel="x", ylim=(-6, 14))


fig, axes = style.subplots(3, height=2.9, width=10.5, sharey=True)
for ax, noise in zip(axes, (0.5, 1.5, 4)):
    plot_noise_level(ax, noise)
axes[0].set_ylabel("y")
fig.tight_layout()
plt.show()

Figure 4.1.6: Three simulated data sets with the same true line (slope 0.5) and increasing noise. The slope comes out close to 0.5 every time (less precisely as the noise grows), while R² falls from 0.9 to 0.1: R² measures how much of the variation the line explains, not how strong or important the effect is.

R² is a description, not a verdict. A low R² can come with a precise, important slope (many things affect health, and one of them can still matter). A high R² can come from a line that is the wrong shape. And R² never goes down when you add predictors, which makes it a bad guide for choosing them (Chapter 6.4).

4.1.5 Predicting one household, or knowing the average

The same line answers two different questions, with very different uncertainty.

def income_grid() -> np.ndarray:
    """Incomes from 0 to the highest in the data, for drawing lines."""
    return np.linspace(0, income.max(), 100)


def plot_households(ax: Axes, alpha: float) -> None:
    """The 235 households as dots, with the axes they share in both panels."""
    ax.scatter(income, food, s=8, color=style.BLUE, alpha=alpha)
    ax.set(xlabel="household income", ylabel="food spending", xlim=(0, 2500), ylim=(0, 1600))


def mark_interval(ax: Axes, x: float, low: float, guess: float, high: float) -> None:
    """A point with an error bar from low to high."""
    yerr = [[guess - low], [high - guess]]
    ax.errorbar(x, guess, yerr=yerr, fmt="o", ms=4, color=style.CORAL, capsize=5, lw=2)


def plot_prediction(ax: Axes, at: float = 1000) -> None:
    """One household's predicted food spending, inside the band 95% of households fall in."""
    grid = income_grid()
    band = fit.get_prediction(sm.add_constant(grid)).summary_frame(alpha=0.05)
    plot_households(ax, alpha=0.35)
    ax.fill_between(
        grid, band.obs_ci_lower, band.obs_ci_upper, color=style.CORAL, alpha=0.2, lw=0
    )
    ax.plot(grid, band["mean"], color=style.GREY, lw=1.5)
    new = fit.get_prediction(np.array([[1, at]])).summary_frame(alpha=0.05).iloc[0]
    mark_interval(ax, at, new.obs_ci_lower, new["mean"], new.obs_ci_upper)
    ax.set_title("Predict: one household")
    style.note(
        ax, f"income {at:,.0f}: food {new.obs_ci_lower:,.0f} to {new.obs_ci_upper:,.0f}\nwide"
    )


def plot_relationship(ax: Axes, at: float = 1000, resamples: int = 200) -> None:
    """The relationship: lines other samples could have given, and where they cross `at`."""
    grid = income_grid()
    plot_households(ax, alpha=0.12)
    lines = [
        fit_food_on_income(engel.sample(len(engel), replace=True, random_state=rng))
        for _ in range(resamples)
    ]
    for line in lines:
        ax.plot(grid, line.predict(sm.add_constant(grid)), color=style.CORAL, lw=0.6, alpha=0.1)
    averages = [line.params.const + line.params.income * at for line in lines]
    low, high = np.percentile(averages, [2.5, 97.5])
    mark_interval(ax, at, low, np.mean(averages), high)
    ax.set_title("Understand: the relationship")
    style.note(ax, f"average at income {at:,.0f}: food {low:,.0f} to {high:,.0f}\nnarrow")


fig, (left, right) = style.subplots(2, height=3.4)
plot_prediction(left)
plot_relationship(right)
fig.tight_layout()
plt.show()

Figure 4.1.7: One line, two questions, both at an income of 1,000. Left: what will one household spend on food? The red band is where 95% of households fall. Right: what do households spend on average, and how sure can we be from one sample? Each thin line is the line from one bootstrap resample. Two households with incomes above 2,500 are outside the picture.

On the left the question is about one household with an income of 1,000. The line says about 630 on food, but real households with that income spend anywhere from about 410 to 860. This prediction interval includes the noise $\varepsilon$, and more data would not make it much narrower, because households really do differ.

On the right the question is about the average household with that income. Each thin line is the line from one bootstrap resample (Chapter 3.3). They lie close together, so the average is known precisely: between about 620 and 650. This confidence interval shrinks with more data.

A single household is hard to predict while the relationship is known precisely. Both are true at once, because they answer different questions. The two standard errors at an income $x_0$ show why:

$$\text{SE}_{\text{average}}(x_0) = \hat\sigma \sqrt{\frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum (x_i - \bar x)^2}}, \qquad \text{SE}_{\text{one household}}(x_0) = \hat\sigma \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum (x_i - \bar x)^2}}$$

The “1 +” under the second root is the household’s own noise, which no amount of data removes. Both grow with the distance from the average income $\bar x$: the bands are narrowest in the middle of the data and flare out towards the edges, and beyond the data they are not to be trusted at all.

4.1.6 Regression to the mean

The word “regression” comes from a puzzle Francis Galton found in the 1880s: very tall parents tend to have children who are tall, but less tall than themselves. The same happens whenever two measurements are related but not perfectly: those who are extreme on the first tend to be less extreme on the second.

PEOPLE = 2_000
skill = rng.normal(100, 10, PEOPLE)
first, second = skill + rng.normal(0, 10, PEOPLE), skill + rng.normal(0, 10, PEOPLE)
top, bottom = first > np.quantile(first, 0.9), first < np.quantile(first, 0.1)


def plot_two_tests(ax: Axes) -> None:
    """Second score against first, with the diagonal and the regression line."""
    ax.scatter(first, second, s=5, color=style.GREY, alpha=0.4)
    ends = np.array([55, 145])
    ax.plot(ends, ends, color=style.INK, ls="--", lw=1)
    slope, intercept = np.polyfit(first, second, 1)
    ax.plot(ends, intercept + slope * ends, color=style.CORAL)
    ax.set(
        title="Two tests",
        xlabel="first test",
        ylabel="second test",
        xlim=(55, 145),
        ylim=(55, 145),
    )


def plot_groups(ax: Axes) -> None:
    """Average first and second score of the top and bottom 10% on the first test."""
    for group, color, name in (
        (top, style.CORAL, "top 10%"),
        (bottom, style.BLUE, "bottom 10%"),
    ):
        scores = [first[group].mean(), second[group].mean()]
        ax.plot([0, 1], scores, "o-", color=color, ms=7, label=name)
        for x, y in zip((0, 1), scores):
            ax.annotate(
                f"{y:.0f}",
                (x, y),
                xytext=(8, -4),
                textcoords="offset points",
                fontsize=9,
                color=color,
            )
    ax.axhline(100, color=style.INK, ls="--", lw=1)
    ax.set_xticks([0, 1], ["first test", "second test"])
    ax.legend(fontsize=8)
    ax.set(title="The extremes move to the middle", ylabel="average score", xlim=(-0.3, 1.4))


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_two_tests(axes[0])
plot_groups(axes[1])
fig.tight_layout()
plt.show()

Figure 4.1.8: Two tests of the same 2,000 people, each score = skill + luck. Left: the regression line (red) is flatter than the diagonal (dashed). Right: the top 10% on the first test (red) score lower on the second, and the bottom 10% (blue) score higher, without anything having changed.

A top score on the first test is partly skill and partly good luck. The skill stays; the luck does not. So the second score falls back towards the average, even though nobody changed.

This catches people out constantly. Patients enrolled because their blood pressure is high will show lower readings next time, treatment or not. A team praised after its best season does worse the next. Without a comparison group, regression to the mean looks like an effect.

4.1.7 In one picture

• Least squares picks the line with the smallest sum of squared residuals: $\hat b = \sum (x_i - \bar x)(y_i - \bar y) / \sum (x_i - \bar x)^2$.
• The slope’s t is the slope divided by its standard error, $\hat\sigma / \sqrt{\sum (x_i - \bar x)^2}$; its p-value comes from a t distribution with n − 2 degrees of freedom.
• The slope is the change in Y per unit of X; its standard error, p-value and interval say how sure we are.
• R² $= 1 - \text{RSS}/\text{TSS}$ is the share of variation the line explains, and depends on the noise as much as on the slope; the residual standard error is the typical gap that is left.
• A prediction interval (one new case) is much wider than a confidence interval (the average).
• Extreme cases regress to the mean on a second measurement.

Next: several predictors at once.