---
title: Trends and Smoothing I
date: 8 September 2020 (Lecture 3)
author: "36-467/667, Fall 2020"
output: slidy_presentation
bibliography: locusts.bib
---

\[
\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}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\dof}{DoF}
\]

## The mathematical setting

- Treat the data as the realization of a **random process** or **stochastic process** = random variables spread over time and space
- At point $r$ and time $t$, we see $X(r,t)$
    + Just $X(r)$ if only spatial, just $X(t)$ if only temporal
- Generally, $X(r,t)$ and $X(q,s)$ are statistically dependent
- Can also think of the process as a _random function_ of the coordinates


## Trend = central tendency of the random process

\[
\TrueRegFunc(r,t) \equiv \Expect{X(r,t)}
\]

- This is what we'll mean by the **trend**
- $\TrueRegFunc$ is a _deterministic_ function of the coordinates; not random



## Trend + Fluctuations

Random process = trend + fluctuations

\begin{eqnarray}
X(r,t) & = & \TrueRegFunc(r,t) + \TrueNoise(r,t)\\
\Expect{\TrueNoise(r,t)} & = & 0
\end{eqnarray}

(Sometimes called "signal plus noise" representation)

$\TrueNoise =$ another random process

- Can't, _in general_, assume more about $\TrueNoise$:
    + It might be correlated
    + It might not have constant variance
    + It might have different distributions everywhere
    + It might not be Gaussian
    + etc., etc.
- But we _can_ assume $\Expect{\TrueNoise(r,t)} = 0$



## Learning the Trend

- We want to learn $\TrueRegFunc$
- Either we have many realizations of the process or we just have one
- Things are a lot easier if we have many realizations

## Averaging independent realizations gives the trend

- We run the experiment many times, getting independent realizations $X^{(1)}, X^{(2)}, \ldots X^{(m)}$
- Extracting the trend is simple, because the law of large numbers applies:
\[
\frac{1}{m}\sum_{i=1}^{m}{X^{(i)}(r,t)} \rightarrow \TrueRegFunc(r,t)
\]
- Used a lot in neuroscience, physics, chemistry...
- Harder to use when we're not doing controlled experiments


## Trends from a _single_ realization

- We only have one history of cherry-blossoms in Kyoto
    + Or US economy post-WWII or shale deposits under Pennsylvania, etc.
- To get a trend, we'll need to assume _something_
- Assume a whole theory
- _or_ Assume a little about $\TrueRegFunc$

## Assume a whole theory {.smaller}

- Sometimes we actually have good theories that tell us what the trend should be
\[
\TrueRegFunc(r, t) = f(r,t;\theta)
\]
with $f$ completely specified except for some constants $\theta$
- Then all we have to do is estimate $\theta$
- Usual approach is least squares:
\[
\hat{\theta} = \argmin_{\theta}{\left(\frac{1}{n}\sum_{i=1}^{n}{(X(r_i, t_i) - f(r_i, t_i;\theta))^2} \right)}
\]
    + Issues with correlations in $\TrueNoise$ that we'll come back to later in the course
- Situations where this applies:
    + Physics, astronomy
        * This is where the method of least squares came from [@Farebrother-on-linear-relationships]
    + Some other natural sciences: chemistry, geology, climatology...
    + Parts of biology: population genetics, ecology, epidemiology
    + Social sciences: demography, parts of economics (maybe)


## Assume a whole theory

- Pros:
    + We know how to estimate parameters!
    + We can interpret estimated parameters using the theory
    + Statistically efficient if the model is right
    + Good for extrapolating if the model is right
- Cons:
    + Doesn't check the theory
    + Systematically misleading if the model is wrong
    + Only works if the theory is "tight" enough to give us a function up to unknown constants; most of science isn't like that



## Assume something about $\TrueRegFunc$

> - Assume $\TrueRegFunc$ is a **smooth** function
> - Concretely: assume $\TrueRegFunc$ changes fairly slowly, $\TrueRegFunc(t-h) \approx \TrueRegFunc(t) \approx \TrueRegFunc(t+h)$
> - We know that $X(t)$ is a noisy observation of $\TrueRegFunc(t)$
> - With this assumption, $X(t-h)$ and $X(t+h)$ are _also_ noisy observations of $\TrueRegFunc(t)$ (plus a little bit more error)
> - $\therefore$ If $\TrueRegFunc$ is smooth, we can estimate it by averaging nearby observations

## Our first smoothing estimator

\[
\EstRegFunc(t_i) = \frac{1}{3}\sum_{k=-1}^{k=+1}{X(t_{i+k})}
\]

- Why do this?

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

