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