33-smoothing.qmd

# Smoothing

```{r, echo=FALSE}
img_path <- "img"
```

* Other names given to this technique are *curve fitting* and *low pass filtering*. 

* It is designed to detect trends in the presence of noisy data in cases in which the shape of the trend is unknown. 

* The *smoothing* name comes from the fact that to accomplish this feat, we assume that the trend is *smooth*, as in a smooth surface. 
* In contrast, the noise, or deviation from the trend, is unpredictably wobbly:

```{r signal-plus-noise-example, message=FALSE, warning=FALSE, fig.height=6, echo=FALSE}
library(tidyverse)
set.seed(1)
n <- 100
x <- seq(-pi*4, pi*4, len = n)
tmp <- data.frame(x = x , f = sin(x) + x/8, e = rnorm(n, 0, 0.5)) 
p1 <- qplot(x, f, main = "smooth trend", ylim = range(tmp$f + tmp$e), data = tmp, geom = "line")
p2 <- qplot(x, e, main = "noise", ylim = range(tmp$f + tmp$e), data = tmp, geom = "line")
p3 <- qplot(x, f + e, main = "data = smooth trend + noise", ylim = range(tmp$f + tmp$e), data = tmp, geom = "line")
gridExtra::grid.arrange(p1, p2, p3)
```

* Part of what we explain in this section are the assumptions that permit us to extract the trend from the noise.



## Simplified MNIST: Is it a 2 or a 7? {#sec-two-or-seven}

* we define the challenge as building an algorithm that can determine if a digit is a 2 or 7 from the proportion of dark pixels in the upper left quadrant ($X_1$) and the lower right quadrant ($X_2$). 

* We also selected a random sample of 1,000 digits, 500 in the training set and 500 in the test set. We provide this dataset in the `dslabs` package:

```{r two-or-seven-scatter, warning=FALSE, message=FALSE, cache=FALSE}
library(tidyverse)
library(caret)
library(dslabs)
mnist_27$train |> ggplot(aes(x_1, x_2, color = y)) + geom_point()
```

* We can immediately see some patterns. For example, if $X_1$ (the upper left panel) is very large, then the digit is probably a 7. 

* Also, for smaller values of $X_1$, the 2s appear to be in the mid range values of $X_2$. 

To illustrate how to interpret $X_1$ and $X_2$, we include four example images.

*  On the left are the original images of the two digits with the largest and smallest values for $X_1$ and on the right we have the images corresponding to the largest and smallest values of $X_2$:

```{r two-or-seven-images-large-x1, echo=FALSE, out.width="100%", fig.height=3, fig.width=6.5}
if (!exists("mnist")) mnist <- read_mnist()
is <- mnist_27$index_train[c(which.min(mnist_27$train$x_1), which.max(mnist_27$train$x_1))]
titles <- c("smallest","largest")
tmp <- lapply(1:2, function(i){
    expand.grid(Row = 1:28, Column = 1:28) |>  
      mutate(label = titles[i],  
             value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
p1 <- tmp |> ggplot(aes(Row, Column, fill = value)) + 
  geom_raster(show.legend = FALSE) + 
  scale_y_reverse() +
  scale_fill_gradient(low = "white", high = "black") +
  facet_grid(.~label) + 
  geom_vline(xintercept = 14.5) +
  geom_hline(yintercept = 14.5) +
  ggtitle("Largest and smallest x_1")

    
is <- mnist_27$index_train[c(which.min(mnist_27$train$x_2), which.max(mnist_27$train$x_2))]
titles <- c("smallest","largest")
tmp <- lapply(1:2, function(i){
    expand.grid(Row = 1:28, Column = 1:28) |>  
      mutate(label = titles[i],  
             value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
p2 <- tmp |> ggplot(aes(Row, Column, fill = value)) + 
    geom_raster(show.legend = FALSE) + 
    scale_y_reverse() +
    scale_fill_gradient(low = "white", high = "black") +
    facet_grid(.~label) + 
    geom_vline(xintercept = 14.5) +
    geom_hline(yintercept = 14.5) +
  ggtitle("Largest and smallest x_2")
gridExtra::grid.arrange(p1, p2, ncol = 2)
```

