Baby's First Linear Regression

SCULLY: What kind of business are you in?
WELL-MANICURED MAN: We predict the future…
— The X Files, S3E1

Problem

How can we use data that we already have to predict what data we might see in the future?

For example, it’s common wisdom in political science that we can predict how large chunks of the United States electorate will vote for president many months before anyone casts a ballot using purely economic data. We already have data about the state of the economy before previous presidential elections; is there a way we can use this information to forecast future elections?

I recently read Democracy for Realists by Christopher Achen and Larry Bartels, where they did just that: they built a model that predicts the presidential popular vote margin to surprising accuracy just from how much Americans’ disposable income has grown (or shrunk) and how long the incumbent party has held the presidency. The most surprising thing is that this relationship is essentially a straight line.

They only had election data from 1952-2012 (violet) and produced the line of best fit (orange) that you see here. I’ve added in data for the 2016, 2020, and 2024 presidential elections (blue). For 2016 and 2024 they’re correct within a reasonable margin of error; only the plague year 2020 (a clear outlier) failed to follow the pattern. (Although it did get it catastrophically wrong, it was a catastrophic time.)

The x-axis here is a single variable that represents how much real disposable income the average voter had gained in the final quarters leading up to the election, adjusted by the length (or tenure) of the incumbent party in office. Presumably, as soon as the necessary economic data is available close to the 2028 (or 2032, or 2036, or …) election, this model will be able to provide a relatively solid prediction again.

How do we build something like that ourselves?

The Vocabulary of Data Science

Data Points: Inputs and Outputs

In the case of the Achen and Bartels’s work, we consider each election to be a single point of data living on a 2D Cartesian grid. In this case we only have two dimensions, because we have only one input variable (tenure-adjusted income growth) and one output variable (popular vote percentage). In general, we can have many input variables: for example, if we’re predicting COVID mortality based on patient characteristics, we would probably want input variables like age, height, weight, and so on. Each of these additional pieces of information aren’t new data points, they’re more like enhancements to each point: each data point represents more information. Each of these additional input variables, representing new pieces of information about each data point, are called characteristics (or features): age is a characteristic, height is a characteristic, etc.

