---
title: Trends and Smoothing II
author: 36-467/36-667
date: 4 September 2018
output: slidy_presentation
---

\[
\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}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\dof}{DoF}
\DeclareMathOperator{\det}{det}
\newcommand{\TrueNoise}{\epsilon}
\newcommand{\EstNoise}{\widehat{\TrueNoise}}
\]


## In our last episode...

> - Data $X(t) = \TrueRegFunc(t) + \TrueNoise(t)$
> - $\TrueRegFunc$ deterministic (=**trend**), $\TrueNoise$ stochastic and mean-zero (=**fluctuations**)
> - Wanted: estimates of $\TrueRegFunc$ and/or $\TrueNoise$ from _one_ data set
> - Hope: $\TrueRegFunc$ is a smooth function $\Rightarrow$ average nearby $X$'s
> - **Linear smoother**: $\EstRegFunc(t) = \sum_{j=1}^{n}{w(t, t_j) x_j}$
> - **Fitted values** on the data $\mathbf{\EstRegFunc} = \mathbf{w}\mathbf{x}$
> - $\mathbf{w}$ is the source of all knowledge

## Expectation of the fitted values

\begin{eqnarray}
\Expect{\mathbf{\EstRegFunc}} & = & \Expect{\mathbf{w}\mathbf{X}}\\
& = & \mathbf{w}\Expect{\mathbf{X}}\\
& = & \mathbf{w} \mathbf{\mu}
\end{eqnarray}

**Unbiased** estimate $\Leftrightarrow \mathbf{w} \mathbf{\mu} = \mathbf{\mu}$

## Expanding in eigenvectors

> - Generally, $\mathbf{w}$ has $n$ linearly-independent eigenvectors $\mathbf{e}_1, \ldots \mathbf{e}_n$, with eigenvalues $\lambda_1, \ldots \lambda_n$
> - So $\mathbf{x} = \sum_{j=1}^{n}{c_j \mathbf{e}_j}$
> - So $\mathbf{w}\mathbf{x} = \mathbf{w}\sum_{j=1}^{n}{c_j \mathbf{e}_j} = \sum_{j=1}^{n}{c_j \lambda_j \mathbf{e}_j}$
> - Components of the data which match large-$\lambda$ eigenvectors are enhanced
> - Components of the data which match small-$\lambda$ eigenvectors are shrunk

## A little example

```{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
```

## A little example

```{r}
w
```

## A little example

```{r}
eigen(w)$values
eigen(w)$vectors[,1]
```

## A little example

```{r, echo=FALSE}
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"))
```

## Variance of the fitted values

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


IF $\Var{\mathbf{\TrueNoise}} = \sigma^2 \mathbf{I}$, THEN
$\Var{\mathbf{\EstRegFunc}} = \sigma^2\mathbf{w}\mathbf{w}^T$

## How much do the fitted values respond to the data?

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

IF $\Var{\mathbf{\TrueNoise}} = \sigma^2 \mathbf{I}$, THEN
this $= \sigma^2\tr{\mathbf{w}} = \sigma^2 \text{(sum of eigenvalues)}$

$\tr{\mathbf{w}} =$ (effective) **degrees of freedom**

## Data = trend + fluctuation

- $X(t) = \TrueRegFunc(t) + \TrueNoise(t)$
- $\Rightarrow$ $\TrueNoise(t) = X(t) - \TrueRegFunc(t)$
- $\Rightarrow$ $\EstNoise(t) \equiv X(t) - \EstRegFunc(t) =$ **residuals**

\begin{eqnarray}
\mathbf{\EstNoise} & = & \mathbf{x} - \mathbf{\EstRegFunc}\\
& = & \mathbf{x} - \mathbf{w}\mathbf{x}\\
& = & (\mathbf{I} - \mathbf{w})\mathbf{x}
\end{eqnarray}

_Convince yourself_: $\mathbf{I}-\mathbf{w}$ has same eigenvectors as
$\mathbf{w}$, but eigenvalues $1-\lambda$

## Expected residuals

\begin{eqnarray}
\Expect{\mathbf{\EstNoise}} & = & \Expect{(\mathbf{I}-\mathbf{w})\mathbf{X}}\\
& = & (\mathbf{I}-\mathbf{w})\mathbf{\TrueRegFunc}
\end{eqnarray}

Biased trend estimate $\Leftrightarrow$ biased fluctuation estimate