* We can start getting a sense for why these predictors are useful, but also why the problem will be somewhat challenging.

* We haven't really learned any algorithms yet, so let's try building an algorithm using multivariable regression.

*  The model is simply:

$$
p(x_1, x_2) = \mbox{Pr}(Y=1 \mid X_1=x_1 , X_2 = x_2) = 
\beta_0 + \beta_1 x_1 + \beta_2 x_2
$$

* We fit it like this:

```{r}
fit <- mnist_27$train |> mutate(y = ifelse(y == 7, 1, 0)) |> lm(y ~ x_1 + x_2, data = _)
```

* We can now build a decision rule based on the estimate of $\hat{p}(x_1, x_2)$:

```{r}
p_hat <- predict(fit, newdata = mnist_27$test)
y_hat <- factor(ifelse(p_hat > 0.5, 7, 2))
confusionMatrix(y_hat, mnist_27$test$y)$overall[["Accuracy"]]
```

* We get an accuracy well above 50%. Not bad for our first try.

*  But can we do better?

* Because we constructed the `mnist_27` example and we had at our disposal 60,000 digits in just the MNIST dataset, we used this to build the _true_ conditional distribution $p(x_1, x_2)$. 

* Keep in mind that in practice we don't have access to the true conditional distribution.

*  We include it in this educational example because it permits the comparison of $\hat{p}(x_1, x_2)$ to the true $p(x_1, x_2)$. This comparison teaches us the limitations of different algorithms.

*  We have stored the true $p(x_1,x_2)$ in the `mnist_27` and can plot it as an image.

* We  draw a curve that separates pairs $(x_1,x_2)$ for which $p(x_1,x_2) > 0.5$ and pairs for which $p(x_1,x_2) < 0.5$:

```{r true-p-better-colors, echo=FALSE}
mnist_27$true_p |> 
  ggplot(aes(x_1, x_2, z = p)) +
  geom_raster(aes(fill = p)) +
  scale_fill_gradientn(colors = c("#F8766D", "white", "#00BFC4")) +
  stat_contour(breaks = c(0.5), color = "black")
```

* To start understanding the limitations of regression, first note that with regression $\hat{p}(x_1,x_2)$ has to be a plane, and as a result the boundary defined by the decision rule is given by: 
$\hat{p}(x_1,x_2) = 0.5$:

$$
\hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 = 0.5 \implies
\hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 = 0.5  \implies
x_2 = (0.5-\hat{\beta}_0)/\hat{\beta}_2  -\hat{\beta}_1/\hat{\beta}_2 x_1
$$

* This implies that for the boundary $x_2$ is a linear function of $x_1$. 

* This implies that our regression approach has no chance of capturing the non-linear nature of the true $p(x_1,x_2)$. 

* Below is a visual representation of $\hat{p}(x_1, x_2)$ which clearly shows how it fails to capture the shape of $p(x_1, x_2)$.

```{r regression-p-hat, echo=FALSE, out.width="100%", fig.height=3, fig.width=7}
p_hat <- predict(fit, newdata = mnist_27$true_p)
p_hat <- scales::squish(p_hat, c(0, 1))
p1 <- mnist_27$true_p |> mutate(p_hat = p_hat) |>
  ggplot(aes(x_1, x_2,  z = p_hat)) +
  geom_raster(aes(fill = p_hat)) +
  scale_fill_gradientn(colors = c("#F8766D","white","#00BFC4")) +
  stat_contour(breaks = c(0.5), color = "black") 

p2 <- mnist_27$true_p |> 
  mutate(p_hat = p_hat) |>
  ggplot() +
  stat_contour(aes(x_1, x_2, z = p_hat), breaks = c(0.5), color = "black") + 
  geom_point(mapping = aes(x_1, x_2, color = y), data = mnist_27$test) 
gridExtra::grid.arrange(p1, p2, ncol = 2)
```

