44 Models and Loss Functions

In semester 1, we estimated means and differences using models we took for granted—horizontal lines, all-functions models. Now we’re going to think more carefully about the choices involved: what model to use and why, and what it means to choose the “best” curve from that model.

Error Criteria: Why Least Squares?

Fitting Curves to Data

I’ve talked a lot about fitting curves to data. Fitting means choosing from some set of options. Like sizing a shoe. You’ve got some choices and you decide which is best.

Vague Concept \[ \hat \mu \quad \text{ minimizes } \quad \text{badness}(m) \quad \text{ among the curves } m \text{ we're choosing from.} \]

Precise Realization \[ \hat \mu = \operatorname{argmin}_{m \in \mathcal{M}} \frac{1}{n}\sum_{i=1}^n \{ m(X_i) - Y_i \}^2. \]

We call the set of curves $\mathcal{M}$ we’re choosing from a regression model. People use all kinds of models. It can be any set of curves at all. But for the most part, $\text{badness}$ is always mean squared error.

Apparently Reasonable Options

Mean Absolute Error

\[ \hat \mu = \operatorname{argmin}_{m \in \mathcal{M}} \frac{1}{n}\sum_{i=1}^n | m(X_i) - Y_i | \]

With this criterion, we’re acting like we’d just as soon have

One point where we fit well and one where we fit badly.
Two points where we fit ok.

Maximum Absolute Error

\[ \hat \mu = \operatorname{argmin}_{m \in \mathcal{M}} \max_{i \le n} | m(X_i) - Y_i | \]

With this one, we act like we’d much rather have two points where we fit ok.

Mean Squared Error

\[ \hat \mu = \operatorname{argmin}_{m \in \mathcal{M}} \frac{1}{n}\sum_{i=1}^n \{ m(X_i) - Y_i \}^2 \]

With this one, we express a slight preference for two points where we fit ok.

Why Mean Squared Error?

People say we use mean squared error because it’s convenient. And that’s true. It isn’t because we always prefer to fit two points ok. But that sells it short. It’s convenient in a profound way. To get a sense of why, we’ll start by fitting a constant—a horizontal line.

Fitting Constants

Least Squares

\[ \hat \mu = \operatorname{argmin}_{m \in \mathbb{R}} \frac{1}{n}\sum_{i=1}^n \{ m - Y_i \}^2 \]

It’s the sample mean. \[ \begin{aligned} m=\frac{1}{n}\sum_{i=1}^n Y_i \quad \text{ solves }\quad 0&=\frac{d}{dm} \frac{1}{n}\sum_{i=1}^n \{ m - Y_i \}^2 \\ &=\frac{1}{n}\sum_{i=1}^n 2\{m - Y_i \} = 2\left\{m - \frac{1}{n}\sum_{i=1}^n Y_i\right\}. \end{aligned} \]

Least Absolute Deviations

\[ \hat \mu = \operatorname{argmin}_{m \in \mathbb{R}} \frac{1}{n}\sum_{i=1}^n | m - Y_i | \]

It’s the sample median. \[ \begin{aligned} m \quad \text{with} \quad \#\{Y_i > m\} = \#\{Y_i < m\} \quad \text{ solves }\quad 0&=\frac{d}{dm} \frac{1}{n}\sum_{i=1}^n | m - Y_i | \\ &=\frac{1}{n}\sum_{i=1}^n \operatorname{sign}(m - Y_i) \\ &=\frac{1}{n}\left(\sum_{i : m > Y_i} 1 - \sum_{i:m < Y_i} 1\right). \end{aligned} \]

Least Maximal Absolute Deviation

\[ \hat \mu = \operatorname{argmin}_{m \in \mathbb{R}} \max_{i \le n} | m - Y_i | \]

It’s the center of the range of $Y_i$. \[ \frac{\max_{i} Y_i + \min_{i} Y_i}{2}. \]

The Key Property of the Mean: Decomposability

The mean, median, and center of the range are all decent summaries. If I had to describe my data with one number, they’d be the candidates. But the mean has a property the others don’t. It’s decomposable.

Suppose that, instead of fitting one horizontal line, we fit two.

The left one is the mean of observations on the left. Same for the right.
We can get the mean of all observations by averaging these two means.
An average with weights proportional to number of observations.

\[ \bar Y = \frac{\text{count}_{\text{left}} \bar Y_{\text{left}} + \text{count}_{\text{right}} \bar Y_{\text{right}}}{\text{count}_{\text{left}} + \text{count}_{\text{right}}}. \]

This means we can also predict the average we would get if we changed the number of observations on the left and the right. e.g., if we changed the distribution of class sizes by imposing a cap. We’d just take an average with different weights. \[ \bar Y(\text{cap}) = \frac{\text{count}_{\text{left}}(\text{cap})\ \bar Y_{\text{left}} + \text{count}_{\text{right}}(\text{cap}) \ \bar Y_{\text{right}}}{\text{count}_{\text{left}}(\text{cap}) + \text{count}_{\text{right}}(\text{cap})}. \]

Medians Don’t Decompose

Using our other error criteria, we can break the data down the same way.

We get the medians on the left and right using mean absolute error.
We get the centers of the ranges on the left and right using max abs error.

But we cannot put them back together again. The overall median is not determined by medians on the left and right.

\[ \begin{aligned} &\text{median}\ \{1,2,3\} \cup \{4,5,5,6\} = 4; \\ &\text{median}\ \{1,2,3\} \cup \{0,5,5,6\} = 3. \end{aligned} \]

Can we approximate the overall median using the left and right medians?

No, not really. If you add a million to the points on the right in our two examples, the overall medians are $1{,}000{,}004$ and $3$ respectively.

What This Means

If you want to break stuff down and put it back together, you’ve only got one option. You’re going to be minimizing mean squared error. That means it’s what you do to estimate average treatment effects.

You break things down into categories in which you can make meaningful comparisons. You’d break schools down by student income, region, etc.
Then you put them back together into something relevant to policy, like the average effect of a class-size cap over all students in the state.

\[ \text{mean improvement} \approx \frac{\sum_{i=1}^n \text{e}_i \cdot \left\{ \hat\mu\left(\frac{\text{e}_i}{\lceil\text{e}_i/20\rceil}, \text{etc}_i\right) - \hat\mu\left(\frac{\text{e}_i}{\lceil\text{e}_i/40\rceil}, \text{etc}_i\right)\right\}}{\sum_{i=1}^n \text{e}_i}. \]

It also means that to estimate the median improvement is much harder. Not just in the sense that it takes more work, but that it’s less possible. Your estimate winds up being based on much stronger assumptions. That’s why we talk about mean improvements most of the time.

What About Other Models?

This means/medians/range-centers stuff is for constant models. And for piecewise constant ones. Two horizontal lines, three, etc. But I’d say that’s enough to favor least squares more generally. My reasoning is about what we’d do if we had all the data we could want. If we did, there’d be no reason to fit a restrictive model.

We’d fit a model that can take any value at any observable $x$. That model— the fully saturated model—is piecewise constant. When we have less data, we’re doing our best to approximate it. And if we’d want to use least squares in the fully saturated model, you’d think we’d want to use least squares in our approximations.

Beyond Linear Models

The least squares criterion $\frac{1}{n}\sum_{i=1}^n \{m(X_i)-Y_i\}^2$ doesn’t care whether $\mathcal{M}$ is a subspace or not. You can minimize over any set of functions.

Monotone regression: $\mathcal{M}$ = all non-decreasing functions. Still least squares, but the model isn’t a subspace (a sum of two non-decreasing functions is non-decreasing, but a negative scalar multiple isn’t).
Exponential models: $\mathcal{M} = \{m(x) = e^{a+bx}\}$. The curve is nonlinear in the parameters, but we can still minimize squared error over $a$ and $b$.

The point: least squares is the error criterion. Linear models are a choice of $\mathcal{M}$. They go together well, but they’re separate choices.

Models

Linear Models

\[ \hat \mu = \operatorname{argmin}_{m \in \mathcal{M}} \frac{1}{n}\sum_{i=1}^n \{ m(X_i) - Y_i \}^2 \]

In this class, we’ll focus on linear models. That’s what’s most common. A linear model is a subspace of the functions of $x$.

That means sums of functions in our model are also in our model. \[m_1, m_2 \in \mathcal{M} \implies m_1 + m_2 \in \mathcal{M}.\]
And scalar multiples of functions in our model are also in our model. \[m \in \mathcal{M} \implies \beta m \in \mathcal{M} \quad \text{ for all } \quad \beta \in \mathbb{R}.\]

Here we’re adding and scaling pointwise. \[(m_1 + m_2)(x) = m_1(x) + m_2(x) \quad \text{ and } \quad (\beta m)(x) = \beta m(x).\]

Basis Expansion

A subspace has a basis.

In particular, a subspace of functions has a basis of functions. \[\phi_0(x), \phi_1(x), \ldots, \phi_p(x)\]
Each function in the model is a linear combination of these basis functions. \[m(x)=\sum_{j=0}^p \beta_j \phi_j(x) \quad \text{ for all } m \in \mathcal{M}.\]
And our model, the subspace, is the set of all linear combinations of them. \[\mathcal{M} = \left\{ \sum_{j=0}^p \beta_j \phi_j(x) \ \mid \ \beta \in \mathbb{R}^p \right\}\]

That means we can find the best function in our model by optimizing over vectors of coefficients $\beta \in \mathbb{R}^p$. \[ \begin{aligned} \hat \mu &= \operatorname{argmin}_{m \in \mathcal{M}}\ \frac{1}{n}\sum_{i=1}^n \left\{ m(X_i) - Y_i \right\}^2 \\ &= \sum_{j=0}^p \hat\beta_j \phi_j \ \text{ where } \ \hat \beta = \operatorname{argmin}_{b \in \mathbb{R}^p} \ \frac{1}{n}\sum_{i=1}^n \left\{ \sum_{j=0}^p b_j \phi_j(X_i) - Y_i \right\}^2. \end{aligned} \]

Examples

All constants

\[\mathcal{M}=\left\{ m(x)= a \quad \mid \quad a \in \mathbb{R} \right\}\]

This has a basis made up of one function: $\phi_0(x)=1$.

All functions that are constant except at $x=0$

\[\mathcal{M} = \left\{ m(x) = \begin{cases} a & \text{ if } x \le 0 \\ b & \text{ if } x > 0 \end{cases} \quad \middle| \quad a,b \in \mathbb{R} \right\}\]

This has a basis made up of two functions. \[\phi_0(x)=\begin{cases} 1 & x \le 0 \\ 0 & x > 0 \end{cases}, \qquad \phi_1(x)=\begin{cases} 0 & x \le 0 \\ 1 & x > 0 \end{cases}\]

$\beta_0$ is the function’s value on the left, $\beta_1$ its value on the right.

Sometimes we prefer this basis: \[\tilde \phi_0(x)=\phi_0(x) + \phi_1(x) = 1, \qquad \tilde \phi_1 = \phi_1(x)-\phi_0(x) = \operatorname{sign}(x).\]

$\beta_0$ is the function’s value on the left, $\beta_1$ the increment to its value on the right.

All lines

\[\mathcal{M}=\left\{ m(x)= a + bx \quad \mid \quad a,b \in \mathbb{R} \right\}\]

This has a basis made up of two functions: $\phi_0(x)=1, \quad \phi_1(x)=x$.

All cubic polynomials

\[\mathcal{M} = \left\{ m(x) = a+bx+cx^2 + dx^3 \quad \mid \quad a,b,c,d \in \mathbb{R} \right\}\]

This has a basis made up of four functions: $\phi_j(x)=x^j$ for $j=0\ldots 3$.

Are These Linear Models?

All non-negative constants

\[\mathcal{M}=\left\{ m(x)= a \quad \mid \quad a \ge 0 \right\}\]

This is not a linear model. If we scale by a negative constant, we leave the model.

All lines through the origin

\[\mathcal{M}=\left\{ m(x)= a + bx \quad \ \mid \quad a,b \in \mathbb{R} \ \ \text{ with } \ \ m(0)=0 \right\}\]

This is a linear model. It has a basis made up of one function: $\phi_0(x) = x$.

All lines through the point $(0,1)$

\[\mathcal{M}=\left\{ m(x)= a + bx \quad \mid \quad a,b \in \mathbb{R} \ \ \text{with} \ \ m(0)=1 \right\}\]

This is not a linear model. If we add two, we get a line through $(0,2)$.

All lines through the point $(1,0)$

\[\mathcal{M}=\left\{ m(x)= a + bx \quad \mid \quad a,b \in \mathbb{R} \ \ \text{with} \ \ m(1)=0 \right\}\]

This is a linear model. It has a basis made up of one function: $\phi_0(x) = x-1$.

All functions that are constant except at a single point

\[\mathcal{M} = \left\{ m(x) = \begin{cases} a & \text{ if } x \le c \\ b & \text{ if } x > c \end{cases} \quad \middle| \quad a,b,c \in \mathbb{R} \right\}\]

This is not a linear model. If we add two, we’ll usually get a function that’s constant except at two points.

All continuous functions that are linear to the left and right of $x=0$

\[\mathcal{M} = \left\{ \text{continuous}\ \ m(x) = \begin{cases} a+bx & \text{ if } x \le 0 \\ c+dx & \text{ if } x > 0 \end{cases} \quad \middle| \quad a,b,c,d \in \mathbb{R} \right\}\]

This is a linear model. \[\phi_0(x)=1,\ \ \phi_1(x)=x, \ \ \phi_2(x)=x \operatorname{sign}(x)\]

Least Squares in Linear Models

Residuals and Orthogonality

When we fit a curve $\hat\mu$, we decompose our observations into two parts.

The predictable part, the value of our fitted curve at the observed $x$.
The residual, the difference between that value and what we observe.

\[Y_i = \underset{\text{predictable}}{\hat \mu(X_i)} + \underset{\text{residual}}{\hat \varepsilon_i}.\]

When we use least squares with a linear model, we can think of these as projections onto the model and its orthogonal complement respectively. This tells us that the residuals are orthogonal to the model.

\[ \begin{aligned} 0 &= \frac{\partial}{\partial b_j} \frac{1}{n}\sum_{i=1}^n \left\{ \sum_{j=0}^p b_j \phi_j(X_i) - Y_i \right\}^2 \bigg\vert_{b=\hat\beta} \\ &= \frac{2}{n}\sum_{i=1}^n \left\{ \sum_{j=0}^p \hat\beta_j \phi_j(X_i) - Y_i \right\}\phi_j(X_i) \\ &= -\frac{2}{n}\sum_{i=1}^n \hat\varepsilon_i \phi_j(X_i) \quad \text{ for all } j \ \text{ and therefore } \\ 0 &= \frac{1}{n}\sum_{i=1}^n \hat\varepsilon_i \left\{\sum_{j=0}^p \beta_j \phi_j(X_i) \right\} \quad \text{ for all $\beta$} \\ &= \frac{1}{n}\sum_{i=1}^n \hat\varepsilon_i m(X_i) \quad \text{ for all $m \in \mathcal{M}$}. \end{aligned} \]

Implications

\[0 = \frac{1}{n}\sum_{i=1}^n \hat\varepsilon_i m(X_i) \quad \text{ for all } \quad m \in \mathcal{M}\]

This generalizes what we said about least squares in the constant model.

If our model includes constants, our residuals will have mean zero.
- That is, our fitted curve will have the same mean as the data.
- It won’t necessarily be the mean, but it’ll get the mean right.
This extends naturally to models including piecewise constant curves.
- Our fitted curve will have the same mean as the data on each piece.
And if our model includes lines, our residuals will have no linear trend.
- That is, our fitted curve will rise or fall the right amount on average.
One way to think about this: if there’s a trend we could fit, we’d fit it.

This doesn’t mean our curve fits the data, but it limits the ways it doesn’t.