--- title: Trends, Smoothing, Detrending author: 36-467/36-667 date: 30 August 2018 bibliography: locusts.bib output: html_document: toc: true number_sections: true toc_depth: 4 --- \[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator{\dof}{DoF} \DeclareMathOperator{\det}{det} \] # In our last episode... Our data will generally be a realization of a **stochastic process**, i.e., a collection of random variables spread over time and space. The random variable we observe at spatial coordinate $r$ and time $t$ will be, generically, $X(r,t)$. If we just care about spatial data, we'll just have $X(r)$, and purely temporal data (like the Kyoto data) will be $X(t)$. In general, the random variables at different points will all be dependent on each other. Because they are random variables, it makes sense to think about their central tendency, and the expected value or mean is often a convenient one: \[ \TrueRegFunc(r,t) = \Expect{X(r,t)} \] This is a _deterministic_ function of the coordinates; it's what we mean by the **trend**. Mathematically, we can always decompose the deterministic random variables into the trend plus random or stochastic **fluctuations**, which always have expectation zero: \begin{eqnarray} X(r,t) & = & \TrueRegFunc(r,t) + \TrueNoise(r,t)\\ \Expect{\TrueNoise(r,t)} & = & 0 \end{eqnarray} We can't assert anything more about $\TrueNoise$ without further assumptions. In particular, we can't make the usual assumptions that are made in regression courses, that $\TrueNoise$ is uncorrelated, that it has constant variance, that it has the same distribution at each point, that it's Gaussian, etc. But we can, without any loss of generality, say $\Expect{\TrueNoise(r,t)} = 0$ everywhere. # Learning the Trend We would like to learn the trend function $\TrueRegFunc$. ## Many realizations If we have _many_ independent realizations of the process, we could just average them. Say we can re-run some experiment many times, getting independent realizations $X^{(1)}, X^{(2)}, \ldots X^{(m)}$. Then \[ \frac{1}{m}\sum_{i=1}^{m}{X^{(i)}(r,t)} \rightarrow \TrueRegFunc(r,t) \] by the ordinary law of large numbers. As my mention of "experiment" suggests, this idea _can_ be applied in the experimental sciences. In neuroscience, for instance, we might repeatedly measure what happens to a particular neuron (or group of neurons, or brain region) after we apply some stimulus to an experimental animal[^whiskers]. This is, however, not available in non-experimental contexts. [^whiskers]: I once analyzed an experiment where a mechanical vibrator repeatedly gave very precise tweaks to a rat's whiskers. (Rats devote a _lot_ of their brain to processing sensations from their whiskers, which help them navigate in the dark.) ## Single realization More typically, we have only _one_ realization of the $X$ process --- only one history of cherry-blossoms in Kyoto, only one map of ground-water under Wisconsin, etc. We will only be able to learn about $\TrueRegFunc$ under assumptions. ### Trends from theory In some situations, we might have a well-confirmed scientific theory which tells us _exactly_ what the trend should be, say that $\TrueRegFunc(r,t) = f(r,t)$ for some completely-specified function $f$. There are very few situations where this applies --- _perhaps_ in parts of astronomy and physics. More common is the situation where theory pins down a specific functional _form_ for the trend, but there are some unknown coefficients in the formula. That is, $\TrueRegFunc(r,t) = f(r,t,\theta)$ for one or more parameters $\theta$. One would then need to estimate $\theta$ from the data $X$ to obtain an estimate of the trend. The dependence of the observations on each other makes the statistical inference more complicated, which we'll come back to later in the course, but in principle this is just about estimating parameters. Of course, for these approaches to be good ideas, you have to actually have reason to think your theory really does give you a good idea of the trend. There are domains where this is true --- physics, chemistry, some parts of biology (especially genetics, evolution and ecology) --- where we have very well-supported theories which are detailed enough to pin down specific functional forms. (Some economists think this is true in economics as well, but I think they're fooling themselves.) Beyond that, though, scientific theories are much looser -- they might tell us that one variable tends to go up with another, but just don't have the detail to give a functional form. This might not be a permanent condition --- it used to be true even of physics! --- but for now, in lots of scientific areas, let alone in business or policy, we just don't have credible theories which tell us about trends. ### Trends from data We might, however, make a different sort of assumption, which is just that $\TrueRegFunc$ doesn't change very rapidly --- that it is a **smooth** function. The basic idea we need is that $\TrueRegFunc(t-h) \approx \TrueRegFunc(t) \approx \TrueRegFunc(t+h)$, at least for sufficiently small $h$, at least at most $t$. (You can write out the corresponding idea for spatio-temporal trends, with some more notation.) Each $X(t)$ is a noisy observation of $\TrueRegFunc(t)$, but if $\TrueRegFunc(t-h)$ and $\TrueRegFunc(t+h)$ are close to $\TrueRegFunc(t)$, then $X(t-h)$ and $X(t+h)$ will be noisy observations of things close to $\TrueRegFunc(t)$, giving us more information. To make this concrete, and see how it can work, let's think of a very crude estimate: \[ \EstRegFunc(t_i) = \frac{1}{3}\sum_{k=-1}^{k=+1}{X(t_{i+k})} \] Each $X$ is trend plus fluctuation, so \begin{eqnarray} \EstRegFunc(t_i) & = & \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueRegFunc(t_{i+k})} + \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueNoise(t_{i+k})}\\ & = & \TrueRegFunc(t_i) + \frac{1}{3}\sum_{k=-1}^{k=+1}{(\TrueRegFunc(t_{i+k}) - \TrueRegFunc(t_i))} + \frac{1}{3}\sum_{k=-1}^{k=+1}{\TrueNoise(t_{i+k})}\\ & = & \text{truth} + \text{bias} + \text{noise} \end{eqnarray} The first term is exactly what we want. The third term is noise, and the average of many noise terms is _generally_ smaller than the individual noises. (If the $\TrueNoise$ were uncorrelated with equal variance $\sigma^2$, this would have variance $\sigma^2/3$.) The middle term is **bias**, systematic error. The more rapidly $\TrueRegFunc$ changes, and the further apart our observations are located, big this bias gets. Conversely, if $\TrueRegFunc$ changes slowly and our observations are closely spaced, the bias will be small. At this point, we _could_ go into considerable detail about exactly what the bias is, under various precise mathematical characterizations of "smooth function". (For instance, we could use Taylor approximation[^taylor] to relate the bias to the first and second derivatives of $\TrueRegFunc$.) But it'll be more productive to think about the more general _type_ of estimator which includes the simple averaging of near-by values. There is, after all, no reason to just use the immediate neighbors, nor to use equal weights. [^taylor]: As you remember from calculus, the Taylor expansion of a function $f$ around a point $t_0$ is $f(t) = f(t_0) + ((t-t_0) f^{\prime}(t_0) + \frac{1}{2}(t-t_0)^2 f^{\prime\prime}(t_0) + \ldots$, where $f^{\prime}$, $f^{\prime\prime}$, etc., are the first, second, etc., derivatives. In a Taylor approximation, we truncate the series at a certain order. The error of a $k^{\mathrm{th}}$ order approximation is $O((t-t_0)^{k+1})$. # Linear smoothers We have data points $x(t_1), \ldots x(t_n)$; for conciseness I will sometimes abbreviate these as $x_1, \ldots x_n$. The general form of a **linear smoother** is \[ \EstRegFunc(t) = \sum_{j=1}^{n}{w(t, t_j) x_j} \] That is, the predictions or fitted values, the estimates of the trend $\TrueRegFunc$, are always linear combinations of the observed values of $x$. Notice that the smoother is linear in $x$, not in $t$ --- the resulting curve can be _very_ nonlinear in $t$. ## Examples of linear smoothers Most of the smoothing, regression and estimation techniques used in practice fall in this class. - _Global constant_: The trivial case, $w=1/n$ no matter what. - $k$ _Nearest neighbors_ (kNN): $w(t, t_j) = 1/k$ if $t_j$ is one of the $k$ points closest to $t$, and otherwise $w=0$. - _Moving average_ (MA): $w(t, t_j) = 0$ outside some distance from $t$, and constant within it. + If data points are always equally spaced in time (or space, etc.), moving averages and nearest neighbors are equivalent, but they can differ for irregularly-spaced data. + You can also do one-sided moving averages, usually for the past. - _Weighted moving averages_: $w(t, t_j) = f(|t-t_j|/\tau)$, for some fixed function $f$. A common choice is the exponential function, giving **exponentially-weighted moving averages**. - _Kernels_: Like moving averages, but with weights normalized, so we're always doing a proper average. That is, one picks a probability-density function $K$, and some scale, or **bandwidth**, $h$, and \[ w(t, t_j) = \frac{K(|t-t_j|/h)}{\sum_{j^{\prime}=1}^{n}{K(|t-t_{j^{\prime}}|/h)}} \] - _Linear and polynomial regression_: It's not obvious, but these too are linear smoothers, though the weights are weird. (See the optional section on linear regression below.) - _Splines_: I teased you with these last time. They're smooth, piecewise polynomials, and so they, too, are linear smoothers. (We'll go over the exact formulas below.) ## Properties of linear smoothers The fact that lots of statistical methods can be seen as linear smoothers is useful because it gives us a unified way of understanding how these methods work, and in particular of when they will do a good job at recovering trends. Before establishing these properties, I need to introduce one more bit of notation. At each $t_i$ where we took an observation, we have a fitted value $\EstRegFunc(t_i)$, or for short $\EstRegFunc_i$. Group these together into a vector, or $n\times 1$ matrix, \[ \mathbf{\EstRegFunc} = \left[\begin{array}{c} \EstRegFunc(t_1) \\ \EstRegFunc(t_2) \\ \vdots \\ \EstRegFunc(t_n) \end{array} \right] \] Similarly, group all the observations of $x$ into a vector $\mathbf{x}$. Then the fitted vectors are the observations _times_ a matrix: \[ \mathbf{\EstRegFunc} = \mathbf{w} \mathbf{x} \] where $w_{ij} = w(t_i, t_j)$. ### Expectation of the fitted values; bias The first thing this lets us get at is the expected value of the fitted values: \begin{eqnarray} \Expect{\mathbf{\EstRegFunc}} & = & \Expect{\mathbf{w}\mathbf{X}}\\ & = & \mathbf{w}\Expect{\mathbf{X}}\\ & = & \mathbf{w}\Expect{\mathbf{\TrueRegFunc} + \mathbf{\TrueNoise}}\\ & = & \mathbf{w}\Expect{\mathbf{\TrueRegFunc}} + \mathbf{w}\Expect{\mathbf{\TrueNoise}}\\ & = & \mathbf{w}\mathbf{\TrueRegFunc} \end{eqnarray} because $\mathbf{\TrueRegFunc}$ is a constant vector, and the noise $\TrueNoise$ always has expectation 0. The **bias** of any estimator is the difference between its expected value and the truth. The estimator is **unbiased** when that difference is zero. So we have a simple equation for when a smoother gives us an unbiased estimate of the trend: \begin{equation} \text{unbiased trend estimates} \Leftrightarrow \mathbf{\TrueRegFunc} = \mathbf{w}\mathbf{\TrueRegFunc} \end{equation} This says that we get an unbiased estimate if, and only if, $\mathbf{\TrueRegFunc}$ is an eigenvector of $\mathbf{w}$ with eigenvalue 1. (See below for some reminders of the necessary linear algebra.) Otherwise, we get some bias in our estimate. ### Shrinkage Suppose (what is generally true) that $\mathbf{w}$ has $n$ linearly independent eigenvectors, say $\mathbf{v}_1, \ldots \mathbf{v}_n$, with corresponding eigenvalues $\lambda_1, \ldots \lambda_n$. These eigenvectors form a basis for space of vectors, so for an arbitrary vector $\mathbf{x}$, \[ \mathbf{x} = \sum_{i=1}^{n}{c_i \mathbf{v}_i} \] for some coordinates $c_1, \ldots c_n$. These coordinates say how much $\mathbf{x}$ projects on to the eigenvectors, i.e., how similar it is to each of the patterns those eigenvectors represent. Now what happens when we multiply by $\mathbf{w}$, i.e., when we apply our favorite linear smoother to $\mathbf{x}$? \[ \mathbf{w}\mathbf{x} = \sum_{i=1}^{n}{\mathbf{w} c_i \mathbf{v}_i} = \sum_{i=1}^{n}{c_i \lambda_i \mathbf{v}_i} \] The components of $\mathbf{x}$ which match eigenvectors $\mathbf{v}_i$ with large eigenvalues $\lambda_i$ get enhance, and those which match eigenvectors with small eigenvalues get **shrunk**. Most of the time, when we're dealing with smoothers, the largest eigenvalue is exactly 1. The eigenvectors with eigenvalue 1 define an **invariant** subspace, one which is left alone by the smoother. Because all the other eigenvalues are smaller, smoothing _shrinks_ the part of the data which doesn't conform to the invariant subspace. If the eigenvalues which are less than 1 are in fact exactly 0, every component of the data which doesn't fall in the invariant subspace is projected away, and the vector of fitted values $\mathbf{\EstRegFunc}$ has to fall in the span of the leading eigenvectors. If the smaller eigenvalues aren't quite 0, the shrinkage isn't quite so extreme, and $\mathbf{\EstRegFunc}$ doesn't have to be in the invariant subspace, but it is always closer to that space than the data was. 1. If all the weights are non-negative, $w_{ij} \geq 0$, and the weights for each fitted value sum to 1, $\sum_{j}{w_{ij}} = 1$ for all $i$, then the largest eigenvalue of $\mathbf{w}$ is exactly 1, and all the other eigenvalues are $\leq 1$ in magnitude. So when the linear smoother really does proper averages to get each fitted value, the vector of fitted values is shrunk towards the invariant subspace. 2. For linear regression, every eigenvalue of $\mathbf{w}$ is either 0 or 1. Hence vector of fitted values is always in the invariant subspace. #### A little example Let's go back to the idea of averaging each observation with the one ahead and behind it. We need to decide what to do with the first and last observations; I'll set them to weights of 1/2 ```{r} n <- 10 w <- matrix(0, nrow=10, ncol=10) diag(w) <- 1/3 for (i in 2:(n-1)) { w[i,i+1] <- 1/3 w[i,i-1] <- 1/3 } w[1,1] <- 1/2 w[1,2] <- 1/2 w[n,n-1] <- 1/2 w[n,n] <- 1/2 ``` (You can check that the `w` this produces has the right entries.) Here are the eigenvalues: ```{r} eigen(w)$values ``` Notice, as promised, that 1 is an eigenvalue, and all the others are smaller in magnitude. The corresponding eigenvector is ```{r} eigen(w)$vectors[,1] ``` i.e., _constant_. This should make sense --- if every entry in $\mathbf{x}$ was the same, averaging wouldn't change anything. The next few eigenvectors are better grasped by visualization than by looking at the numbers. ```{r} plot(1:10, eigen(w)$vectors[,2], lty="solid", type="l", xlab="Position", ylab="Eigenvectors of smoothing matrix") lines(1:10, eigen(w)$vectors[,3], lty="dashed") lines(1:10, eigen(w)$vectors[,4], lty="dotted") legend("bottomleft", legend=c(expression(j==2), expression(j==3), expression(j==3)), lty=c("solid","dashed","dotted")) ``` These look somewhat like sine waves, because they _are_ somewhat like sine waves[^movingaveragenadsinewaves]. Because the eigenvalues which go with these waves are pretty close to 1, the components of the original data which match these sine waves will be shrunk just a little, while the other components of the data will be shrunk a lot. So these (plus the constant) are the sorts of patterns which this moving average will preserve; everything else it tends to ignore. [^movingaveragenadsinewaves]: We will come back to this connection, between moving averages and sine waves, when we look at Fourier analysis later in the course. For now, you can take it on trust. --- If you are unwilling to take it on trust, think about how you'd approximate a 2nd derivative with finite differences: $\frac{\frac{x(t+h) - x(t)}{h} - \frac{(x(t) - x(t-h)}{h}}{h} = \frac{x(t+h) - 2 x(t) + x(t-h)}{h^2}$. So replacing $x(t)$ with the moving average is approximately the same as replacing $x(t)$ with $x(t) + \kappa x^{\prime\prime}(t)$ for some scaling factor $\kappa$. And you can verify for yourself that every sine or cosine function is an eigenfunction of this differential operator, i.e., that $\left(1+\kappa\frac{d^2}{dt^2}\right)\sin{\omega t} = \lambda\sin{\omega t}$. (Can you find $\lambda$?) So it should make some sense that the eigenvectors of the smoothing matrix are close to the eigenfunctions of the differential operator. ### Variance of the fitted values Because the data contains fluctuations, which are random, it makes sense to worry about the variance and covariance of the fitted values. These are easy: \begin{eqnarray} \Var{\mathbf{\EstRegFunc}} & = & \Var{\mathbf{w}\mathbf{X}}\\ & = & \mathbf{w}\Var{\mathbf{X}}\mathbf{w}^T\\ & = & \mathbf{w}\Var{\mathbf{\TrueRegFunc} + \mathbf{\TrueNoise}}\mathbf{w}^T\\ & = & \mathbf{w}\Var{\mathbf{\TrueNoise}}\mathbf{w}^T \end{eqnarray} Of course, we don't necessarily know what $\Var{\mathbf{\TrueNoise}}$ is. In other classes, we often assume that all of the noise terms have equal variance $\sigma^2$, and are uncorrelated with each other. That means $\Var{\mathbf{\TrueNoise}} = \sigma^2 \mathbf{I}$, and we get the special case \[ \Var{\EstRegFunc} = \sigma^2 \mathbf{w}\mathbf{w}^T \] ### More about shrinkage; degrees of freedom As the data randomly vary, the fitted values are going to randomly vary in response. One way to get a sense of how sensitive the fitted values are to the data is to look at the sum of the covariances, \begin{eqnarray} \sum_{i=1}^{n}{\Cov{\EstRegFunc_i, X_i}} & = & \sum_{i=1}^{n}{\Cov{\sum_{j=1}^{n}{w_{ij} X_j}, X_i}}\\ & = & \sum_{i=1}^{n}{\sum_{j=1}^{n}{w_{ij} \Cov{X_i, X_j}}}\\ & = & \sum_{i=1}^{n}{\sum_{j=1}^{n}{w_{ij} \Cov{\TrueNoise_i, \TrueNoise_j}}} \end{eqnarray} Again, in the special case where $\Cov{\TrueNoise_i, \TrueNoise_j} = 0$ if $i\neq j$, and $=\sigma^2$ if $i=j$, this simplifies to \[ \sum_{i=1}^{n}{\Cov{\EstRegFunc_i, X_i}} = \sigma^2 \tr{\mathbf{w}} \] Remembering that the trace of a matrix equals the sum of its eigenvalues, we use this to _define_ the **effective number of degrees of freedom**: \[ \dof{\mathbf{w}} \equiv \tr{\mathbf{w}} \] If $\mathbf{w}$ is an idempotent projection matrix, this is exactly equal to the number of dimensions we're projecting down to. If $\mathbf{w}$ isn't a projection matrix, we pick up the dimension of the invariant space (= the number of $\lambda=1$ eigenvalues), plus a bit more. ## Splines: Direct control over smoothness I have teased you several times by now by mentioning **splines**, or more formally **smoothing splines**[^splines]. These are smooth curves which we find by balancing coming close to the data points, against not having too much wiggliness or curvature. More specifically, we pick a $\lambda > 0$, and then solve the optimization problem \[ \EstRegFunc = \argmin_{m}{\frac{1}{n}\sum_{i=1}^{n}{(x_i - m(t_i))^2} + \lambda\int{(m^{\prime\prime}(t))^2 dt}} \] The first term we're optimizing over is the mean squared error. The second term is the average curvature of the function. It would be zero for any linear function (whose second derivative is 0 everywhere), and it grows the more the function bends and twists around. One way to read this is that we are **penalizing** curvature. In fact, you can think of $\lambda$ as a price --- we're willing to trade one unit more of curvature for at least $\lambda$ units of reduced MSE. (Conversely, we're willing to increase the MSE by one unit, if in doing so we reduce the curvature by at least $1/\lambda$ units.) As $\lambda \rightarrow 0$, we get incredibly wiggly functions which go through the data points exactly. As $\lambda \rightarrow \infty$, we approach the best-fitting straight line, i.e., linear regression. [^splines]: People have defined many variants on splines. When I use the word without an adjective, I mean smoothing splines, as defined by the optimization problem immediately following this footnote. The optimization here is over all possible functions with second derivatives. You might well ask how we could possibly do this, since it's an infinite-dimensional space. The answer is that the solution is _always_ a piecewise cubic polynomial, with pieces at the data points $t_i$, and that the solution is also continuous, with continuous first and second derivatives. (See the optional section below, if you care.) All of this means that we only have $n$ unknown coefficients to find, rather than having to search over all possible functions. (Again, see the optional section below.) One upshot of being a piecewise polynomial is that the spline is a linear smoother, with weights which do a kind of local averaging. Another is that there is a direct correspondence between $\lambda$ and the effective degrees of freedom, where increasing $\lambda$ lowers the degrees of freedom and vice versa. (The exact relationship depends on the values of $t_i$ --- again, see the optional section.) ### Multi-dimensional splines All of the above is for one-dimensional smoothing. If we're dealing with data over two dimensions, so it's really $x(r,t)$ rather than just $x(t)$, we can define the analogous problem: \[ \argmin_{m}{\frac{1}{n}\sum_{i=1}^{n}{(x_i - m(r_i, t_i))^2} + \lambda\int{\left[ {\left(\frac{\partial^2 m}{dr^2}\right)}^2 + 2{\left(\frac{\partial^2 m}{dr dt}\right)}^2 + {\left(\frac{\partial^2 m}{dt^2}\right)}^2\right] dr dt}} \] (Three or more coordinates is similar, but with even more notation.) These are called **thin-plate** splines. An alternative approach is to use what are called a "tensor-product splines". If you dive into the math (see the optional section), you find that one-dimensional splines are linear combinations of "basis functions" $B_1(t), \ldots B_n(t)$ --- the exact shape of the basis functions depends on the $t_i$ but not on the $x_i$. In two dimensions, we can find the basis functions for the two coordinates, say $B_i$ for $r$ and $C_j$ for $t$, and take products, so \[ m(r,t) = \sum_{i=1}^{n}{\sum_{j=1}^{n}{\beta_{ij} B_i(r) C_j(t)}} \] - Thin-plate splines have the nice interpretation of trading curvature against ME, just like in 1D. Tensor-product splines are harder to interpret. + On the other hand, this itself makes more sense when the coordinates have similar scales and meaning --- if $r$ and $t$ are both map coordinates, total curvature is a lot more meaningful than if $r$ is space and $t$ is time. - Tensor-product splines are usually easier to calculate than thin-plate splines. ### Splines in R The simplest way to work with splines in R is to use the `smooth.spline` command. This takes two vectors, `x` (corresponding to our $t$ variable) and `y` (corresponding to our $x$ variable), and returns a slightly complicated object. The `x` and `y` components contain the `x` values, and the _fitted_ `y` values, so it can be used for plotting. (See the example in the slides for lecture 1.) You can also calculate predictions at _new_ values, not corresponding to the data. For more sophisticated uses, especially multi-dimensional splines, you need to use a package. I recommend the `mgcv` package [@Wood-mgcv-2004,@Wood-on-GAMs], which is primarily designed for smooth, additive models, and as such includes a lot of different smoothers. A command like ```{r, eval=FALSE} library(mgcv) my.fit <- gam(x ~ s(t), data=df) ``` estimates a regression model where `x` (drawn from data frame `df`) is a smooth spline function of `t` (drawn from the same source), and stores the result in the object `my.fit`. On the other hand, ```{r, eval=FALSE} gam(x ~ s(r,t), data=df) ``` models `x` as a thin-plate spline over the two variables `r` and `t`, while ```{r, eval=FALSE} gam(x ~ te(r, t), data=df) ``` uses a tensor-product spline. If you want to model `x` as the _sum_ of two smooth functions, that would be ```{r, eval=FALSE} gam(x ~ s(r) + s(t), data=df) ``` There are also a large number of more advanced packages exclusively for dealing with splines --- too many for me to list here. ### Hodrick-Prescott filtering In econometrics, it is very common to use what is called a "Hodrick-Prescott filter" to extract the trend of a time series. (The canonical reference is @Hodrick-Prescott-on-HP-filter, but this circulated as a working paper for many years before being formally published.) As shown by @Paige-Trindade-on-HP-filter, this is _exactly_ the same as applying a smoothing spline to the data. The only difference is that Hodrick and Prescott supplied a (fallacious) argument for why $\lambda$ should always be exactly 1600. ### Some mathematical details for splines > This section is optional, but a lot of it is cool. The text is largely ripped off from @CRS-ADAfaEPoV. #### Solutions to the spline problem are piecewise-cubic polynomials The spline problem asks us to find the function which minimize the sum of the MSE and a certain integral. Even the MSE can be brought inside the integral, using Dirac delta functions[^dirac]: \begin{equation} \mathcal{L} = \int{\left[\lambda(m^{\prime\prime}(t))^2 + \frac{1}{n} \sum_{i=1}^{n}{(x_i - m(t_i))^2 \delta(t-t_i)}\right] dt} \end{equation} In what follows, without loss of generality, assume that the $t_i$ are ordered, so $t_1 \leq t_2 \leq \ldots t_i \leq t_{i+1} \leq \ldots t_n$. With some loss of generality but a great gain in simplicity, assume none of the $t_i$ are equal, so we can make those inequalities strict. [^dirac]: The Dirac delta function (named after the physics P. A. M. Dirac) is defined by how it integrates: for any function $f(t)$, $\int{f(t) \delta(t) dt} = f(0)$. You can think of this as the limit of a sequence of Gaussian pdfs, all centered at 0, with variances shrinking towards zero. The limit would be something like a function $\delta$ which was $0$ except at the origin, but had an infinitely skinny, infinitely tall spike at the origin, with the total area under the curve being exactly 1. There isn't really a _function_ like that, but it's definitely true that $\lim_{\sigma\rightarrow 0}{\int{f(t) \frac{1}{\sigma}phi(t/\sigma) dt}} = f(0)$, so $\delta(t)$ makes sense _in integrals_. One can in fact build a whole theory of "generalized functions" which make sense in integrals in this way, including actual functions as a special case. If this interests you, @Kolmogorov-Fomin-analysis is still one of the best ways in to the subject. (Dirac, being a physicist, didn't let something like $\delta$ not really being a function stop him from calculating.) The subject which deals with maximizing or minimizing integrals of functions is the calculus of variations, and one of its basic tricks is to write the integrand as a function of $x$, the function, and its derivatives: \begin{equation} \mathcal{L} = \int{L(t,m,m^{\prime},m^{\prime\prime}) dt} \end{equation} where, in our case, \begin{equation} L = \lambda(m^{\prime\prime}(t))^2 + \frac{1}{n} \sum_{i=1}^{n}{(x_i - m(t_i))^2 \delta(t-t_i)} \end{equation} This sets us up to use a general theorem of the calculus of variations, to the effect that any function $\hat{m}$ which minimizes $L$ must also solve $L$'s **Euler-Lagrange equation**: \[ \left.{\frac{\partial L}{\partial m} - \frac{d}{dt}\frac{\partial L}{\partial m^{\prime}} + \frac{d^2}{dt^2}\frac{\partial L}{\partial m^{\prime\prime}}}\right|_{m=\hat{m}} = 0 \] In our case, the Euler-Lagrange equation reads \[ -\frac{2}{n}\sum_{i=1}^{n}{(x_i-\hat{m}(t_i))\delta(t-t_i)} + 2\lambda\frac{d^2}{dt^2}\hat{m}^{\prime\prime}(t) = 0 \] Remembering that $\hat{m}^{\prime\prime}(t) = d^2\hat{m}/dt^2$, \[ \frac{d^4}{dt^4}\hat{m}(t) = \frac{1}{n\lambda}\sum_{i=1}^{n}{(x_i-\hat{m}(t_i))\delta(t-t_i)} \] The right-hand side is zero at any point $t$ other than one of the $t_i$, so the fourth derivative has to be zero in between the $t_i$. This in turn means that the function must be piecewise cubic. Now fix an $t_i$, and pick any two points which bracket it, but are both greater than $t_{i-1}$ and less than $t_{i+1}$; call them $l$ and $u$. Integrate our Euler-Lagrange equation from $l$ to $u$: \begin{eqnarray} \int_{l}^{u}{\frac{d^4}{dt^4}\hat{m}(t) dt} & = & \int_{l}^{u}{\frac{1}{n\lambda}\sum_{i=1}^{n}{(x_i-\hat{m}(t_i))\delta(t-t_i)}}\\ \hat{m}^{\prime\prime\prime}(u) - \hat{m}^{\prime\prime\prime}(l) &= & \frac{x_i - \hat{m}(t_i)}{n\lambda} \end{eqnarray} That is, the third derivative makes a jump when we move across $t_i$, though (since the fourth derivative is zero), it doesn't matter which pair of points above and below $t_i$ we compare third derivatives at. Integrating the equation again, \begin{equation} \hat{m}^{\prime\prime}(u) - \hat{m}^{\prime\prime}(l) = (u-l)\frac{x_i - \hat{m}(t_i)}{n\lambda} \end{equation} Letting $u$ and $l$ approach $x_i$ from either side, so $u-l \rightarrow 0$, we see that $\hat{m}^{\prime\prime}$ makes no jump at $t_i$. Repeating this trick twice more, we conclude the same about $\hat{m}^{\prime}$ and $\hat{m}$ itself. In other words, $\hat{m}$ must be continuous, with continuous first and second derivatives, and a third derivative that is constant on each $(t_i, t_{i+1})$ interval. Since the fourth derivative is zero on those intervals (and undefined at the $t_i$), the function must be a piecewise cubic, with the piece boundaries at the $t_i$, and continuity (up to the second derivative) across pieces. To see that the optimal function must be linear below $t_1$ and above $t_n$, suppose that it wasn't. Clearly, though, we could reduce the curvature as much as we want in those regions, without altering the value of the function at the boundary, or even its first derivative there. This would yield a better function, i.e., one with a lower value of $\mathcal{L}$, since the MSE would be unchanged and the average curvature would be smaller. Taking this to the limit, then, the function must be linear outside the observed data range. #### Finding the spline coefficients Splines, I said, are piecewise cubic polynomials. To see how to fit them, let's think about how to fit a global cubic polynomial. We would define four **basis functions**, \begin{eqnarray} B_1(x) & = & 1\\ B_2(x) & = & t\\ B_3(x) & = & t^2\\ B_4(x) & = & t^3 \end{eqnarray} and chose to only consider regression functions that are linear combinations of the basis functions, \begin{equation} m(t) = \sum_{j=1}^{4}{\beta_j B_j(t)} \end{equation} Such regression functions would be linear in the _transformed_ variables $B_1(t), \ldots B_4(t)$, even though it is nonlinear in $t$. To estimate the coefficients of the cubic polynomial, we would apply each basis function to each data point $t_i$ and gather the results in an $n\times 4$ matrix $\mathbf{b}$, \begin{equation} b_{ij} = B_j(t_i) \end{equation} Then we would do ordinary least squares using the $\mathbf{b}$ matrix in place of the usual data matrix $\mathbf{t}$: \begin{equation} \hat{\beta} = (\mathbf{b}^T\mathbf{b})^{-1}\mathbf{b}^T \mathbf{x} \end{equation} Since splines are piecewise cubics, things proceed similarly, but we need to be a little more careful in defining the basis functions. Recall that we have $n$ values of $t$, $t_1, t_2, \ldots t_n$, and we can assume they're in increasing order. These $n$ **knots** define $n+1$ pieces or segments: $n-1$ of them between the knots, one from $-\infty$ to $t_1$, and one from $t_n$ to $+\infty$. A third-order polynomial on each segment would seem to need a constant, linear, quadratic and cubic term per segment. So the segment running from $t_i$ to $t_{i+1}$ would need the basis functions \begin{equation} \mathbf{1}_{(t_i,t_{i+1})}(t), ~ (t-t_i)\mathbf{1}_{(t_i,t_{i+1})}(t), ~ (t-t_i)^2\mathbf{1}_{(t_i,t_{i+1})}(t), ~ (t-t_i)^3\mathbf{1}_{(t_i,t_{i+1})}(t) \end{equation} where as usual the indicator function $\mathbf{1}_{(t_i,t_{i+1})}(t)$ is 1 if $t \in (t_i,t_{i+1})$ and $0$ otherwise. This makes it seem like we need $4(n+1) = 4n+4$ basis functions. However, we know from linear algebra that the number of basis vectors we need is equal to the number of dimensions of the vector space. The number of adjustable coefficients for an arbitrary piecewise cubic with $n+1$ segments is indeed $4n+4$, but splines are constrained to be smooth. The spline must be continuous, which means that at each $t_i$, the value of the cubic from the left, defined on $(t_{i-1}, t_i)$, must match the value of the cubic from the right, defined on $(t_{i}, t_{i+1})$. This gives us one constraint per data point, reducing the number of adjustable coefficients to at most $3n+4$. Since the first and second derivatives are also continuous, we are down to just $n+4$ coefficients. Finally, we know that the spline function is linear outside the range of the data, i.e., on $(-\infty,t_1)$ and on ($t_n,\infty)$, lowering the number of coefficients to $n$. There are no more constraints, so we end up needing only $n$ basis functions. And in fact, from linear algebra, any set of $n$ piecewise cubic functions which are linearly independent can be used as a basis. One common choice is \begin{eqnarray} B_1(t) & = & 1\\ B_2(t) & = & t\\ B_{i+2}(x) & = & \frac{(t-t_i)_{+}^3 - (t-t_n)_{+}^3}{t_n - t_i} - \frac{(t-t_{n-1})_{+}^3 - (t-t_n)_{+}^3}{t_n - t_{n-1}} \end{eqnarray} where $(a)_{+} = a$ if $a > 0$, and $=0$ otherwise. This rather unintuitive-looking basis has the nice property that the second and third derivatives of each $B_j$ are zero outside the interval $(t_1, t_n)$. Now that we have our basis functions, we can once again write the spline as a weighted sum of them, \begin{equation} m(t) = \sum_{j=1}^{m}{\beta_j B_j(t)} \end{equation} and put together the matrix $\mathbf{b}$ where $b_{ij} = B_j(t_i)$. We can write the spline objective function in terms of the basis functions, \begin{equation} n\mathcal{L} = (\mathbf{x}-\mathbf{b}\beta)^T (\mathbf{x}-\mathbf{b}\beta) + n\lambda \beta^T \Omega \beta \end{equation} where the matrix $\Omega$ encodes information about the curvature of the basis functions: \begin{equation} \Omega_{jk} = \int{B^{\prime\prime}_j(t) B^{\prime\prime}_k(t) dt} \end{equation} Notice that only the quadratic and cubic basis functions will make non-zero contributions to $\Omega$. With the choice of basis above, the second derivatives are non-zero on, at most, the interval $(t_1, t_n)$, so each of the integrals in $\Omega$ is going to be finite. $\Omega$ is something we (or, realistically, R) can calculate _once_, no matter what $\lambda$ is. Now we can find the smoothing spline by differentiating with respect to $\beta$: \begin{eqnarray} 0 & = & -2\mathbf{b}^T\mathbf{x} + 2 \mathbf{b}^T\mathbf{b} \hat{\beta} + 2n\lambda\Omega\hat{\beta}\\ \mathbf{b}^T\mathbf{x} & = & \left(\mathbf{b}^T \mathbf{b} + n\lambda\Omega\right)\hat{\beta}\\ \hat{\beta} & = & {\left(\mathbf{b}^T \mathbf{b} + n\lambda\Omega\right)}^{-1} \mathbf{b}^T\mathbf{x} \label{eqn:spline-coefficients-as-ridge-regression} \end{eqnarray} Notice, incidentally, that we can now show splines are linear smoothers: \begin{eqnarray} \mathbf{\EstRegFunc} & = & \mathbf{b}\hat{\beta}\\ & = & \mathbf{b}{\left(\mathbf{b}^T \mathbf{b} + n\lambda\Omega\right)}^{-1} \mathbf{b}^T\mathbf{x} \end{eqnarray} Once again, if this were ordinary linear regression, the OLS estimate of the coefficients would be $(\mathbf{t}^T\mathbf{t})^{-1} \mathbf{t}^T \mathbf{x}$. In comparison to that, we've made two changes. First, we've substituted the basis-function matrix $\mathbf{b}$ for the original matrix of predictor variables, $\mathbf{t}$ --- a change we'd have made already for a polynomial regression. Second, the "denominator" is not $\mathbf{t}^T\mathbf{t}$, or even $\mathbf{b}^T \mathbf{b}$, but $\mathbf{b}^T \mathbf{b} + n\lambda\Omega$. Since $\mathbf{t}^T\mathbf{t}$ is $n$ times the covariance matrix of the independent variables, we are taking the covariance matrix of the spline basis functions and adding some extra covariance --- how much depends on the shapes of the functions (through $\Omega$) and how much smoothing we want to do (through $\lambda$). The larger we make $\lambda$, the less the actual data matters to the fit. In addition to explaining how splines can be fit quickly (do some matrix arithmetic), this illustrates two important tricks. One, which we won't explore further here, is to turn a nonlinear regression problem into one which is linear in another set of basis functions. This is like using not just one transformation of the input variables, but a whole library of them, and letting the data decide which transformations are important. There remains the issue of selecting the basis functions, which can be quite tricky. In addition to the spline basis, most choices are various sorts of waves --- sine and cosine waves of different frequencies, various wave-forms of limited spatial extent ("wavelets"), etc. The ideal is to chose a function basis where only a few non-zero coefficients would need to be estimated, but this requires some understanding of the data... The other trick is that of stabilizing an unstable estimation problem by adding a penalty term. This reduces variance at the cost of introducing some bias. #### Effective degrees of freedom Since the effective degrees of freedom of a linear smoother is $\tr{\mathbf{w}}$, we've just shown that the a spline has \[ \tr{\left(\mathbf{b}{\left(\mathbf{b}^T \mathbf{b} + n\lambda\Omega\right)}^{-1} \mathbf{b}^T\right)} \] degrees of freedom. Holding the data (and so $\mathbf{b}$ and $\Omega$) fixed, as $\lambda$ goes up, this will go down, and vice versa. # Choosing How, or How Much, to Smooth --- Cross-validation Pretty much any smoothing method we might use has some adjustable control settings[^hyperparameters] --- the size of the moving-average window, say, or the "bandwidth" of a kernel --- which we need to set. How should we tune these control settings? [^hyperparameters]: Some people call these "hyperparameters", rather than "control settings". I think this is a misleading name, because they're not in any way parameters of the process we're studying, they're just aspects of _our_ method. One of the simplest, but also most practical, ways of picking control settings is to **make the smoother predict well on data it didn't see**. The uses to which we're putting our smoother --- interpolating to fill in gaps in the data, extrapolating beyond the range of the data, filtering to remove noise from the data --- are all, in the end, about predicting unseen values. The most natural way to make sure the smoother can predict unseen values is to _have it predict unseen values_. This brings us to the idea of **cross-validation**. There are many variants, but the key notion is this: - Take the data set, and divide the data points into a **training set** and a **testing set**. - Fit each smoother (or other model) to the training set. - Evaluate each smoother on the testing set, say by taking the squared error. - Average the scores for each model over different choices of training and testing sets. This gives a score for each smoother, which says how well it was able to predict data which it _didn't_ get to see during the fitting or training process. This is exactly what we want. ## Leave-one-out Cross-validation The most straight-forward form of this is **leave-one-out cross-validation** (LOOCV).[^adaripoff] We do the following for each smoother we're considering: [^adaripoff]: Some of the text in this section is modified from @CRS-ADAfaEPoV. - For each of the $n$ data points: + Fit each smoother using every data point _except_ data point $i$; call this function $\EstRegFunc^{(-i)}$; + Calculate the prediction at data point $i$, $\EstRegFunc^{(-i)}(t_i)$; + Calculate squared[^squarederror] prediction error $(x_i - \EstRegFunc^{(-i)}(t_i))^2$. - Average the squared prediction errors for all models over all $n$ data points, to get $n^{-1}\sum_{i=1}^{n}{(x_i - \EstRegFunc^{(-i)}(t_i))^2}$. [^squarederror]: The same idea works if you'd rather use the absolute error, a "robust" loss function, etc. This makes it sound as though, with LOOCV, we need to estimate each possible smoother $n$ times, which would be slow for a large data set. There is in fact a trick for linear smoothers which lets us fit the smoother just once. What happens when we hold back data point $i$, and then make a prediction at $t_i$? The observed response at $i$ can't contribute to the prediction, but otherwise the linear smoother should work as before, so \[ \EstRegFunc^{(-i)}(t_i)= \frac{({\mathbf{w} \mathbf{x})}_i - w_{ii} x_i}{1-w_{ii}} \] The numerator just removes the contribution to $\EstRegFunc(t_i)$ that came from $x_i$, and the denominator just re-normalizes the weights in the smoother. Now a little algebra says that \[ x_i - \EstRegFunc^{(-i)}(t_i) = \frac{x_i - \EstRegFunc(x_i)}{1-w_{ii}} \] The quantity on the left of that equation is what we want to square and average to get the leave-one-out CV score, but everything on the right can be calculated from the fit we did to the whole data. The leave-one-out CV score is therefore \begin{equation} \frac{1}{n}\sum_{i=1}^{n}{\left(\frac{x_i-\EstRegFunc(t_i)}{1-w_{ii}}\right)^2} \end{equation} Again, this formula requires _both_ leave-one-out-CV _and_ a linear smoother. If we had independent observations, you might think that LOOCV would give us $n$ independent, unbiased estimates of how well a smoother predicts new data, and so an unbiased, low-variance estimate of model performance. The problem is that even if the _data points_ are independent, $\EstRegFunc^{(-i)}$ and $\EstRegFunc^{(-j)}$ are based on $n-2$ shared data points, so the leave-one-out errors for each data point are pretty correlated with each other. (Can you work out the covariance between $x_i - \EstRegFunc^{(-i)}(t_i)$ and $x_j - \EstRegFunc^{(-j)}(t_j)$?) One consequence of this is that LOOCV tends to over-fit --- to pick models with more degrees of freedom than are strictly necessary. We sometimes therefore prefer to use _many_ data points in each testing set, so the predictions are less correlated with each other. The commonest approach is **$k$-fold CV**, where we divide all the data points into $k$ equal groups, or "folds", and each fold gets to be the testing set in turn. (Typically $k=5$ or $10$.) This is much less prone to over-fitting than leave-one-out. ## Special issues on cross-validation for dependent data These forms of cross-validation work very, very well for independent data, to the point where they're pretty much the standard for any situation involving statistical modeling of independent data. But what about dependent data? If each data point tends to be like its neighbors, you might well worry that a smoother which _just_ copies the neighbors will have an unfair advantage in cross-validation. At the very least, this would seem to lead to more noise in our estimate of how well each smoother predicts. A natural way out of this is to do leave-one-out cross-validation, but to leave out a buffer of size (say) $h$ around each point [@Burman-Chow-Nolan-cross-validation-for-dependent]. This buffer is not part of the training set, but we also don't use it in evaluating the prediction of the smoother. (We _do_ use values in the buffer when we _calculate_ the prediction at the held-out point, however.) If remote observations are nearly independent of each other, and the buffer is big enough, then the testing set is nearly independent of the training set, and everything is well-behaved. You can combine this idea with $k$-fold cross-validation, where you omit a buffer around each block [@Racine-consistent-cv-for-dependent-data]. The obvious question is "how big a buffer?", to which the obvious answer is "big enough to make the training and test points nearly independent". Since we usually don't know that, we typically rely on rules of thumb. @Burman-Chow-Nolan-cross-validation-for-dependent suggest using a constant fraction of $n$, say $h=0.5 n$. @Racine-consistent-cv-for-dependent-data suggests that this rule is prone to over-fitting, and that a better one would be $h = n^{0.25}$. I don't know of a really great, automatic resolution of this issue. In practice, it can often be easiest to just live with the extra imprecision and over-fitting that come from ignoring dependence, but we'll come back to this later when we look in more detail at time-series forecasting. ## Smoothing ratios and proportions Before moving on from smoothers, I should mention one last special consideration. When our data consists of a ratio or proportion, say $X(t) = Y(t)/Z(t)$, we _could_ smooth it the same way we'd smooth any other numerical data. However, if the denominator is small, any noise in it can have a massive effect on the ratio. It can make a lot more sense to _separately_ smooth $Y$ and $Z$, and then take the ratio of the smoothed values. @Kafadar-smoothing-geographical, from whom I learned this point, explores how best to do this with splines and related smoothers. # Detrended values and residuals After we have chosen how to smooth the data, and gotten an estimate of the trend, we also have an estimate of the fluctuations. Since we can always write \[ X(r,t) = \TrueRegFunc(r,t) + \TrueNoise(r,t) \] where $\TrueRegFunc$ is a deterministic (but unknown) function and $\TrueNoise$ is a stochastic process, the natural estimate of the latter is \[ \EstNoise(r,t) \equiv X(r,t) - \EstRegFunc(r,t) \] These natural estimates of the noise are the **residuals** of the smoother. When we use a linear smoother, we can find the residuals in matrix form: \begin{eqnarray} \mathbf{\EstNoise} & = & \mathbf{x} - \mathbf{\EstRegFunc}\\ & = & \mathbf{x} - \mathbf{w}\mathbf{x}\\ & = & (\mathbf{I} - \mathbf{w})\mathbf{x} \end{eqnarray} Thus the residuals, like the fitted values, also linear in $\mathbf{x}$. We can therefore work out properties of the residuals much as we did properties of the fitted values. The expected value of the residual vector is \begin{eqnarray} \Expect{\mathbf{\EstNoise}} & = & \Expect{(\mathbf{I}-\mathbf{w})\mathbf{X}}\\ & = & (\mathbf{I}-\mathbf{w})\mathbf{\TrueRegFunc} \end{eqnarray} This will be zero --- i.e., match the true expectation of the noise --- if and only if $\TrueRegFunc$ is an eigenvalue of $\mathbf{I}-\mathbf{w}$ with eigenvalue 0. --- It should not be hard to convince yourself that if $\mathbf{v}$ is an eigenvector of $\mathbf{w}$ with eigenvalue $\lambda$, then $\mathbf{v}$ is also an eigenvector of $\mathbf{I}-\mathbf{w}$, with eigenvalue $1-\lambda$. It should also not be hard to check that if $\mathbf{w}$ is symmetric and idempotent, then so is $\mathbf{I} - \mathbf{w}$. The variance of the residuals is \[ \Var{\mathbf{\EstNoise}} = (\mathbf{I}-\mathbf{w}) \Var{\mathbf{\epsilon}} (\mathbf{I}-\mathbf{w})^T \] Under the special assumption that the true noise has constant variance and is uncorrelated, $\Var{\mathbf{\epsilon}} = \sigma^2 \mathbf{I}$, and this simplifies to \[ \sigma^2 (\mathbf{I}-\mathbf{w})(\mathbf{I}-\mathbf{w})^T \] Notice that in general this is _not_ $\sigma^2 \mathbf{I}$. That is, the residuals generally have somewhat different variances, and correlation structure, than the true noise. Homework 2 will, among other things, explore some consequences of this. # Linear regression as a linear smoother > This section can be skipped, if you're willing to take the fact that linear regression is really, or also, a linear smoother on trust. We will be coming back to some of this stuff later, when we look at linear predictors more generally. ## Simple linear regression as a linear smoother Think of the simple linear regression model: \[ m(t) = \beta_0 + \beta_1 t \] We usually estimate this by minimizing the mean squared error (MSE) --- what's called the **ordinary least sqaures** or OLS estimate: \begin{eqnarray} (\hat{\beta}_0, \hat{\beta}_1) = \argmin_{b_0, b_1}{\frac{1}{n}\sum_{i=1}^{n}{(x_i - (b_0 + b_1 t_i))^2}} \end{eqnarray} Take the derivatives with respect to $b_0$ and $b_1$, and set them to zero. \begin{eqnarray} \frac{1}{n}\sum_{i=1}^{n}{2(x_i - (\hat{\beta}_0 + \hat{\beta}_1 t_i))(-1)} & = & 0\\ \frac{1}{n}\sum_{i=1}^{n}{2(x_i - (\hat{\beta}_0 + \hat{\beta}_1 t_i))(-t_i)} & = & 0 \end{eqnarray} We can drop the over-all factor of $-2$ in what follows. First, solve for $\hat{\beta}_0$ in terms of $\hat{\beta}_1$: \begin{eqnarray} \overline{x} - \hat{\beta}_0 - \hat{\beta}_1 \overline{t} & = & 0\\ \overline{x} - \hat{\beta}_1 \overline{t} & = & \hat{\beta}_0 \end{eqnarray} Now solve for $\hat{\beta}_1$: \begin{eqnarray} \overline{xt} - \hat{\beta}_0 \overline{t} - \hat{\beta}_1 \overline{t^2} & = & 0\\ \overline{xt} - \overline{x}\overline{t} + \hat{\beta}_1 (\overline{t})^2 - \hat{\beta}_1 \overline{t^2} & = & 0\\ \frac{\overline{xt} - \overline{x}\overline{t}}{\overline{t^2} - \overline{t}^2} & = & \hat{\beta}_1 \end{eqnarray} It's convenient to abbreviate the denominator as $s_t^2$. Putting this together, \[ \EstRegFunc(t) = \overline{x} + \frac{t - \overline{t}}{s_t^2}(\overline{xt} - \overline{x}\overline{t}) \] This is linear in $x$, so we know it can be put in the linear-smoother form. In fact, \[ \EstRegFunc(t) = \sum_{j=1}^{n}{\left(\frac{1}{n} + \frac{(t-\overline{t})(t_j - \overline{t})}{n s_t^2}\right) x_j} \] This particular form of weights ensures both that the estimated curve is a straight line, and that it comes as close to the data as possible, on average, among all straight lines. Beyond that, however, the weights are really weird. The first part is just the sample mean, but the variable part is $\propto (t-\overline{t})(t_j-\overline{t})$. This means we give more weight to data points which are far from the center of the data, no matter what. We also ramp up all the weights as the point $t$ where we're making a prediction moves away from the center of the data. ## Linear regression with multiple predictors Suppose we want to make our prediction a linear combination of multiple variables: \[ m(t,z, \ldots u) = \beta_0 + \beta_ 1 + \beta_2 z + \ldots \beta_p u \] We fit this by least squares again. It simplifies the math immensely if we define two matrices, \[ \mathbf{t} = \left[ \begin{array}{ccccc} 1 & t_1 & z_1 & \ldots & u_1\\ 1 & t_2 & z_2 & \ldots & u_2\\ \vdots & \vdots & \vdots & \ldots & \vdots\\ 1 & t_n & z_n & \ldots & u_n \end{array} \right] \] and \[ \mathbf{b} = \left[ \begin{array}{c} b_0 \\ b_1 \ldots b_p \end{array}\right] \] Then the fitted values are \[ \mathbf{m} = \mathbf{t}\mathbf{b} \] and the mean-squared error is \[ \frac{1}{n}\sum_{i=1}^{n}{(x_i - (b_0 + b_1 t_i + b_2 z_i + \ldots + b_p u_i))^2} = \frac{1}{n} (\mathbf{x} - \mathbf{t}\mathbf{b})^T (\mathbf{x} - \mathbf{t}\mathbf{b}) \] Taking the derivative and setting it to zero \[ \frac{2}{n}\mathbf{t}^T (\mathbf{x} - \mathbf{t}\mathbf{\widehat{\beta}}) = 0 \] Solving, \[ \mathbf{\widehat{\beta}} = (\mathbf{t}^T\mathbf{t})^{-1} \mathbf{t}^T \mathbf{x} \] The fitted valued come from \[ \mathbf{\EstRegFunc} = \mathbf{t} (\mathbf{t}^T\mathbf{t})^{-1} \mathbf{t}^T \mathbf{x} \] so the weight matrix in this case is \[ \mathbf{w} = \mathbf{t} (\mathbf{t}^T\mathbf{t})^{-1} \mathbf{t}^T \] You can check that this matrix is symmetric, and that it is idempotent. You can also check that when we just deal with one predictor, these matrix formulas are equivalent to the ones for simple linear regression above. These formula handles _all_ of the coefficients at once. You can check, by explicitly taking the partial derivative just with respect to $b_0$ and setting that to 0, that it is still the case that \[ \hat{\beta}_0 = \overline{y} - (\hat{\beta}_1 \overline{t} + \hat{\beta}_2 \overline{z} + \ldots + \hat{\beta}_p \overline{u}) \] This tells us that the individual weights in $\mathbf{w}$ still have the same form --- they're a constant, plus terms that give more weight to points far from the center of the data. ## Transformations of one predictor; polynomial regression Nothing in the logic for the previous section requires that the extra predictors be new variables. If we want to do polynomial regression, \[ m(t) = \beta_0 + \beta_1 t + \beta_2 t^2 + \ldots \beta_p t^p \] then everything carries through exactly as before. The matrix $\mathbf{t}$ would just have columns containing the constant, $t_i$, $t_i^2$ and so on. What's important is that we're just taking an additive combination of the functions of $t$ on the right-hand side of the model formula, and that all of the functions on the RHS are linearly independent of each other. (If there was linear dependence, then the matrix $\mathbf{t}^T \mathbf{t}$ wouldn't be invertible.) # Some linear algebra reminders: eigen-stuff > This section can be skipped if you are comfortable working with the eigenvalues and eigenvectors of matrices. Even if you are, it won't hurt to at least skim it. A (non-zero) vector or $n\times 1$ matrix $\mathbf{v}$ is an **eigenvector** of the $n\times n$ matrix $\mathbf{w}$, with **eigenvalue** $\lambda$, when \[ \mathbf{w}\mathbf{v} = \lambda\mathbf{v} \] We can re-write this as \[ (\mathbf{w} - \lambda\mathbf{I})\mathbf{v} = 0 \] and demanding $\mathbf{v} \neq 0$ tells us that the $\lambda$ which make this true have to obey the **characteristic polynomial** or **eigenpolynomial** of the matrix \[ \det{(\mathbf{w} - \lambda\mathbf{I})} = 0 \] In R, the command to find the eigenvalues and eigenvectors of a matrix is `eigen`: ```{r} theta <- runif(n=2) # Why "mark"? mark <- matrix(c(theta[1], 1-theta[1], 1-theta[2], theta[2]), byrow=TRUE, nrow=2, ncol=2) mark eigen(mark) ``` This returns a list with two components, one of them a vector of eigenvalues, the other a matrix whose columns are the eigenvectors. ## General properties of eigenvalues and eigenvectors There are a number of salient properties of eigenvalues and eigenvectors. 1. There are at most $n$ distinct eigenvalues of $\mathbf{w}$, which might be complex[^neigenvalues]. When there are fewer than $n$ distinct eigenvalues, we nonetheless write the eigenvalues as $\lambda_1, \lambda_2, \ldots \lambda_n$, keeping track of how often each one is a root of the characteristic polynomial. [^neigenvalues]: $\det{(\mathbf{w} - \lambda\mathbf{I})} = 0$ is called the characteristic _polynomial_ because working out the determinant gives us an $n^{\mathrm{th}}$ order polynomial in $\lambda$. This has at most $n$ distinct roots, some of which might be complex; the same $\lambda$ can appear multiple times as a root. 2. The determinant of a matrix is the product of its eigenvalues, "counted with multiplicity": \[ \det{\mathbf{w}} = \prod_{j=1}^{n}{\lambda_j} \] The trace of a matrix is the sum of its eigenvalues: \[ \tr{\mathbf{w}} = \sum_{j=1}^{n}{\lambda_j} \] 3. For any matrix $\mathbf{w}$, and any power $k > 0$, $\mathbf{v}$ is an eigenvector of $\mathbf{w}$ with eigenvalue $\lambda$ if and only if $\mathbf{v}$ is an eigenvector of $\mathbf{w}^k$ with eigenvalue $\lambda^k$. If $\mathbf{w}$ is invertible, then $\mathbf{v}$ is an eigenvector of $\mathbf{w}^{-1}$ with eigenvalue $1/\lambda$. 4. If $\mathbf{v}$ is an eigenvector with eigenvalue $\lambda$, then for any scalar $a \neq 0$, $a\mathbf{v}$ is also an eigenvector, with the same eigenvalue[^scaleeigenvectors]. This means that the eigenvectors, unlike the eigenvalues, aren't uniquely defined. Instead, we usually _choose_ eigenvectors to have nice properties, e.g., length 1. Also, $-\mathbf{v}$ is always another eigenvector, so two calculations of eigenvectors might end up with different signs. [^scaleeigenvectors]: By the ordinary rules for matrix multiplication, $\mathbf{w}(a\mathbf{v}) = a\mathbf{w}\mathbf{v} = a\lambda \mathbf{v}$, because $\mathbf{v}$ is an eigenvector. 5. If $\mathbf{v}$ and $\mathbf{u}$ are eigenvectors with the same eigenvalue $\lambda$, then so is $a\mathbf{v}+b\mathbf{u}$, for any scalars $a, b$. (Check this yourself!) 6. For fixed $\lambda$, the eigenvectors solve the equation $(\mathbf{w} - \lambda\mathbf{I})\mathbf{v} = 0$. This matrix equation is equivalent to a system of $n$ linear equations in the $n$ unknowns $(v_1, v_2, \ldots v_n)$. If there are no redundant (linearly dependent) equations in this system, then the solution is unique, up to the over-all scaling factor we just saw in (4). But if there are redundant equations, we can have a whole space of solutions, the **eigenspace** of $\lambda$. The number of redundant equations gives the dimension of this eigenspace, a.k.a. the **geometric multiplicity** of $\lambda$. When there are multiple, linearly-independent eigenvectors with eigenvalue $\lambda$, we can _choose_ orthogonal eigenvectors to represent the whole eigenspace. 7. If $\mathbf{w}$ is symmetric, and $\mathbf{v}$ and $\mathbf{u}$ are eigenvectors with distinct eigenvalues $\lambda_1 \neq \lambda_2$, then $\mathbf{v}^T\mathbf{u} = 0$, and we say the vectors are orthogonal[^orthgonaleigenvectors]. _Remark_: The orthogonality of eigenvectors actually holds for many asymmetric matrices as well, so long as they are **normal**, meaning $\mathbf{w}^T\mathbf{w} = \mathbf{w}\mathbf{w}^T$. The proof is similar, but more involved. I am afraid I mis-spoke in lecture when I asserted that it is always true for all matrices[^counterexample]. [^orthgonaleigenvectors]: Since $\mathbf{w}$ is symmetric, $\mathbf{w} = \mathbf{w}^T$, we have $\mathbf{u}^T \mathbf{w}\mathbf{v} = \mathbf{v}^T \mathbf{w}\mathbf{u}$. Because these are eigenvectors, this means \begin{eqnarray} \mathbf{u}^T \lambda_1 \mathbf{v} & = & \mathbf{v}^T \lambda_2 \mathbf{u}\\ \lambda_1 \mathbf{u}^T\mathbf{v} & = & \lambda_2 \mathbf{v}^T \mathbf{u} \end{eqnarray} but since $\mathbf{v}^T\mathbf{u} = \mathbf{u}^T\mathbf{v}$ (it's a $1\times 1$ matrix, which is always symmetric), we get \[ (\lambda_1 - \lambda_2)( \mathbf{u}^T\mathbf{v}) = 0 \] and we assumed $\lambda_1 \neq \lambda_2$. ```{r, include=FALSE} w <- matrix(c(2,2,0,3), nrow=2) eigen(w)$values e <- eigen(w)$vectors e[,1] %*% e[,2] ``` [^counterexample]: Here's a counter-example. The matrix $\left[ \begin{array}{cc} 2 & 0 \\ 2 & 3 \end{array} \right]$, has eigenvalues are $`r eigen(w)$values[1]`$ and $`r eigen(w)$values[2]`$, and the corresponding eigenvectors are $\left[ `r e[,1]` \right]$ and $\left[ `r e[,2]` \right]$, which aren't orthogonal. 8. What is _generally_ true, of most matrices, is that the eigenvectors span the whole $n$-dimensional space. (There are "defective" matrices where this is not true.) In those cases, we can stack the eigenvectors column-wise into a matrix, say $\mathbf{e}$, and have $\mathbf{w}\mathbf{e} =\mathbf{e} \mathbf{\lambda}$, where $\mathbf{\lambda}$ is the diagonal matrix of the eigenvalues. But, because the eigenvectors are (by assumption) linearly independent, $\mathbf{e}^{-1}$ exists, and we have the **eigendecomposition** \[ \mathbf{w} = \mathbf{e}\mathbf{\lambda}\mathbf{e}^{-1} \] ## Eigen-properties of some special types of matrices ### Symmetric matrices If $\mathbf{w} = \mathbf{w}^T$, then all the eigenvalues are real. Because eigenvectors with different eigenvalues are distinct in this case, and we can always normalize all the eigenvectors, we have $\mathbf{e}^{-1} = \mathbf{e}^T$, and $\mathbf{w} = \mathbf{e}\mathbf{\lambda}\mathbf{e}^T$. ### "Stochastic" matrices A matrix is called **stochastic** if $w_{ij} \geq 0$ for all $i, j$, and if $\sum_{j=1}^{n}{w_{ij}} = 1$ for all $i$. Some very important properties follow from this: 1. The largest eigenvalue is exactly 1. 2. The eigenvectors corresponding to the eigenvalue 1 have only non-negative entries. 3. All other eigenvalues have magnitude $\leq 1$. These are all consequences of a result in linear algebra called the Frobenius-Perron theorem, which we will take on trust. ### Projection / idempotent matrices A matrix is called **idempotent**[^idempotent] when $\mathbf{w}^2 = \mathbf{w} \mathbf{w} = \mathbf{w}$. It follows[^idempotenteigenvalues] that the only possible eigenvalues are 0 and 1. Projecting on to a linear subspace (line, plane, hyper-plane, etc.) always corresponds to an idempotent matrix. (Points already in the space aren't affected by the projection, so projecting twice is the same as projecting once.) The dimension of the space corresponds to the number of eigenvalues equal to 1. [^idempotent]: From Latin _idem_, "the same", and _potens_, "power", hence something which is the same as its powers. [^idempotenteigenvalues]: If $\mathbf{v}$ is an eigenvector of $\mathbf{w}$ with eigenvalue $\lambda$, then clearly $\mathbf{w}^2 \mathbf{v} = \lambda^2\mathbf{v}$, so $\lambda^2$ is an eigenvalue of $\mathbf{w}$. But, since we're assuming $\mathbf{w}$ is idempotent, $\mathbf{w}^2 = \mathbf{w}$, so all eigenvalues must obey the equation $\lambda^2=\lambda$, and the only two solutions are 0 and 1. # References