* We need something more flexible: a method that permits estimates with shapes other than a plane.

*  Smoothing techniques permit this flexibility.

*  We will start by describing nearest neighbor and kernel approaches.

*  To understand why we cover this topic, remember that the concepts behind smoothing techniques are extremely useful in machine learning because conditional expectations/probabilities can be thought of as *trends* of unknown shapes that we need to estimate in the presence of uncertainty.


## Signal plus noise model

* To explain these concepts, we will focus first on a problem with just one predictor.

*  Specifically, we try to estimate the time trend in the 2008 US popular vote poll margin (difference between Obama and McCain). Later we will learn about methods such as k-nearest neighbors that can be used to smooth with higher dimensions. 


```{r polls-2008-data, warning=FALSE, message=FALSE, cache=FALSE}
polls_2008 |> ggplot(aes(day, margin)) + geom_point()
```


* For the purposes of the popular vote example, do not think of it as a forecasting problem.

*  Instead, we are simply interested in learning the shape of the trend *after* the election is over.

* We assume that for any given day $x$, there is a true preference among the electorate $f(x)$, but due to the uncertainty introduced by the polling, each data point comes with an error $\varepsilon$. A mathematical model for the observed poll margin $Y_i$ is:

$$
Y_i = f(x_i) + \varepsilon_i
$$

* To think of this as a machine learning problem, consider that we want to predict $Y$ given a day $x$. If we knew the conditional expectation $f(x) = \mbox{E}(Y \mid X=x)$, we would use it.

*  But since we don't know this conditional expectation, we have to estimate it.

*  Let's use regression, since it is the only method we have learned up to now.

```{r linear-regression-not-flexible, echo=FALSE}
resid <- ifelse(lm(margin~day, data = polls_2008)$resid > 0, "+", "-")
polls_2008 |> 
  mutate(resid = resid) |> 
  ggplot(aes(day, margin)) + 
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  geom_point(aes(color = resid), size = 3)
```

* The fitted regression line does not appear to describe the trend very well.

*  For example, on September 4 (day -62), the Republican Convention was held and the data suggest that it gave John McCain a boost in the polls.

*  However, the regression line does not capture this potential trend.

*  To see the *lack of fit* more clearly, we note that points above the fitted line (blue) and those below (red) are not evenly distributed across days.

*  We therefore need an alternative, more flexible approach.

## Bin smoothing

* The general idea of smoothing is to group data points into strata in which the value of $f(x)$ can be assumed to be constant.

*  We can make this assumption when we think $f(x)$ changes slowly and, as a result, $f(x)$ is almost constant in small windows of $x$. An example of this idea for the `poll_2008` data is to assume that public opinion remained approximately the same within a week's time.

*  With this assumption in place, we have several data points with the same expected value.

* If we fix a day to be in the center of our week, call it $x_0$, then for any other day $x$ such that $|x - x_0| \leq 3.5$, we assume $f(x)$ is a constant $f(x) = \mu$. This assumption implies that: 

$$
E[Y_i | X_i = x_i ] \approx \mu \mbox{   if   }  |x_i - x_0| \leq 3.5
$$

* In smoothing, we call the size of the interval satisfying $|x_i - x_0| \leq 3.5$ the *window size*, *bandwidth* or *span*. Later we will see that we try to optimize this parameter.

* This assumption implies that a good estimate for $f(x_0)$ is the average of the $Y_i$ values in the window.

*  If we define $A_0$ as the set of indexes $i$ such that $|x_i - x_0| \leq 3.5$ and $N_0$ as the number of indexes in $A_0$, then our estimate is:

$$
\hat{f}(x_0) = \frac{1}{N_0} \sum_{i \in A_0}  Y_i
$$

* We make this calculation with each value of $x$ as the center.

*  In the poll example, for each day, we would compute the average of the values within a week with that day in the center.

*  Here are two examples: $x_0 = -125$ and $x_0 = -55$. The blue segment represents the resulting average.