- Bias term:
    + Big if $\TrueRegFunc$ changes rapidly and observations are far apart
    + Small if $\TrueRegFunc$ changes slowly and observations are close
- Noise: mean 0, usually smaller than the noise variance in any one observation

## Linear smoothers

- Data: $x(t_1), \ldots x(t_n)$; abbreviate as $x_1, \ldots x_n$.
- A **linear smoother** is an estimator that looks like this:
\[
\EstRegFunc(t) = \sum_{j=1}^{n}{w(t, t_j) x_j}
\]
- Notice:
    + Linear in the observations
    + Can be _very_ nonlinear in $t$
    + Same idea with space but more notation
- Different linear smoothers $=$ different choices of weights $w$


## Examples of linear smoothers {.smaller}

- _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.
    + Also one-sided moving averages (usually only from the past)
- _Weighted moving averages_: $w(t, t_j) = g(|t-t_j|/\tau)$, for some fixed function $g$
    + Often $g$ is the (negative) exponential function, so **exponentially-weighted moving averages**.
- _Kernels_: like a weighted moving average, but weights are **normalized** (non-negative, always sum to 1), like picking a pdf $K$ and a **bandwidth** $h$ and setting
\[
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_: Linear smoothers but with very weird weights
- _Splines_: Smooth piecewise polynomials that actually do what you hoped polynomial regression would


## Properties of linear smoothers

- **Fitted values**: at each $t_i$ where we got data, we have $\EstRegFunc(t_i)$, for short $\EstRegFunc_i$
- Group the fitted values in to a vector ($n\times 1$ matrix):
\[
\mathbf{\EstRegFunc} = \left[\begin{array}{c} \EstRegFunc(t_1) \\ \EstRegFunc(t_2) \\ \vdots \\ \EstRegFunc(t_n) \end{array} \right]
\]
- Group all the observations of $x$ into a vector $\mathbf{x}$
- The fitted values are the observations _times_ a matrix:
\[
\mathbf{\EstRegFunc} = \mathbf{w} \mathbf{x}
\]
where $w_{ij} = w(t_i, t_j)$.
- The matrix $\mathbf{w}$ is the **smoothing**, **influence** or **hat** matrix
- Important: the fitted values $\mathbf{\EstRegFunc}$ are _random_

## Expectation of the fitted values; bias

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

- **Bias** = expected estimate minus true value
    + **Unbiased** = bias of 0
- When is a linear smoother unbiased?
\begin{equation}
\text{unbiased trend estimates} \Leftrightarrow \mathbf{\TrueRegFunc} = \mathbf{w}\mathbf{\TrueRegFunc}
\end{equation}
- A linear smoother is unbiased if, and only if, $\mathbf{\TrueRegFunc}$ is an eigenvector of $\mathbf{w}$ with eigenvalue 1

## Eigenvectors and eigenvalues (reminders)

- A vector $\mathbf{v}$ is an **eigenvector** of the matrix $\mathbf{w}$, with **eigenvalue** $\lambda$, when
\[
\mathbf{w}\mathbf{v} = \lambda \mathbf{v}
\]
- An $n\times n$ matrix has (at most) $n$ eigenvalues, $\lambda_1, \lambda_2, \ldots \lambda_n$, with eigenvectors $\mathbf{v}_1, \ldots \mathbf{v}_n$
- For most matrices, the eigenvectors form a **basis**, so for any vector $\mathbf{x}$,
\[
\mathbf{x} = \sum_{i=1}^{n}{c_i \mathbf{v}_i}
\]
for some coefficients/coordinates $c_1, \ldots c_n$
    + The $c_i$ says how much $\mathbf{x}$ projects on to $\mathbf{v}_i$
    + i.e., how similar is $\mathbf{x}$ to $\mathbf{v}_i$?

## Shrinkage

- Apply this last part to linear smoothing:
\[
\mathbf{x} = \sum_{i=1}^{n}{c_i \mathbf{v}_i}
\]
- We break the data up in to different "components" that come from the influence matrix $\mathbf{w}$
- After smoothing:
\[
\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}
\]
- Components of the data which match eigenvectors $\mathbf{v}_i$ with big eigenvalues $\lambda_i$ get enhanced/amplified by smoothing
- Components of the data which match eigenvectors with small eigenvalues get **shrunk**
    + If $\lambda_i = 0$, that component of the data is completely eliminated by smoothing
- Usually, the biggest eigenvalue of $\mathbf{w}$ will be 1
    + For a linear trend (=linear regression), every eigenvalue is either 0 or 1


## A little example {.smaller}

- Smoothing by averaging to either side
    + Need to handle the end-points...