## Break for the in-class exercise

- $X(t) = \TrueRegFunc(t) + \TrueNoise(t)$, and $\Var{\TrueNoise(t)} = \sigma^2$,
$\Cov{\TrueNoise(t_1), \TrueNoise(t_2)} = 0$
- Set $\EstRegFunc(t) = \frac{1}{3}\sum_{s=t-1}^{t+1}{X(s)}$
    + Ignore the ends of the data where we don't have neighbors on both sides
- What is $\Cov{\EstRegFunc(t), \EstRegFunc(t+1)}$?
- What is $\Cov{\EstRegFunc(t), \EstRegFunc(t+2)}$?
- What is $\Cov{\EstRegFunc(t), \EstRegFunc(t+3)}$?
- Why aren't all of these 0?

## Variance and covariance of the residuals

\[
\Var{\mathbf{\EstNoise}} =  (\mathbf{I}-\mathbf{w}) \Var{\mathbf{\epsilon}} (\mathbf{I}-\mathbf{w})^T
\]

IF $\Var{\mathbf{\epsilon}} = \sigma^2 \mathbf{I}$, THEN this $= \sigma^2 (\mathbf{I}-\mathbf{w})(\mathbf{I}-\mathbf{w})^T$

**NB**: Correlations from off-diagonal entries in $\mathbf{w}$




## Splines

\[
\EstRegFunc = \argmin_{m}{\frac{1}{n}\sum_{i=1}^{n}{(x_i - m(t_i))^2} + \lambda\int{(m^{\prime\prime}(t))^2 dt}}
\]

> - This $\lambda$ not an eigenvalue (sorry)
> - Fit the data points vs. over-all curvature
> - Minimization is over _all_ functions
> - Solution is always a piecewise cubic polynomial, but continuous, with continuous 1st and 2nd derivatives
> - $\lambda \rightarrow 0$ $\Rightarrow$ Straight lines between data points
> - $\lambda \rightarrow \infty$ $\Rightarrow$ Global linear fit
> - $\downarrow$ degrees of freedom as $\uparrow \lambda$

## How do we pick $\lambda$?

> - Want trend to predict not-yet-seen stuff (interpolate, extrapolate, filter)
> - A good $\lambda$ predicts new stuff well
> - Hold out part of the data and try to predict that from the rest

## Leave-one-out cross-validation (LOOCV)

- For each of the $n$ data points:
    + Fit using every data point _except_ $i$, get $\EstRegFunc^{(-i)}$;
    + Find $\EstRegFunc^{(-i)}(t_i)$;
    + Find $(x_i - \EstRegFunc^{(-i)}(t_i))^2$.
- Average over all data points, $n^{-1}\sum_{i=1}^{n}{(x_i - \EstRegFunc^{(-i)}(t_i))^2}$

- Low LOOCV $\Leftrightarrow$ good ability to predict new data
- This is what `smooth.spline` does automatically

## Leave-one-out cross-validation (LOOCV)

Don't have to re-fit linear smoothers $n$ times

\begin{eqnarray}
\EstRegFunc^{(-i)}(t_i) &= & \frac{({\mathbf{w} \mathbf{x})}_i - w_{ii} x_i}{1-w_{ii}}\\
x_i - \EstRegFunc^{(-i)}(t_i) & = & \frac{x_i - \EstRegFunc(t_i)}{1-w_{ii}}\\
LOOCV & = & \frac{1}{n}\sum_{i=1}^{n}{\left(\frac{x_i-\EstRegFunc(t_i)}{1-w_{ii}}\right)^2}
\end{eqnarray}

## Many variants

- $h$-block CV: omit a buffer of radius $h$ around the hold-out point from the training set
- $k$- or $v$-fold CV: divide data into $k$ equal-sized "folds", try to predict each fold using the rest of the data
- $hv$-block CV: $v$-fold with a buffer
- etc., et.

## The moral

- _Never_ care about how good the in-sample fit is ($R^2$, $R^2_{adj}$, etc.)
- _Always_ care about ability to predict new data



## Summing up

> - If the trend is smooth, we can estimate it by smoothing
> - Every smoother is biased towards some patterns and against others
> - Properties of the fitted values come from the weights
> - Fluctuations are residuals after removing a trend
> - De-trending can create correlations
> - We decide how to smooth by cross-validation
