5
$\begingroup$

Recently I was comparing the output of LOWESS regressions performed in R (and using Python's statsmodels module) and Stata. I realized that some of the values obtained by Stata seem to be off; specifically, it's the tails that seem to be estimated incorrectly.

I dove into the source code of the R's lowess() function (which seems to be based on Cleveland's original Fortran code found here) and the legacy Stata code for the ksm program found in the Stata 7 ado update file (found here). Note that post-v7 Statas implemented the lowess command using the _LOWESS C routine that is not exposed to the user. I verified that running Stata 7's ksm command using the optional lowess argument in Stata 14 generates the same ("incorrect") output as running lowess directly. That is, ksm Y X, lowess is the same as lowess Y X.

After inspecting the respective source codes I think the problem with Stata's implementation lies in fact that the subset size used in each local regression depends on the ordinal position of the X values. That is, for extreme and near-extreme values of X (in other words, for values close to the tails), Stata uses a smaller subset than for more central X's. Intuitively, the problem can be illustrated using a simple example with 100 data points where the bandwidth parameter is chosen to be 0.4 so that each subset is of the size 0.4*100=40. In R, the size of the subsets used for estimating $Y_1$, $Y_{10}$, $Y_{20}$, $Y_{30}$ and $Y_{40}$ would look like this:

R's subsets

On the other hand, Stata seems to reduce the size of the subsets for all X's from (roughly) 1 to 20 and from 80 to 100:

Stata's subsets

When I rewrote the ksm program such that the subsets size was held fixed for all X's, I got the same results as in R (or Python). Below are the pseudo-codes of both implementations:

# bwidth = bandwidth
# count = N
# i = index of local regression considered

# Cleveland's Fortran implementation
size = max(min(floor(bwidth*count), count), 2) #size of neighborhood
lower_bound = 1
upper_bound = size
for each i
    x = X[i]

    diff_left = x - X[lower_bound]
    diff_right = X[upper_bound + 1] - x
    if diff_left > diff_right and diff_right < count
        increment lower_bound and upper_bound by 1

    h = max(x - X[lower_bound], X[upper_bound] - x)
    h9 = 0.999*h and h1 = 0.001*h
    r = abs(X - x)
    W = 1 - (abs(X - x)/h)^3)^3 if r > h1 and r <= h9, W = 1 if r <= h1, W = 0 otherwise
    run WLS of Y on X using weights W if within lower_bound/upper_bound, get Y_hat[i]

# Stata implementation
k = int((bwidth*count - 0.5)/2) #half-bandwidth
for each i
    x = X[i]

    lower_bound = max(1, i - k)
    upper_bound = min(i + k, count)

    h = 1.0001*max(X[upper_bound] - x, x - X[lower_bound])
    W = 1 - (abs(X - x)/h)^3)^3 if h != 0, W = 1 otherwise
    run WLS of Y on X using weights W if within lower_bound/upper_bound, get Y_hat[i]

I understand that sometimes the implementations differ across software tools (or even packages made for the same software). That being said, I would expect to find at least some online references discussing such non-standard approaches. I tried to find as many other implementations of the LOWESS regression as I could and they all seem to be based on the Fortran code mentioned above and/or at least follow the same approach. I was not able to find any references discussing the specific approach used by Stata.

I would like to know whether such an implementation of LOWESS is correct, whether it's some "proprietary", quasi-correct method or simply a bug. In the former case, I would also appreciate if someone could please point me to a (preferably academic) reference. Thanks!

EDIT I: I wonder if this is because historically Stata's lowess command was based on the ksm program, which (when used with its default settings) was not meant to estimate the LOWESS regression. It appears that since Stata 8, the decision was made to abandon ksm and instead implement lowess where the default was chosen to be the LOWESS regression (whereas ksm's default options were implemented as lowess's non-defaults). The treatment of the tails, however, remained the same as in ksm.

EDIT II: Note that all of these comparisons relied on matching all other parameters of the LOWESS regressions. As Stata doesn't allow for multiple iterations and also doesn't implement the interpolation to cut down on the number of local regressions required (as do R and the original Fortran code), I made sure R's command used both iter=0 and delta=0 parameters. Specifically, I was comparing lowess(X, Y, f=0.4, delta=0, iter=0) (R) with lowess Y X, bwidth(0.4) (Stata) and statsmodels.api.nonparametric.lowess(Y, X, frac=0.4, it=0, delta=0) (Python). Both X and Y were well-behaved.