```{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 {.smaller}

- Eigenvalues:
```{r}
eigen(w)$values
```

- Leading eigenvector:
```{r}
eigen(w)$vectors[,1]
```
- This is _constant_
    + If every entry in $\mathbf{x}$ is the same, averaging changes nothing!

##

```{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")
rounded.eigenvalues <- signif(eigen(w)$values[2:4], 2)
legend("bottomleft", legend=c(paste("lambda=", rounded.eigenvalues[1]), paste("lambda=", rounded.eigenvalues[2]),
                           paste("lambda=", rounded.eigenvalues[3])),
       lty=c("solid","dashed","dotted"))
```


## Variance of the fitted values

- Remember, fitted values are random:
\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}
- Unfortunately we don't know $\Var{\mathbf{\TrueNoise}}$
- **IF** we can assume $\Var{\mathbf{\TrueNoise}} = \sigma^2 \mathbf{I}$, then
\[
\Var{\EstRegFunc} = \sigma^2 \mathbf{w}\mathbf{w}^T
\]
as a special case

## Degrees of freedom {.smaller}

- How much do the fitted values co-vary with 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** $\Cov{\TrueNoise_i, \TrueNoise_j} = 0$ if $i\neq j$, and $=\sigma^2$ if $i=j$, _then_ as a special case
\[
\sum_{i=1}^{n}{\Cov{\EstRegFunc_i, X_i}} = \sigma^2 \tr{\mathbf{w}}
\]
- People use this to _define_ the **effective number of degrees of freedom**:
\[
\dof{\mathbf{w}} \equiv \tr{\mathbf{w}}
\]
    + Remember $\tr{\mathbf{w}} \equiv \sum_{i=1}^{n}{w_{ii}}$
    + Also remember $\tr{\mathbf{w}} =$ sum of the eigenvalues of $\mathbf{w}$

## Using linear smoothers

- Pros:
    + Don't need a theory about what the trend should be
    + Can _discover_ the underlying trend
    + Can be used to _check_ theories which make predictions about the trend
    + Mostly straightforward calculations
- Cons:
    + Doesn't use the data as efficiently as a theory-based model
        * ...when the model is right
    + Need to decide _how much_ smoothing/averaging to do
    + Need to decide _what kind_ of smoothing/averaging to do
        * Usually less important than the "how much" question
    + Extrapolation is hazardous

## Linear and polynomial trends not recommended {.smaller}

- We could just linearly regress $X(t)$ on $t$
- This'd give a **linear time trend**
     + Same ideas for space but more notation
- This is in fact a linear smoother but with weird weights:
\[
w(t, t_j) \propto (t-\overline{t})(t_j-\overline{t})
\]
- Unless you have an actual theory saying things should be linear in time, this is usually not a great idea
- Linearly regressing on $t$ plus $t^2$ plus $t^3$ etc. gives a **polynomial time trend**
- Usually an even worse idea
    + Polynomials behave badly away from the training data
- Extrapolating these trends is not something responsible professionals do

## Linear and polynomial trends not recommended {.smaller}

![](https://pbs.twimg.com/media/EXQlsS6WAAAXy3W?format=jpg&name=medium)

[(Official White House Council of Economic Advisers tweet, 2020-05-05)](https://twitter.com/WhiteHouseCEA/status/1257680257332674560)

## Linear and polynomial trends not recommended {.smaller}

- The CEA got a lot of push-back for fitting a cubic polynomial to the log of deaths
    + A polynomial _has_ to go to $\pm \infty$ as $t\rightarrow\infty$
    + This one was going to $-\infty$ so after exponentiating it says number of deaths $\rightarrow 0$
    + Sample push-back: ["haha cubic fitter go brrr"](https://www.facebook.com/Rmemes0/posts/cubic-fit-goes-brrr-brrrsource-httpstwittercomkjhealystatus1257733138370957315ra/3184714764914126/)
- To be quite fair, all the Institute for Health Metrics and Evaluation at the University of Washington was doing was fitting a quadratic polynomial!
    + Equivalent to saying that number of deaths on day $t$ is $\propto e^{-(t-t_0)^2/\tau}$
    + This has no more basis than the cubic but was somehow taken more seriously

## Summing up

- Every random process can be written as deterministic trend plus mean-zero random noise
- If we see the process multiple times, we can estimate the trend by averaging
- If we see the process only once, we can only estimate the trend with assumptions
- If we assume a particular theoretical model for the trend, we can estimate the parameters in that model and extrapolate
- If we assume the trend is smooth, we can use linear smoothers (= averaging) to estimate the trend
- We can understand a lot of what a linear smoother does by understanding the eigenvalues and eigenvectors of its influence matrix
- Extrapolating linear smoothers is rarely a good idea


## References
