2.1 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¹. This is, however, not available in non-experimental contexts.

2.2 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.

3.1 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.)

3.2 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)\).

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

eigen(w)$values

##  [1]  1.00000000  0.96261129  0.85490143  0.68968376  0.48651845
##  [6] -0.31012390  0.26920019 -0.23729622 -0.11137134  0.06254301

eigen(w)$vectors[,1]

##  [1] 0.3162278 0.3162278 0.3162278 0.3162278 0.3162278 0.3162278 0.3162278
##  [8] 0.3162278 0.3162278 0.3162278

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"))

3.3 Splines: Direct control over smoothness

I have teased you several times by now by mentioning splines, or more formally smoothing 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.

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.)

library(mgcv)
my.fit <- gam(x ~ s(t), data=df)

gam(x ~ s(r,t), data=df)

gam(x ~ te(r, t), data=df)

gam(x ~ s(r) + s(t), data=df)

4.1 Leave-one-out Cross-validation

The most straight-forward form of this is leave-one-out cross-validation (LOOCV).⁷ We do the following for each smoother we’re considering:

• 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⁸ 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}\).

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.

4.2 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, and Nolan 1994). 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 2000).

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, and Nolan (1994) suggest using a constant fraction of \(n\), say \(h=0.5 n\). Racine (2000) 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.

4.3 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 (1996), from whom I learned this point, explores how best to do this with splines and related smoothers.

6.1 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.

6.2 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.

6.3 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.)

7.1 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⁹. 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.

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¹⁰. 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.

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¹¹. 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¹².

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} \]

7.2 Eigen-properties of some special types of matrices