UPDATE: OK, it seems I wasn't the first to notice this. I just found an old user-written adjksm command, which "is identical to ksm except that the bandwidth of the smoother is constant along the x-axis." See here. I also confirmed that adjksm Y X, bwidth(0.4) lowess generates the same output as R's lowess() using iter=0 and delta=0, up to the 4th decimal point.

$\endgroup$
  • $\begingroup$ Interesting. I'm not actually convinced that Stata's implementation is incorrect (incorrect being not the exact definition of the LOWESS algorithm), even if R's version is more stable. I found that R's STL algorithm, I think also relying FORTRAN code written by Cleveland, does not follow the exact definition of the STL algorithm, but it's open to discussion whether that is a bad thing. $\endgroup$ – Cliff AB Feb 19 '17 at 21:53
  • 1
    $\begingroup$ Cleveland's own definitions went through a couple of stages of evolution (including a number of changes when he changed the name from LOWESS to LOESS). I don't think he ever did that bit at the ends with it, but it's probably a mistake to treat LOWESS as a specific unchanging thing rather than an evolving collection of related things, unless you specify something like "I mean LOWESS as defined in this particular paper" (and in some cases you need to further specify which options you mean to include) $\endgroup$ – Glen_b -Reinstate Monica Feb 19 '17 at 23:29
  • 1
    $\begingroup$ Thanks, point taken. As mentioned above, there's nothing wrong with different implementations of algorithms. But as someone who was taught to be critical of other people's work (be it research, code, etc.), I do mind when it cannot be easily scrutinized due to its closed-source nature yet it departs from the canonical form (and is not very transparent about it). But I guess that's already a topic for another discussion. $\endgroup$ – akavalar Feb 21 '17 at 23:29
2
$\begingroup$

Interesting investigation about details that I never paid attention to.

It's not a bug, because it is documented this way in the Stata documentation.

I never looked at the details and references for boundary problems in kernel regression, but it looks like Stata uses a different trade-off between bias and variance at the boundary than R or statsmodels. My guess is that R and statsmodels oversmooth compared to Stata if there are fluctuations or trends close to the boundaries. With constant bandwidth the effective window size would always shrink when going to the boundary, i.e. the Stata behavior results.

So maybe it's a partial misnomer in Stata if there is a strict definition of LOWESS, but Stata's behavior still looks common for windowed convolution filter or similar local polynomial kernel regression with fixed bandwidth.

In this case fixed bandwidth is in terms of neighbors and not in terms of metric distance, where Stata does not trade off neighbors from one side to the other side by extending the one-sided window.

For references I would look at Fan and Gijbels and related literature on boundary problems of kernel or local polynomial regression.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Interesting, thanks for chiming in. It was indeed the "boundary oversmoothing" that caught my attention, and I agree that calling this a bug is probably not entirely correct. I do find it quite frustrating that they've implemented this in such a non-standard way as it makes it harder to replicate one's results in other software. I still wonder if this was really intentional, or merely a byproduct/oversight of switching from "running-mean and running-line least-squares smoothing" to LOWESS - see Edit I and the subtle note (p9) of this. $\endgroup$ – akavalar Feb 19 '17 at 20:49
  • $\begingroup$ I don't know what the motivation of Stata developers is for this choice. As statsmodels developer I don't always like the choices of other packages, and then I don't feel compelled to follow them, at least not for the defaults. Sometimes there are also some computational shortcuts that differ across packages. $\endgroup$ – Josef Feb 19 '17 at 21:27
  • $\begingroup$ Related to design choices for LOWESS: I don't know what the tie handling is supposed to be. statsmodels' lowess has problems when there are a large number of observations with identical x, which do not show up if the underlying data is continuous, e.g. should neighbors include all observations that are tied at the boundary? $\endgroup$ – Josef Feb 19 '17 at 21:29
  • $\begingroup$ Another possibility about Stata's choice. Maybe it was initially an unintentional oversight or wrong analogy, but if it works better in some cases, it might be an intentional choice to keep it that way. $\endgroup$ – Josef Feb 19 '17 at 21:53
  • $\begingroup$ Thanks Josef! I don’t think Stata takes ties into account either. $\endgroup$ – akavalar Feb 21 '17 at 23:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.