5.3 Local Regression

Kernel Smoothing¶

We’ve encountered data smoothing in the form of a moving average. Figure 1 and Figure 2 present a synthetic white noise process and the same process after smoothing with a moving average of length 10.

We can think of the smoothing process depicted in Figure 2 as consisting of sliding a kernel (also known as a filter) across the data. In this case, the kernel is defined as

which we can envision as being slid across the data. At each point, we take the inner (or dot) product of the data with the kernel. For the kernel used in Eq. (1) this gives the unweighted average of 10 points. In later chapters we will learn that this can be thought of as a convolution in the frequency domain.

The box (or boxcar) kernel used in Figure 2 is only one of many possible choices. We will learn in future chapters that a box filter has substantial shortcomings in the frequency domain that are mitigated (though not eliminated) by use of tapered kernels, such as kernels based on Gaussians, the cosine function, or the tricube function introduced in the next section.

Tricube Function¶

For this chapter, the most important kernel method is generated from the tricube function $W$:

where $u$ is the distance between the the focal point $x_0$ and the point $x_i$ calculated as $u=\frac{x_0-x_i}{\max_j{(|x_0-x_j|)}}$^[1]. The maximum $j$ is determined by an adjustable parameter giving the length of the kernel.

For example, for $\max_j{(|x_0-x_j|)}=5$, the tricube function $W$ generates a sharply peaked kernel with $2j-1=9$ non-zero values out of $2j+1=11$ total values

with all points outside the window weighted as 0 (note that for the first and last elements $u=\frac{\pm5}{5}$, resulting in zero for $W$). In actual implementation we would simply use a kernel of length $2j-1$ (9, in our example) instead of using a kernel of length $2j+1$ (11, in our example) with zeros in the first and last positions. We explicitly included the zeros in Eq. (3) to emphasize that the kernel falls to zero at $|u|=1$.

Weighting Points Near the Edges¶

In order to avoid shortening the time series (as occurs with moving averages), any point less than $\max_j{(|x_0-x_j|)}$ from either end has an asymmetric selection of points in order to maintain the same kernel size at the expense of a larger value for $\max_j{(|x_0-x_j|)}$. For example, continuing with $\max_j{(|x_0-x_j|)}=5$ from before, the asymmetric window used for the first point in the series will now have $\max_j{(|x_0-x_j|)}=9$ with values given as:

LOWESS Algorithm¶

Local regression is perhaps best thought of as a hybrid between kernel smoothing and linear regression (or more broadly polynomial regression). The most common family of methods is LOWESS (locally weighted scatterplot smoothing), a type of LOESS (locally estimated scatterplot smoothing)^[2]. The LOWESS algorithm is generally credited as having been formally introduced in Cleveland (1979) and Cleveland (1981) (though much of the concept had been proposed independently by Savitzky & Golay (1964) about fifteen years beforehand).

LOWESS Demo¶

Before diving into the underlying algorithm, let’s get a feeling for how LOWESS behaves. The following interactive figure plots the function

for some noise $w$ in gray dots. For reference, the original (presumably unknown) clean signal without noise that we are attempting to reconstruct is denoted by a gray dashed line. The toggle button LOWESS CURVE turns on or off the LOWESS curve fitted through the points. The button LOCAL FITS toggles on and off a selection of the individual local regression lines (more on those soon). Finally, the slider SPAN (f) tells the algorithm what fraction of the overall data to consider for each local fit. Note that at 1.00, the line is almost a flat linear regression line with large local regression lines. As you shrink the fraction down, the LOWESS fit becomes more jagged with smaller local regression lines (denoting that the local line is calculated with a smaller portion of the data).

If the above fails to render correctly in your browser you can also open the demo as a new browser window using the Open Demo in a New Tab ↗ button at the top of the frame. Note that you may need to enable popups for this to work.

Local Regression¶

The LOWESS algorithm begins by dividing the data into overlapping windows. The size of each window is governed by the parameter $f$, the fraction of the total observations $n$. In statsmodels.nonparametric.smoothers_lowess.lowess, the fraction is set with the argument frac (default frac=2.0/3.0). In the demo above, $f$ can be smoothly varied from 0.05 to 1. Once $f$ is chosen, we create $n$ windows of length $\lceil f\times n \rceil$ centered at each data point where $\lceil f\times n \rceil$ is the ceiling function indicating the smallest integer greater than or equal to $f\times n$. For each window about a given central point (denoted as $x_{center}$), we calculate a weighted linear regression defined as

where $W(x_{center}, x_i)$ is a weight function, most commonly the tricube function, that weights the contribution of each $x_i$ based on its distance from $x_{center}$ and where we have expressed the coefficients as $\beta_0(x_{center})$ and $\beta_1(x_{center})$ to emphasize that distinct $\beta$’s are calculated for each data point. For points closer than $\frac{f\times n}{2}$ to either end, the window no longer places the point in the center of the window, and uses the modification to $W$ discussed above.

LOWESS can be extended from linear regression to polynomial regression using expressions such as

though we will not need to go beyond Eq. (6) for our purposes.

Having calculated $n$ distinct linear regressions, we use each linear regression about a given $x_{center}$ to calculate the $\hat{y}_{center}$ at that point, leaving us with $n$ $\hat{y}$ values (one per data point). The $\hat{y}$ points are then connected into a single curve—though note in the demo above that the fitted curve becomes notably more jagged as $f$ is decreased.

Robust Smoothing¶

In the demo above, we calculated the LOWESS curve solely by implementing Eq. (6). This method works well for data that has relatively few outliers. Cleveland (1979) describes an optional second stage of smoothing using the bisquare function $B(x)$ defined as

where as in Eq. (2) $B(x)$ falls to zero at $x=\pm1$. After the initial LOWESS curve calculation using Eq. (6), the algorithm then goes through a predetermined number of iterations of the following procedure:

1. Calculate the residual between the predicted and observed points

   $$e_i = y_i -\hat{y}_i.$$

   (9)

2. Determine new weightings $\delta$ for each point by first calculating the median $s$ of the $|e_i|$'s and then applying the formula

   $$\delta_k \stackrel{\triangle}{=} B\Big(\frac{e_k}{6s}\Big).$$

   (10)

3. Recalculate the local regressions using Eq. (6) using $B\times W$ instead of $W$.

Cleveland refers to this algorithm as “robust locally weighted regression.” In statsmodels, the number of iterations of the above procedure is controlled by the argument it (default it=3). Setting it=0 will result in the algorithm only using Eq. (6) in its original form.

In general, robust smoothing is primarily required for datasets with large outliers, such as those produced by Cauchy distributed noise. Figure 3 demonstrates LOWESS smoothing with different numbers of iterations on a synthetic dataset with Cauchy distributed noise. Note in particular that for 0 iterations (non-robust LOWESS) the fitted LOWESS line “chases” the outliers out of the plot. Once we reach 3 iterations, the LOWESS line almost perfectly captures the original signal.