```{r binsmoother-expained, echo=FALSE}
span <- 3.5
tmp <- polls_2008 |>
  crossing(center = polls_2008$day) |>
  mutate(dist = abs(day - center)) |>
  filter(dist <= span) 

tmp |> filter(center %in% c(-125, -55)) |>
  ggplot(aes(day, margin)) +   
  geom_point(data = polls_2008, size = 3, alpha = 0.5, color = "grey") +
  geom_point(size = 2) +    
  geom_smooth(aes(group = center), 
              method = "lm", formula = y~1, se = FALSE) +
  facet_wrap(~center)
```

* By computing this mean for every point, we form an estimate of the underlying curve $f(x)$. Below we show the procedure happening as we move from the -155 up to 0. 

* At each value of $x_0$, we keep the estimate $\hat{f}(x_0)$ and move on to the next point:

```{r binsmoother-animation, echo=FALSE, warning=FALSE}
knitr::include_graphics("https://rafalab.dfci.harvard.edu/dsbook-part-2/ml/img/binsmoother-animation.gif")
```

* The final code and resulting estimate look like this:

```{r binsmoother-final}
span <- 7 
fit <- with(polls_2008, ksmooth(day, margin, kernel = "box", bandwidth = span))

polls_2008 |> mutate(smooth = fit$y) |>
  ggplot(aes(day, margin)) +
    geom_point(size = 3, alpha = .5, color = "grey") + 
  geom_line(aes(day, smooth), color = "red")
```

## Kernels

* The final result from the bin smoother is quite wiggly.

*  One reason for this is that each time the window moves, two points change.

*  We can attenuate this somewhat by taking weighted averages that give the center point more weight than far away points, with the two points at the edges receiving very little weight.

* You can think of the bin smoother approach as a weighted average:

$$
\hat{f}(x_0) = \sum_{i=1}^N w_0(x_i) Y_i
$$

in which each point receives a weight of either $0$ or $1/N_0$, with $N_0$ the number of points in the week.

*  In the code above, we used the argument `kernel="box"` in our call to the function `ksmooth`. This is because the weight function looks like a box.

*  The `ksmooth` function provides a "smoother" option which uses the normal density to assign weights.

```{r gaussian-kernel, echo=FALSE, out.width="80%", fig.height=3, fig.width=6}
x_0 <- -125

p1 <- data.frame(x = polls_2008$day) |> mutate(w_0 = 1*I(abs(x - x_0) <= span/2)) |>
  mutate(w_0 = w_0/sum(w_0)) |>
  ggplot(aes(x, w_0)) +
  geom_step() +
  ggtitle("Box")

tmp <- with(data.frame(day = seq(min(polls_2008$day), max(polls_2008$day), .25)), 
                       ksmooth(day, 1*I(day == x_0), kernel = "normal", x.points = day, bandwidth = span))
p2 <- data.frame(x = tmp$x, w_0 = tmp$y) |>
  mutate(w_0 = w_0/sum(w_0)) |>
  ggplot(aes(x, w_0)) +
  geom_line() +
  ggtitle("Normal")
gridExtra::grid.arrange(p1, p2, ncol = 2)
```

* The final code and resulting plot for the normal kernel look like this:

```{r final-ksmooth-normal-kernel}
span <- 7
fit <- with(polls_2008, ksmooth(day, margin, kernel = "normal", bandwidth = span))

polls_2008 |> mutate(smooth = fit$y) |>
  ggplot(aes(day, margin)) +
  geom_point(size = 3, alpha = .5, color = "grey") + 
  geom_line(aes(day, smooth), color = "red")
```

* Notice that this version looks smoother.

* There are several functions in R that implement bin smoothers.

*  One example is `ksmooth`, shown above.

*  In practice, however, we typically prefer methods that use slightly more complex models than fitting a constant.

*  The final result above, for example, is still somewhat wiggly in parts we don't expect it to be (between -125 and -75, for example). Methods such as `loess`, which we explain next, improve on this.