This changes the number of dimensions that our data “lives in”: Achen and Bartels’s work only lived in two dimensions because they had one input variable and one output variable. In general, we can have arbitrarily many characteristics, and so arbitrarily many dimensions. We use the variable $d$ to denote the number of input dimensions we’re working in: we only typically target one output variable $\textcolor{#fb7185}{y}$ at a time, and will create one model per output variable.

Dimensions and samples

In general, the more dimensions $d$ we have then the more data we want, because a higher $d$ means we expect a more detailed relationship between the inputs and outputs and a higher $n$ gives us enough information to resolve those details. Trying to capture a detailed picture, like the Mona Lisa, with only a few pixels simply doesn’t give us a clear idea of the picture!

The Mona Lisa downsampled to 15x22 pixels (about 0.2% of its original resolution).

This is why Achen and Bartels, with fewer than 20 elections to work with, compressed the tenure of incumbent parties and income growth into a single variable, rather than keeping them separate. This technique of reducing $d$ without “throwing information away” is called dimensionality reduction, and is something we will probably discuss in another article. For now, the big idea is that more data is always better, doubly so if we have many dimensions.

This kind of spatial perspective on our data is why we use the term “data point”: we’re talking about a single point in space. We write the \(i^{th}\) point as a vector: \(\textcolor{#4ade80}{\vec x^{(i)}} = \{\textcolor{#4ade80}{x_1}, \textcolor{#4ade80}{x_2}, \ldots, \textcolor{#4ade80}{x_d}\}\) where $\textcolor{#4ade80}{x_1}$ is one piece of info, $\textcolor{#4ade80}{x_2}$ is the second piece of info, and so on. (The order doesn’t matter so long as we’re being consistent across data points.) For our purposes, a vector is just a list of numbers that we can do math with.

Hypothesizing a Relationship

We call a variable that we want to predict a target variable, and write it as $\textcolor{#fb7185}{y^{(i)}}$. Given $n$ data points, we denote our data $D = {(\textcolor{#4ade80}{\vec{x}^{(i)}}, \textcolor{#fb7185}{y^{(i)}})}_{i=1}^n$, meaning we have $n$ pairs of data: $\textcolor{#4ade80}{x}$ is an input, $\textcolor{#fb7185}{y}$ is an output. Where the first pair is written $(\textcolor{#4ade80}{\vec{x}^{(1)}}, \textcolor{#fb7185}{y^{(1)}})$, the second as $(\textcolor{#4ade80}{\vec{x}^{(2)}}, \textcolor{#fb7185}{y^{(2)}})$, and so on.

We assume that there is some relationship \(\textcolor{#38bdf8}{f}\) roughly correlating these variables with the election result:

\[\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x}) = \textcolor{#fb7185}{y}\]

We put in our inputs $\textcolor{#4ade80}{\vec x}$, and then we get the data $\textcolor{#fb7185}{y}$ that we want to predict. Our goal is to create a model that accurately represents what this unknown function $\textcolor{#38bdf8}{f}$ does.

Unfortunately, there’s almost always some kind of random interference: the points kind of make a line, but they don’t make an exact line. For Achen and Bartels to make a perfect prediction, they’d practically have to know what every single person in America was thinking on election night and whether they could make it to the polling stations or not. That’s just not feasible. We call this random interference error or noise, and denote it \(\epsilon\). We typically assume that this noise is normally distributed: it tends to clump up around the true value of $\textcolor{#38bdf8}{f}$, is just as likely to overshoot $\textcolor{#38bdf8}{f}$ as to undershoot, its behavior is totally independent of the input variables, and so on. This is all describing noise that is normally distributed: if we were to take a huge number of samples and plot the noise, it would look like

So we can revise our earlier equation: we assume that this data is the result of sampling some function

\[\textcolor{#fb7185}{y^{(i)}} = \textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(i)}}) + \epsilon^{(i)}\]

where $\textcolor{#fb7185}{y^{(i)}}$ is the output, $\epsilon^{(i)}$ the random noise that happened to intervene this time, $\textcolor{#4ade80}{\vec x^{(i)}}$ are the values of all the dependent variables, and $\textcolor{#38bdf8}{f}$ is some unknown function. The superscript $(i)$ denotes that we’re just looking at a single data point: the $i$th time we got a data point. For example, since the 1952 election is the first election that Achen and Bartels’s model considers, it’d have $i = 1$, then 1956 would be $i = 2$, and so on.

We very often assume the relationship \(\textcolor{#38bdf8}{f}\) that we’re studying is linear: it’s essentially the simplest relationship that can exist. Nonlinear relationships tend to have some very complex phenomena involved: feedback loops, compounding effects, and very intricate or sensitive relationships between the different inputs. When we don’t anticipate anything too wild going on, linearity is a reasonable assumption. This is why we would choose linear regression as our model.

Linearity

A linear function $\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x})$ is simply a linear combination of the inputs: that is, we just multiply each of the inputs by a constant and add them all together:

\[\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x}) = \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2} + \cdots + \textcolor{#38bdf8}{\theta_d}\textcolor{#4ade80}{x_d}\]

For example, if I go to the grocery store and buy $A$ apples and $B$ bananas, and apples cost $0.25 and bananas cost $0.50, then the total cost $C$ of my shopping cart will be

\[0.25A + 0.5B = C\]

This kind of relationship can be used to explain things as complicated as the stock market, the efficacy of medical treatment, or whether skyscrapers will fall apart. These relationships are shockingly common.

Notice that we don’t multiply any of the variables together, or raise them to powers, or anything fancy. Now, under our assumption that this is what the true value of $\textcolor{#38bdf8}{f}$ actually is, we just need to find each of the $d$ unknown $\textcolor{#38bdf8}{\theta_i}$ constants. The field of machine learning is essentially about creating algorithms to find these kind of constants, called parameters. We think of these parameters as being elements of an ordered list of numbers, called a vector.

Vectors and the Dot Product

You’ve already seen a vector: $\textcolor{#4ade80}{\vec x^{(i)}}$, the list of a single data point’s characteristics, is our vector of inputs.

The equation we’ve seen above, $0.25A + 0.5B = C$ can actually be expressed as a relationship between two vectors: the vector containing the prices for the fruits, and the vector containing the number of fruits that we purchased:

\[\begin{bmatrix} 0.25 \\ 0.5 \end{bmatrix} \cdot \begin{bmatrix} A \\ B \end{bmatrix} = 0.25A + 0.5B = C\]

This operation, where we multiply the first element from each vector, the second item from each vector, etc. and all them all up is called thedot product. If you’re a programmer, this implementation by Peter Norvig is one of my favorite lines of code ever:

def dot(v, u):
    return sum(a * b for a, b in zip(v, u))

We take two vectors, u and v, and simply sum up the products of each of their entries. Obviously, the vectors must have the same number of elements in order for the dot product to work.

So, our assumption is that a $\textcolor{#38bdf8}{f}$ (the true function, with no noise involved) takes on the form

\[\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(i)}}) = \begin{bmatrix} \textcolor{#38bdf8}{\theta_1} \\ \textcolor{#38bdf8}{\theta_2} \end{bmatrix} \cdot \begin{bmatrix} \textcolor{#4ade80}{x_1^{(i)}} \\ 1 \end{bmatrix} = \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(i)}} + \textcolor{#38bdf8}{\theta_2} 1 = \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(i)}} + \textcolor{#38bdf8}{\theta_2}\]

This is just $\textcolor{#fb7185}{y} = m\textcolor{#4ade80}{x} + b$ but where $m = \textcolor{#38bdf8}{\theta_1}$ and $b = \textcolor{#38bdf8}{\theta_2}$, and again we use the $(i)$ superscripts to indicate that this function is just taking the $i$-th data point to the $i$-th output. Notice that we also cleverly included that extra $1$ on the input column to account for the constant $\textcolor{#fb7185}{y}$-intercept term, $b$.

If we have many input features $\textcolor{#4ade80}{x_1}, \textcolor{#4ade80}{x_2}, \ldots, \textcolor{#4ade80}{x_d}$, then the general form for a linear function is:

\[\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(i)}}) = \underbrace{\begin{bmatrix} \textcolor{#38bdf8}{\theta_1} \\ \textcolor{#38bdf8}{\theta_2} \\ \vdots \\ \textcolor{#38bdf8}{\theta_d} \\ \textcolor{#38bdf8}{\theta_{d+1}} \end{bmatrix}}_{\textcolor{#38bdf8}{\Large \vec \theta}} \cdot \underbrace{\begin{bmatrix} \textcolor{#4ade80}{x_1^{(i)}} \\ \textcolor{#4ade80}{x_2^{(i)}} \\ \vdots \\ \textcolor{#4ade80}{x_d^{(i)}} \\ 1 \end{bmatrix}}_{\textcolor{#4ade80}{\Large \vec x^{(i)}}} = \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(i)}} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2^{(i)}} + \cdots + \textcolor{#38bdf8}{\theta_d}\textcolor{#4ade80}{x_d^{(i)}} + \textcolor{#38bdf8}{\theta_{d+1}}\]

The goal of linear regression is to guess the values in the vector $\textcolor{#38bdf8}{\mathbf{\vec \theta}}$. With some clever tricks and enough samples, this would be easy if not for the noise, $\epsilon$, which gets added in.

Before we move on, I have one final technical note about the dot product: notice that if we take the dot product of a vector with itself and then take the square root, we get the $n$-dimensional (Euclidean) distance formula for the distance between the point in space $(x_1, x_2, \ldots, x_n)$ and the origin $(0, 0, \ldots, 0)$ (which we denote $\vec 0$):

\[\sqrt{\vec x \cdot \vec x} = \sqrt{x_1^2 + x_2^2 + \cdots + x_n^2} = \sqrt{(x_1 - 0)^2 + (x_2 - 0)^2 + \cdots + (x_n - 0)^2} = \text{dist}(\vec 0, \vec x)\]

This will be very useful to us as we figure out how we can actually get a good estimate for $\textcolor{#38bdf8}{\mathbf{\vec \theta}}$.

\(\mathbf{\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x}) + \epsilon = \textcolor{#fb7185}{\vec y}}\) as a Matrix Equation

The equation

\[\textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x}) = \underbrace{\begin{bmatrix} \textcolor{#38bdf8}{\theta_1} \\ \textcolor{#38bdf8}{\theta_2} \\ \vdots \\ \textcolor{#38bdf8}{\theta_d} \\ \textcolor{#38bdf8}{\theta_{d+1}} \end{bmatrix}}_{\textcolor{#38bdf8}{\Large \vec \theta}} \cdot \underbrace{\begin{bmatrix} \textcolor{#4ade80}{x_1} \\ \textcolor{#4ade80}{x_2} \\ \vdots \\ \textcolor{#4ade80}{x_d} \\ 1 \end{bmatrix}}_{\textcolor{#4ade80}{\Large \vec x}} = \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2} + \cdots + \textcolor{#38bdf8}{\theta_d}\textcolor{#4ade80}{x_d} + \textcolor{#38bdf8}{\theta_{d+1}}\]

is useful for describing the relationship we want to learn about, but the problem is we don’t have a similar way of representing the data we actually have in a similar form.

This is where we can introduce a major strength of linear algebra: its ability to represent a whole bunch of inputs being matched to their outputs at once. This is called a system of equations. For example, if we have $n = 4$ data points and $d = 2$ features:

\[\begin{aligned} \textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(1)}}) + \epsilon^{(1)} &= \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(1)}} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2^{(1)}} + \textcolor{#38bdf8}{\theta_3} + \epsilon^{(1)} = \textcolor{#fb7185}{y^{(1)}}\\ \textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(2)}}) + \epsilon^{(2)} &= \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(2)}} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2^{(2)}} + \textcolor{#38bdf8}{\theta_3} + \epsilon^{(2)} = \textcolor{#fb7185}{y^{(2)}}\\ \textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(3)}}) + \epsilon^{(3)} &= \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(3)}} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2^{(3)}} + \textcolor{#38bdf8}{\theta_3} + \epsilon^{(3)} = \textcolor{#fb7185}{y^{(3)}}\\ \textcolor{#38bdf8}{f}(\textcolor{#4ade80}{\vec x^{(4)}}) + \epsilon^{(4)} &= \textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{x_1^{(4)}} + \textcolor{#38bdf8}{\theta_2}\textcolor{#4ade80}{x_2^{(4)}} + \textcolor{#38bdf8}{\theta_3} + \epsilon^{(4)} = \textcolor{#fb7185}{y^{(4)}} \end{aligned}\]

We want to represent all of our observed inputs and outputs in a single formula, and somehow go from that formula for our observations to solve for $\textcolor{#38bdf8}{\vec \theta}$. If you imagine packing each $\textcolor{#fb7185}{y^{(i)}}$ into a vector, like

\[\textcolor{#fb7185}{\vec y} = \textcolor{#fb7185}{\begin{bmatrix} y^{(1)} \\ y^{(2)} \\ y^{(3)} \\ y^{(4)} \end{bmatrix}}\]

and then doing the same for each observed value of $\textcolor{#4ade80}{x_1}, \textcolor{#4ade80}{x_2}$, and $\epsilon$:

\[\textcolor{#38bdf8}{\theta_1}\textcolor{#4ade80}{\underbrace{\begin{bmatrix} x_1^{(1)} \\ x_1^{(2)} \\ x_1^{(3)} \\ x_1^{(4)} \end{bmatrix}}_{\Large \vec x_1}} + \textcolor{#38bdf8}{\theta_2}\underbrace{\textcolor{#4ade80}{\begin{bmatrix} x_2^{(1)} \\ x_2^{(2)} \\ x_2^{(3)} \\ x_2^{(4)} \end{bmatrix}}}_{\textcolor{#4ade80}{\Large \vec x_2}} + \textcolor{#38bdf8}{\theta_3}\underbrace{\begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}}_{\Large \vec 1} + 1\underbrace{\begin{bmatrix} \epsilon^{(1)} \\ \epsilon^{(2)} \\ \epsilon^{(3)} \\ \epsilon^{(4)} \end{bmatrix}}_{\Large \vec \epsilon} = \textcolor{#fb7185}{\underbrace{\begin{bmatrix} y^{(1)} \\ y^{(2)} \\ y^{(3)} \\ y^{(4)} \end{bmatrix}}_{\Large \vec y}}\]

Then our system of equations really looks like we again have a dot product, but where one vector is a regular vector of numbers and the other is a vector of vectors. This notation makes sense if we define two operations, known as scalar multiplication and vector addition:

\[\underbrace{\textcolor{#f5a524}{s}\textcolor{#94a3b8}{\begin{bmatrix} v_1 \\ \vdots \\ v_d \end{bmatrix}} = \textcolor{#94a3b8}{\begin{bmatrix} \textcolor{#f5a524}{s}v_1 \\ \vdots \\ \textcolor{#f5a524}{s}v_d \end{bmatrix}}}_{\text{scalar multiplication}} \quad \quad \underbrace{\textcolor{#f5a524}{\begin{bmatrix} v_1 \\ \vdots \\ v_d \end{bmatrix}} + \textcolor{#e879f9}{\begin{bmatrix} u_1 \\ \vdots \\ u_d \end{bmatrix}} = \begin{bmatrix} \textcolor{#f5a524}{v_1} + \textcolor{#e879f9}{u_1} \\ \vdots \\ \textcolor{#f5a524}{v_d} + \textcolor{#e879f9}{u_d} \end{bmatrix}}_{\text{vector addition}}\]

So we once again express this entire formula in terms of a dot product, but this time between a “regular” vector of numbers $\textcolor{#38bdf8}{\vec \theta}$ and some weird “vector of vectors” that we don’t have any term for.

Well, that vector of vectors is called a matrix; rather than a list of numbers, a matrix acts more like a table of numbers: the vectors \(\textcolor{#4ade80}{\vec x_1}, \textcolor{#4ade80}{\vec x_2}, \textcolor{#4ade80}{\vec x_3}\) are called the columns of the matrix, because they form columns exactly the same the entries in a spreadsheet do. This behavior where we take the dot product of the vector with the columns of the matrix are called matrix multiplication. In matrix notation, this equation then becomes

\[\textcolor{#4ade80}{\underbrace{\begin{bmatrix} x_1^{(1)} & x_2^{(1)} & \textcolor{var(--text-color)}{1} \\ x_1^{(2)} & x_2^{(2)} & \textcolor{var(--text-color)}{1} \\ x_1^{(3)} & x_2^{(3)} & \textcolor{var(--text-color)}{1} \\ x_1^{(4)} & x_2^{(4)} & \textcolor{var(--text-color)}{1} \\ \end{bmatrix}}_{\Large X}} \textcolor{#38bdf8}{\underbrace{\begin{bmatrix} \theta_1 \\ \theta_2 \\ \theta_3 \end{bmatrix}}_{\Large \vec \theta}} + \underbrace{\begin{bmatrix} \epsilon^{(1)} \\ \epsilon^{(2)} \\ \epsilon^{(3)} \\ \epsilon^{(4)} \end{bmatrix}}_{\Large \vec \epsilon} = \textcolor{#fb7185}{\underbrace{\begin{bmatrix} y^{(1)} \\ y^{(2)} \\ y^{(3)} \\ y^{(4)} \end{bmatrix}}_{\Large \vec y}}\]

We call the vector of observations (with an extra column of 1’s tacked on) $\textcolor{#4ade80}{X}$. This whole thing is written succinctly as

\[\textcolor{#4ade80}{X}\textcolor{#38bdf8}{\vec\theta} + \vec \epsilon = \textcolor{#fb7185}{\vec y}\]

Summary: The Problem Statement

So, with all of this matrix notation established, we can write our formula:

\[\textcolor{#4ade80}{X}\textcolor{#38bdf8}{\vec\theta} + \vec\epsilon = \textcolor{#fb7185}{\vec y}\]

where $\textcolor{#4ade80}{X}$ is the matrix of our observed data points (with an extra column of 1’s on the end to allow for a constant term), $\textcolor{#fb7185}{\vec y}$ are the observed target variables, $\vec\epsilon$ is the “error” or “noise” that causes the data to not be perfectly linear, and $\textcolor{#38bdf8}{\vec\theta}$ is the vector that provides the linear relationship between the independent variables and $\textcolor{#fb7185}{\vec y}$.

Solving for the Line of Best Fit

So, the noise prevents us from solving the equation directly. What we can do instead is define a cost function, often called a loss function, that takes in a vector and outputs a positive real number representing how closely the formula for the line associated with $\textcolor{#38bdf8}{\vec \theta}$ matches the data. Smaller output corresponds to less “error” (indicating a better fit), and a higher output corresponds to a worse one. We often denote this function $J(\textcolor{#38bdf8}{\vec\theta})$.

So, our goal is to find the input (or “argument”) $\textcolor{#38bdf8}{\vec\theta^*}$ that yields the smallest possible value $J(\textcolor{#38bdf8}{\vec\theta^*})$. We call this desired value the argument minimum or argmin, formally written $\text{argmin}_{\textcolor{#38bdf8}{\vec\theta}} J(\textcolor{#38bdf8}{\vec \theta})$, and it corresponds to the equation for the line that best fits our data.

Now, our attention simply turns to finding a loss function that

1. Produces a line of best fit that gives us useful information, and
2. Isn’t too difficult to compute.

The standard and often most logical choice for a loss function is called residual sum of squares, denoted $\RSS(\textcolor{#38bdf8}{\vec \theta})$. The method of using $\RSS$ to find the line of best fit is called ordinary least squares.

Using Ordinary Least Squares

Residual sum of squares is basically just like the Euclidean distance function: the point is that we try to find a line that minimizes the vertical distances, called residuals, between a proposed line of best fit (which we denote $h_{\textcolor{#38bdf8}{\theta}}$) and the data points. This basically just means that for any input $x^{(i)}$, we try to make sure the total of all the Euclidean distance between our prediction \(h_{\textcolor{#38bdf8}{\theta}}(\textcolor{#4ade80}{x^{(i)}})\) and the associated observation $y^{(i)}$ is as small as possible.

Let’s set this all up using the language of matrices.

Deriving the Matrix Form for Residual Sum of Squares

So, first, a formula that measures the sum of all the Euclidean distance between our prediction \(h_{\textcolor{#38bdf8}{\theta}}(\textcolor{#4ade80}{x^{(i)}})\) and the associated observation $\textcolor{#fb7185}{y^{(i)}}$ would look something like

\[\sqrt{\displaystyle{\sum_{i = 1}^n} \, (h_{\textcolor{#38bdf8}{\theta}}( \textcolor{#4ade80}{x^{(i)}}) - \textcolor{#fb7185}{y^{(i)}})^2}\]

But taking the square root here actually doesn’t do anything: a large sum will correlate with a larger square root, and a smaller sum correlates with a smaller square root; minimizing the inside sum minimizes the square root of the sum, so there’s no reason to do the extra work of taking the square root. Instead, we just use

\[J(\textcolor{#38bdf8}{\theta}) = \displaystyle{\sum_{i = 1}^n} \, (h_{\textcolor{#38bdf8}{\theta}}( \textcolor{#4ade80}{x^{(i)}}) - \textcolor{#fb7185}{y^{(i)}})^2\]

Where $h_{\textcolor{#38bdf8}{\theta}}(\textcolor{#4ade80}{\vec x})$ is the hypothesized function for some value of $\textcolor{#38bdf8}{\theta}$. Some algebra reveals that we can actually write this whole thing as a matrix:

\[\def\arraystretch{2.2} \begin{array}{rcl} \RSS(\textcolor{#38bdf8}{\vec \theta}) & = & \displaystyle{\sum_{i = 1}^n} \, (h_{\textcolor{#38bdf8}{\theta}}( \textcolor{#4ade80}{x^{(i)}}) - \textcolor{#fb7185}{y^{(i)}})^2 \\ & = & ||\textcolor{#4ade80}{X} \textcolor{#38bdf8}{\vec \theta} - \textcolor{#fb7185}{\vec y} ||^2 \\ & = & (\textcolor{#4ade80}{X} \textcolor{#38bdf8}{\vec \theta} - \textcolor{#fb7185}{\vec y})^T (\textcolor{#4ade80}{X} \textcolor{#38bdf8}{\vec \theta} - \textcolor{#fb7185}{\vec y}) \\ & = & \underbrace{\textcolor{#38bdf8}{\vec{\theta}}^{\, T}\textcolor{#4ade80}{X}^T\textcolor{#4ade80}{X}\textcolor{#38bdf8}{\vec{\theta}}}_{\text{quadratic}} - \underbrace{2 \textcolor{#fb7185}{\vec{y}}^{\,T} \textcolor{#4ade80}{X} \textcolor{#38bdf8}{\vec \theta}}_{\text{linear}} + \underbrace{\textcolor{#fb7185}{\vec y}^{\, T} \textcolor{#fb7185}{\vec y}}_{\text{constant}} \\ \end{array}\]

This is the surface that corresponds to all the possible values that residual sum of squares can take on for any value of $\textcolor{#38bdf8}{\vec \theta}$, and because it’s quadratic, it forms a nice bowl shape just like the graph $y = x^2$:

The surface forms a smooth bowl, with no valleys or hills. Importantly, there’s a single point at the bottom, which corresponds to the values of $\textcolor{#38bdf8}{\vec \theta}$ that produce the least error (and therefore best fit). This kind of property is called convexity, and it makes our work extremely convenient.

From here, we do the same thing we do in a calculus class to find the minimums of other functions: we take the derivative, set it equal to zero, and solve for $\textcolor{#38bdf8}{\vec \theta}$; the precise solution we want is $\nabla \RSS(\textcolor{#38bdf8}{\vec \theta}) = 0$. The calculations are pretty straightforward if you know some matrix calculus, so I’m going to spare you the details.

After these calculations are done, we find out that the critical point that minimizes $\RSS$ is

\[\textcolor{#38bdf8}{\vec \theta} = (\textcolor{#4ade80}{X}^T\textcolor{#4ade80}{X})^{-1} \textcolor{#4ade80}{X}^T \textcolor{#fb7185}{\vec y}\]

If you remember the formula $x = \dfrac{-b}{2a}$ from algebra for finding the minimum (or maximum) point on a quadratic polynomial $y = ax^2 + bx + c$, this is basically the same thing for the matrix equation for $\RSS(\textcolor{#38bdf8}{\vec \theta})$. These values for $\textcolor{#38bdf8}{\vec \theta}$ are the missing coefficients for our line of best fit!

A big idea here is that finding our line of best fit isn’t an iterative process like deep learning is: for linear regression, we have a formula to compute the optimal solution directly. For computers, this is virtually instantaneous.

Demos: Linear Regression In Action

With all that math and the major concepts out of the way, let’s actually watch ordinary least squares in action rather than just trusting the algebra. The demo below draws a random sample around a true line and fits one back through it with exactly the formula above. The dashed blue line is the true function that we sample points around, the solid orange line is our fitted line of best fit, and the thin verticals are the residuals whose squares get summed up in $\RSS(\textcolor{#38bdf8}{\vec\theta})$.

Noise (σ): 2.0 Sample size (n): 40

True line: · Fitted line: · R² =

Experimenting with the value of $n$ brings us to another really big idea that spans all of data science: our estimate gets better as we get more points. Check out how much closer the line of best fit gets to the hidden true function in this simulation:

Sample & fit

n = 1

Distance to true θ vs. n

Notice how as $n$ gets bigger, the Euclidean difference with the true function gets way smaller.

Evaluating the Model

It’s possible that we might be mistaken in assuming a linear relationship between $\textcolor{#4ade80}{X}$ and $\textcolor{#fb7185}{y}$. It’s also possible that the noise $\epsilon$ is so great that the line doesn’t fit the data very well until we have an enormous number of samples $n$. Even though the linear regression algorithm gives us the fit with the smallest error, we don’t actually know how “small” the smallest error is! We don’t have a way to compare the fits of linear regression between two data sets from the numbers alone, but the demo above already gives a direct visual way to compare a tight fit against a loose one.

But if there are too many features (input variables) for us to create a human-readable graph, we can’t rely on visual intuition. We might also wonder how a scientist knows if a colleague’s linear regression has turned up an accurate, highly predictive relationship, or if it’s just noise. After all, it’s possible to grab any data set, choose some variables to be $\textcolor{#4ade80}{X}$ and one to be $\textcolor{#fb7185}{y}$ completely at random, and fit a linear regression to it.

Fitted line: · R² =

In the example above, there’s no correlation at all between the $\textcolor{#4ade80}{x}$ and $\textcolor{#fb7185}{y}$ variables; they’re completely independent. Hit the generate button a couple times and notice how the line ends up in wildly different places despite the underlying process being exactly the same each time. However, our linear regression algorithm doesn’t prevent us from pretending that there is and trying to fit on it, but we don’t know of any objective way to argue that linear regression is a bad fit here.

What we’re looking for here is a way to perform regression validation: the process of determining how accurate a linear regression is, and how well it describes the data it’s meant to fit. We can develop several different ways to differentiate the good fits from the bad.

Residual Analysis

A natural first place to look would be the residuals — the same thin vertical lines the demo above draws from each point down to the fitted line. We would hope that the distribution of residuals is heavily skewed towards zero, and might look like a Gaussian (if we’re treating distance below the line as negative) or a half-Gaussian (if we’re just looking at absolute distance). We can see exactly this relationship below:

Sample & fit (n = 100)

Residuals

Noise (σ): 2.0

Residual mean: · Residual std. dev.:

At low noise the bars huddle tightly around zero and roughly trace the dashed curve. The histogram spreads as noise increases, but because the underlying model is truly linear it keeps the same rough bell shape and just ends up being wider. The shape itself tells us that the model’s assumptions are valid.

If the distribution has several spikes, we might wonder whether our data set is a mixture of several distributions. This is a typical sign of Simpson’s paradox, where several distributions all individually exhibit one kind of trend or property which disappears or changes when looking at the distributions together.

Two constructed groups (violet, teal), each with its own downward trend (dashed) — but averaged together, the combined trend (solid) points the other way entirely. Not real data; built to make the effect obvious.

If there are very few residuals close to zero, this might indicate that your data set has “holes” in it (and you will likely need some heavier machinery if you want to say more) and that you will want to search for a way to split up the data set a bit depending on where each data point falls. (The search for ways to do this is called data mining, but that’s a huge topic for another time.)

The benefits of residual analysis are that it can very quickly give you an intuitive idea of how well your model is fitting and explaining the data set. The drawbacks are that it’s not very “scientific”; we’re just eyeballing here. If you want some hard numbers, you can calculate differences between the distribution of the residuals and the expected (Gaussian) distribution, and that will give you some information on how good your fit is, but this is entering the realm of descriptive statistics.

Descriptive Statistic: the Coefficient of Determination $R^2$

There are several descriptive statistics that we can calculate for a linear regression; the most common is the coefficient of determination, denoted $R^2$ (sometimes $r^2$), which measures how much of the variance our model explains. In other words, it measures roughly how close all the points in the data set are to our linear regression. If they are very close to the line with only small deviations, the $R^2$ will be very close to 1. If the points are scattered haphazardly with seemingly little relationship to the line, then the $R^2$ will be close to 0. That’s exactly the number the demo above prints out next to the fitted line as you drag the noise slider — watch it fall as the scatter gets messier.

The $R^2$ isn’t the only statistic you should look at, but it’s beyond the scope of this article to describe them all in meaningful depth, so I’ll just use the $R^2$ as a case study. The point is that this is usually the first statistic that we look at when evaluating a model. The benefits of looking at the $R^2$, as with other descriptive statistics, are that they have unambiguous and meaningful statistical interpretations. A low $R^2$ is an almost sure sign of a poor fit — a high one is definitely more reassuring, but not always a sure sign of a good fit. Let’s see how it can get complicated.

Out-of-Sample Validation and Overfitting

Go back to the Achen-Bartels model from the very beginning of this article. Fit on the 1952-2012 elections, it posts an in-sample $R^2$ of about 0.766. By the rule of thumb above that’s a pretty strong fit, but $R^2$ alone can’t tell you whether a model will keep working on data that we haven’t seen yet (called out of sample data). After all, the in-sample $R^2$ is intrinsically biased towards being better because the fitting procedure that produced $\textcolor{#38bdf8}{\vec\theta}$ is the very same procedure that minimizes the error that $R^2$ tracks. It’s the best case that particular sample of noise ${\epsilon^{(1)}, \epsilon^{(2)}, \ldots, \epsilon^{(n)}}$ can possibly produce, not a promise about data the model hasn’t seen.

This is a general problem with models in data science and machine learning: the more flexible a model is, the more capable it is of just memorizing the specific noise in that one data set rather than learning the real underlying pattern. We call this very common pitfall overfitting. A great in-sample $R^2$ is consistent with finding a real, generalizable relationship, but it’s equally consistent with having built an elaborate noise-fitting machine, and the number alone can’t tell you which you have.

It’s a bit tricky to overfit with a linear model (especially on large $n$), so I’m going to demonstrate this with a nonlinear regression that has more flexibility to fit noise. (Fitting nonlinear models, as well as a rigorous definition of model “flexibility”, is something we’ll discuss in a future article.)

Above, both models fit on the same 50 training points (violet). The flexible model wraps through every single one and achieves a perfect $R^2 = 1.00$; the plain line manages only $R^2 = 0.21$ against the same. Dashed line is the true underlying relationship, invisible to both models during fitting. As you might anticipate from the blue line:

…the flexible model’s $R^2$ collapses to an abysmal $-6.93$, worse than just guessing the test set’s average. The plain line’s $R^2$ barely moves to $R^2 = 0.28$, close to its training performance, because it was never flexible enough to memorize noise in the first place.

The only way to actually tell the difference is to test the model on data it never saw while fitting. That’s exactly what I’ve done with the Achen-Bartels model: 2016, 2020, and 2024 are genuine out-of-sample data for a model made after the 2012 election. The model’s real performance there was measurably worse than its 0.766 in-sample $R^2$ alone would suggest: 2016 and 2024 missed by several points, not quite as tight as “explains 77% of the variance” might lead you to expect, and the prediction for 2020 was simply atrocious. (Again, debatable if 2020 can be held against them — but the record will settle with time, and more data.) The final $R^2$ was $-93.90$ (or $-13.67$ omitting 2020)

None of these flaws show up if we look at in-sample data, but the point of creating these models in the first place is to look at data outside our initial sample. “How well does this model explain data?” is not a well-posed question: the point is that in-sample and out of sample data are fundamentally different. Only out of sample data tells us whether the model is trustworthy, or whether it fails to see the forest for the trees.

Practitioners usually formalize this by splitting their data up front, in a practice called a train-test split, into a training set that the model learns from and a (usually smaller) test set, which is used to determine whether the model plays well with data it hasn’t seen before.

Descriptive Statistic: Pearson correlation coefficient

I’ll end here with talking about the Pearson correlation coefficient because we can actually calculate it before we even fit a linear regression at all. In other words, this is a statistic about the data, rather than about the regression. The Pearson correlation coefficient (PCC), usually written $r$, is defined as follows:

\[r = \dfrac{\sum_i \, (\textcolor{#4ade80}{x_i} - \bar{\textcolor{#4ade80}{x}})(\textcolor{#fb7185}{y_i} - \bar{\textcolor{#fb7185}{y}})}{\sqrt{\sum_i \, (\textcolor{#4ade80}{x_i} - \bar{\textcolor{#4ade80}{x}})^2}\sqrt{\sum_i \, (\textcolor{#fb7185}{y_i} - \bar{\textcolor{#fb7185}{y}})^2}}\]

The PCC ranges from 1 to -1. A PCC close to 1 means that the data tightly fits a line with positive slope, and a PCC close to -1 means the data fits a line with a negative slope. A PCC close to 0 means the data doesn’t fit a line. For simple linear regression (a single input variable, exactly the case the demo above is fitting) $r$ and $R^2$ aren’t just related, they’re the same statistic: $R^2 = r^2$. Every time you watched the $R^2$ readout above fall as the scatter got noisier, you were watching $r$ fall too — just squared, which is also why $R^2$ can never tell you the direction of the relationship the way $r$ can.

r=1

r=0.8

r=0.4

r=0

r=-0.4

r=-0.8

r=-1

Sixty points per panel, generated to hit each labeled r exactly (to six decimal places) rather than approximately. Notice the point cloud tightens from a diffuse blob toward a perfect line as the correlation strengthens toward ±1, and that the tilt (not just the tightness) flips with the sign.

There are two edge cases, where the line is either exactly horizontal or exactly vertical; the lines in these cases have neither a positive nor negative slope, which causes the PCC to be undefined. Fortunately, these situations are both extremely rare in practice and immediately obvious to spot.

Data scientists usually perform exploratory data analysis before deciding what model to fit. The PCC is one statistic you’d likely look at before deciding whether to use a linear regression or not in the first place.

The Problems with Descriptive Statistics

Descriptive statistics like the $R^2$ or PCC can be “fooled” or fail to provide a complete picture of the data set. After all, boiling down hundreds or millions of data points into just a small handful of numbers is usually going to fail to capture at least some information.

For example, Anscombe’s quartet is a collection of four datasets that have the same descriptive statistics (including $R^2$ values) despite obviously being wildly different:

I

II

III

IV

Anscombe’s quartet — F.J. Anscombe’s actual 1973 data. All four datasets share the same fitted line ($\textcolor{#fb7185}{y} \approx 0.5\textcolor{#4ade80}{x} + 3$) and the same $R^2$ ($\approx 0.67$).

The data sets on the left are clearly fit well by a linear regression, with the exception of a single outlier in the lower-left set that skews the regression line. The upper right set clearly exhibits a quadratic rather than linear relationship, and the lower right data set is thrown off by a specific kind of outlier called a high leverage point.

More data sets, known as the Datasaurus dozen, have been created to give more examples of data sets with identical descriptive statistics but very different fits — thirteen data sets (including the original Anscombe-inspired “dino”) that all share the same mean, standard deviation, and Pearson correlation to two decimal places, generated by Justin Matejka and George Fitzmaurice via simulated annealing rather than by hand:

Dino

Away

H Lines

V Lines

X Shape

Star

High Lines

Dots

Circle

Bullseye

Slant Up

Slant Down

Wide Lines

Matejka & Fitzmaurice (2017)’s Datasaurus Dozen. Every one of these thirteen datasets shares essentially the same averages, standard deviations, and Pearson correlation coefficients.

The point is that none of these statistics can replace a grounded understanding of what is actually generating these data sets: we started out by assuming that our underlying process was linear and that the noise was normally distributed. If these assumptions are broken, then regular old linear regression is simply not the right tool for the job.

Conclusion

Now that we have reached the end of this article, you are officially licensed to use

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X, y)
model.predict(new_data)

guilt-free. Congratulations.

But we made several assumptions in this article: we assumed that $\textcolor{#38bdf8}{f}$ was linear, we assumed that our noise $\epsilon$ was normally distributed with mean 0 and constant variance, and that all of our input variables were independent. Each of these assumptions is likely to be true for most of the use cases you encounter in real life, but it’s also very likely that at least one of them won’t be true.

Fortunately, the techniques of linear regression generalize beautifully to very exotic processes — even cases where $\textcolor{#38bdf8}{f}$ isn’t linear at all! Revisiting and playing around with our assumptions allows us to mix and match techniques to construct new kinds of models that can predict much wilder and more varied phenomena. Stay tuned for Baby’s Second Linear Regression.