The Feynman experience · Part 4

Simple linear regression

evergreen · Mar 16, 2020

Introduction

We want to model a linear relationship between a feature $x$ and a scalar target $y$. For a friendly visual intuition before the math, Luis Serrano’s explanation is excellent.¹

The model

We model the target as

where $\epsilon \sim \text{Normal}(0, \sigma^2)$. Because $Y$ is defined through $\epsilon$, it is itself a random variable, with expected value

and variance

Note the variance of $Y$ is the same for every $x$ — a property called homogeneity of variance, or homoscedasticity.²

Let’s generate some data to work with. As throughout this series, we work from scratch in plain Python with NumPy and SciPy:

import numpy as np

np.random.seed(51)

real_beta_0, real_beta_1, real_sigma2 = 2, 3, 1
n = 100

x = np.linspace(0, 5, n)
y = real_beta_0 + real_beta_1 * x + real_sigma2 * np.random.randn(n)

Parameter estimation

Given a sample $\{(x_1, y_1), \dots, (x_n, y_n)\}$, we estimate $\beta_0$ and $\beta_1$ by least squares. The error metric is the residual sum of squares:³

and the estimates minimize it:

Differentiating with respect to each parameter and setting the result to zero — first for $\beta_0$,

and then for $\beta_1$,

where $\bar{x} = \frac{1}{n}\sum_i x_i$ and $\bar{y} = \frac{1}{n}\sum_i y_i$. These two normal equations form a linear system,

whose solution is

Plugging in our data:

beta_1 = np.sum(y * (x - np.mean(x))) / np.sum((x - np.mean(x)) ** 2)
beta_0 = np.mean(y) - beta_1 * np.mean(x)

np.round([beta_0, beta_1], 3)
# → array([1.925, 3.028])

Close to the true $(2, 3)$. We could also reach the same place with a different algorithm — gradient descent,⁴ say — or hand the objective to SciPy’s optimizer:

from scipy.optimize import minimize

def rss(beta, x, y):
    return np.sum((y - beta[0] - beta[1] * x) ** 2)

beta_0, beta_1 = minimize(rss, [0, 0], args=(x, y)).x

np.round([beta_0, beta_1], 3)
# → array([1.925, 3.028])

Properties of the estimators

Since $\hat{\beta}_0$ and $\hat{\beta}_1$ are functions of the random sample, they are themselves random variables, with their own expected values and variances.⁵⁶

Slope $\hat{\beta}_1$

Writing $\hat{\beta}_1 = \sum_i c_i Y_i$ with $c_i = \frac{x_i - \bar{x}}{\sum_j (x_j - \bar{x})^2}$,

Now $\sum_i c_i = 0$ (the numerator $\sum_i (x_i - \bar{x})$ vanishes), and

so $E(\hat{\beta}_1) = \beta_1$ — the estimator is unbiased. Its variance, using $V(Y_i) = \sigma^2$, is

In practice we rarely know $\sigma^2$, so we substitute the estimator $\hat{\sigma}^2 = \frac{\sum_i (y_i - \hat{y})^2}{n - 2}$. Computing the variance for our data:

b1_var = real_sigma2 / np.sum((x - np.mean(x)) ** 2)
np.round(b1_var, 5)
# → 0.0047

To check it, we can estimate $\beta_1$ across a hundred fresh samples and look at the spread:

np.random.seed(51)
y_100 = real_beta_0 + real_beta_1 * x + real_sigma2 * np.random.randn(100, n)

b1 = [minimize(rss, [0, 0], args=(x, yi)).x[1] for yi in y_100]

Intercept $\hat{\beta}_0$

The same approach applies to the intercept. From $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}$,

and, since $\mathrm{Cov}(\bar{Y}, \hat{\beta}_1) = 0$,

b0_var = (1 / n + np.mean(x) ** 2 / np.sum((x - np.mean(x)) ** 2)) * real_sigma2
np.round(b0_var, 5)
# → 0.03941

A confidence band for the mean response

The fitted mean at a point $x_0$, $\hat{\mu}_{Y|x_0} = \hat{\beta}_0 + \hat{\beta}_1 x_0$, is also a random variable. Its expected value is

and, writing it as $\bar{y} + \hat{\beta}_1(x_0 - \bar{x})$ and using $\mathrm{Cov}(\bar{y}, \hat{\beta}_1) = 0$, its variance is

Turning that into a band around the fit:

y_hat = beta_0 + beta_1 * x

y_hat_var = (np.sum((y - y_hat) ** 2) / (n - 2)) * \
    (1 / n + (x - np.mean(x)) ** 2 / np.sum((x - np.mean(x)) ** 2))

y_hat_upper = y_hat + 1.96 * np.sqrt(y_hat_var)
y_hat_lower = y_hat - 1.96 * np.sqrt(y_hat_var)

The band is not parallel to the line — and the variance above tells us why: it grows with $(x_0 - \bar{x})^2$, so our estimate of the mean is most precise near the centre of the data and least precise out at the extremes.

Coefficient of determination

The coefficient of determination, $R^2$, is the share of the target’s variance explained by the feature. It lies between 0 and 1:

r2 = 1 - np.sum((y - y_hat) ** 2) / np.sum((y - np.mean(y)) ** 2)
np.round(r2, 2)
# → 0.95

So the model explains about 95% of the variance in $y$ — unsurprising, since we built the data from a line.

Interpretation

The coefficients read directly: a one-unit change in the predictor shifts the mean target by $\beta_1$, and when the predictor is 0 the mean target is $\beta_0$.

Conclusion

Linear regression is the standard first candidate for a regression problem, and a good choice when you want to understand how each predictor affects the target. The natural next step is more than one feature — see multiple linear regression, the next part of this series.