## Local weighted regression (loess)

* A limitation of the bin smoother approach just described is that we need small windows for the approximately constant assumptions to hold.

*  As a result, we end up with a small number of data points to average and obtain imprecise estimates $\hat{f}(x)$. Here we describe how *local weighted regression* (loess) permits us to consider larger window sizes.

*  To do this, we will use a mathematical result, referred to as Taylor's theorem, which tells us that if you look closely enough at any smooth function $f(x)$, it will look like a line.

*  To see why this makes sense, consider the curved edges gardeners make using straight-edged spades:

![](https://rafalab.dfci.harvard.edu/dsbook-part-2/ml/img/garden.png)


* Instead of assuming the function is approximately constant in a window, we assume the function is locally linear.

*  We can consider larger window sizes with the linear assumption than with a constant.

*  Instead of the one-week window, we consider a larger one in which the trend is approximately linear.

*  We start with a three-week window and later consider and evaluate other options:

$$
E[Y_i | X_i = x_i ] = \beta_0 + \beta_1 (x_i-x_0) \mbox{   if   }  |x_i - x_0| \leq 21
$$

* For every point $x_0$, loess defines a window and fits a line within that window.

*  Here is an example showing the fits for $x_0=-125$ and $x_0 = -55$:

```{r loess, echo=FALSE}
span <- 21/diff(range(polls_2008$day))

tmp <- polls_2008 |>
  crossing(center = polls_2008$day) |>
  mutate(dist = abs(day - center)) |>
  filter(rank(dist) / n() <= span) |>
  mutate(weight = (1 - (dist / max(dist)) ^ 3) ^ 3)

tmp |> 
  filter(center %in% c(-125, -55)) |>
  ggplot(aes(day, margin)) +   
  scale_size(range = c(0, 3)) +
  geom_smooth(aes(group = center, weight = weight), 
              method = "lm", se = FALSE) +
  geom_point(data = polls_2008, size = 3, alpha = .5, color = "grey") +
  geom_point(aes(size = weight)) +
  facet_wrap(~center)
```

* The fitted value at $x_0$ becomes our estimate $\hat{f}(x_0)$. Below we show the procedure happening as we move from the -155 up to 0.

```{r loess-animation, echo=FALSE, warning=FALSE}
knitr::include_graphics("https://rafalab.dfci.harvard.edu/dsbook-part-2/ml/img/loess-animation.gif")
```

* The final result is a smoother fit than the bin smoother since we use larger sample sizes to estimate our local parameters:

```{r final-loess}
total_days <- diff(range(polls_2008$day))
span <- 21/total_days
fit <- loess(margin ~ day, degree = 1, span = span, data = polls_2008)
polls_2008 |> mutate(smooth = fit$fitted) |>
  ggplot(aes(day, margin)) +
  geom_point(size = 3, alpha = .5, color = "grey") +
  geom_line(aes(day, smooth), color = "red")
```

* Different spans give us different estimates.

*  We can see how different window sizes lead to different estimates:

```{r loess-multi-span-animation, echo=FALSE, warning=FALSE}
spans <- c(.66, 0.25, 0.15, 0.10)

fits <- tibble(span = spans) |> 
  group_by(span) |> 
  do(broom::augment(loess(margin ~ day, degree = 1, span = .$span, data = polls_2008)))

tmp <- fits |>
  crossing(center = polls_2008$day) |>
  mutate(dist = abs(day - center)) |>
  filter(rank(dist) / n() <= span) |>
  mutate(weight = (1 - (dist / max(dist)) ^ 3) ^ 3)

knitr::include_graphics("https://rafalab.dfci.harvard.edu/dsbook-part-2/ml/img/loess-multi-span-animation.gif")
```

* Here are the final estimates:

```{r loess-final, echo=FALSE}
tmp |> ggplot(aes(day, margin)) +
  geom_point(size = 2, alpha = .5, color = "grey") +
  geom_line(aes(day, .fitted), data = fits, color = "red") +
  facet_wrap(~span)
```

* There are three other differences between `loess` and the typical bin smoother.

1\. Rather than keeping the bin size the same, `loess` keeps the number of points used in the local fit the same.This number is controlled via the `span` argument, which expects a proportion.
For example, if `N` is the number of data points and `span=0.5`, then for a given $x$, `loess` will use the `0.5 * N` closest points to $x$ for the fit.

2\. When fitting a line locally, `loess` uses a *weighted* approach.Basically, instead of using least squares, we minimize a weighted version:

$$
\sum_{i=1}^N w_0(x_i) \left[Y_i - \left\{\beta_0 + \beta_1 (x_i-x_0)\right\}\right]^2
$$

However, instead of the Gaussian kernel, loess uses a function called the Tukey tri-weight:

$$
W(u)= \left( 1  - |u|^3\right)^3 \mbox{ if } |u| \leq 1 \mbox{ and } W(u) = 0 \mbox{ if } |u| > 1
$$

To define the weights, we denote $2h$ as the window size and define:

$$
w_0(x_i) = W\left(\frac{x_i - x_0}{h}\right)
$$

This kernel differs from the Gaussian kernel in that more points get values closer to the max:

```{r triweight-kernel, echo=FALSE, out.width="80%", fig.height=3, fig.width=6}
x_0 <- -125

tmp <- with(data.frame(day = seq(min(polls_2008$day), max(polls_2008$day), .25)), 
                       ksmooth(day, 1*I(day == x_0), kernel = "normal", x.points = day, bandwidth = 7))
p1 <- data.frame(x = tmp$x, w_0 = tmp$y) |>
  mutate(w_0 = w_0/sum(w_0)) |>
  ggplot(aes(x, w_0)) +
  geom_line() +
  ggtitle("Normal")

p2 <- data.frame(x = seq(min(polls_2008$day), max(polls_2008$day), length.out = 100)) |>
  mutate(w_0 = (1 - (abs(x - x_0)/21)^3)^3*I(abs(x - x_0) <= 21)) |>
  ggplot(aes(x, w_0)) +
  geom_line()  +
  ggtitle("Tri-weight")
gridExtra::grid.arrange(p1, p2, ncol = 2)
```

3\. `loess` has the option of fitting the local model *robustly*. An iterative algorithm is implemented in which, after fitting a model in one iteration, outliers are detected and down-weighted for the next iteration. To use this option, we use the argument `family="symmetric"`.

### Fitting parabolas

* Taylor's theorem also tells us that if you look at any mathematical function closely enough, it looks like a parabola.

*  The theorem also states that you don't have to look as closely when approximating with parabolas as you do when approximating with lines.

*  This means we can make our windows even larger and fit parabolas instead of lines.

$$
E[Y_i | X_i = x_i ] = \beta_0 + \beta_1 (x_i-x_0) + \beta_2 (x_i-x_0)^2 \mbox{   if   }  |x_i - x_0| \leq h
$$

* You may have noticed that when we showed the code for using loess, we set `degree = 1`. This tells loess to fit polynomials of degree 1, a fancy name for lines.

*  If you read the help page for loess, you will see that the argument `degree` defaults to 2. By default, loess fits parabolas not lines.

*  Here is a comparison of the fitting lines (red dashed) and fitting parabolas (orange solid):

```{r polls-2008-parabola-line-loess}
total_days <- diff(range(polls_2008$day))
span <- 28/total_days
fit_1 <- loess(margin ~ day, degree = 1, span = span, data = polls_2008)
fit_2 <- loess(margin ~ day, span = span, data = polls_2008)

polls_2008 |> mutate(smooth_1 = fit_1$fitted, smooth_2 = fit_2$fitted) |>
  ggplot(aes(day, margin)) +
  geom_point(size = 3, alpha = .5, color = "grey") +
  geom_line(aes(day, smooth_1), color = "red", lty = 2) +
  geom_line(aes(day, smooth_2), color = "orange", lty = 1) 
```

* The `degree = 2` gives us more wiggly results.

*  In general, we actually prefer `degree = 1` as it is less prone to this kind of noise.

### Beware of default smoothing parameters

`ggplot` uses loess in its `geom_smooth` function:

```{r ggplot-loess-default, warning=FALSE, message=FALSE}
polls_2008 |> ggplot(aes(day, margin)) +
  geom_point() + 
  geom_smooth(method = loess)
```

* But be careful with default parameters as they are rarely optimal.

*  However, you can conveniently change them:

```{r ggplot-loess-degree-1, warning=FALSE, message=FALSE}
polls_2008 |> ggplot(aes(day, margin)) +
  geom_point() + 
  geom_smooth(method = loess, method.args = list(span = 0.15, degree = 1))
```

## Connecting smoothing to machine learning {#sec-smoothing-ml-connection}

* To see how smoothing relates to machine learning with a concrete example, consider again our @sec-two-or-seven example.

*  If we define the outcome $Y = 1$ for digits that are seven and $Y=0$ for digits that are 2, then we are interested in estimating the conditional probability:

$$
p(x_1, x_2) = \mbox{Pr}(Y=1 \mid X_1=x_1 , X_2 = x_2).
$$ 

* with $X_1$ and $X_2$ the two predictors defined above. 

* In this example, the 0s and 1s we observe are "noisy" because for some regions the probabilities $p(x_1, x_2)$ are not that close to 0 or 1. So we need to estimate $p(x_1, x_2)$. 

* Smoothing is an alternative to accomplishing this.

* We saw that linear regression was not flexible enough to capture the non-linear nature of $p(x_1, x_2)$, thus smoothing approaches provide an improvement.

*  In a future lecture we describe a popular machine learning algorithm, k-nearest neighbors, which is based on bin smoothing.


## Case Study: Estimating indirect effects of Natural Dissasters


```{r, eval=FALSE}
head(pr_death_counts)
```

Explore the data and notice we should remove missing or incomplete values

```{r}
pr_death_counts |>
  ggplot(aes(date, deaths)) +
  geom_point()

dat <- pr_death_counts |> filter(date < make_date(2018, 4, 1))
```


We start by fitting a seasonal model. We will fit the model to 2015-2016 and then obtain residulas for 2017-2018

```{r}
dat <- mutate(dat, tt = as.numeric(date - make_date(2017, 9, 20))) |>
  filter(!is.na(deaths)) |>
  mutate(x1 = sin(2*pi*tt/365), x2 = cos(2*pi*tt/365),
         x3 = sin(4*pi*tt/365), x4 = cos(4*pi*tt/365)) 

seasonal_fit <- lm(deaths ~ x1 + x2 + x3 + x4, data = filter(dat, year(date) < 2017))
```

Check the fit

```{r}
with(dat, plot(tt, deaths))
s <- predict(seasonal_fit, newdata = dat)
lines(dat$tt, s, col = 2, lwd = 2)
```


Let's now fit smooth curve to the residuals

```{r}
dat <- mutate(dat, resid = deaths - s)
fit <- loess(resid ~ tt, span = 1/20, degree = 1, data = dat)
fit <- predict(fit, se = TRUE)
dat <- mutate(dat, f = fit$fit) |>
  mutate(lower = f - qnorm(0.995)*fit$se.fit, upper = f + qnorm(.995)*fit$se.fit)
dat |> ggplot(aes(x = date)) +
  geom_point(aes(y = resid)) +
  geom_line(aes(y = f), color = "red")
```

Is there some evidence of indirect effect after the first week after landfall?

```{r}
dat |> 
  filter(year(date) >= 2017) |>
  ggplot(aes(date, f, ymin = lower, ymax = upper)) +
  geom_ribbon(alpha = 0.5) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0) +
  theme_bw()
```

Note that this is only an exploratory analysis. A formal analysis would take into account changing demographics and perhaps fit a model with a discontinuity on the day the hurricane made landfall.