LOESS that allows discontinuities

• Is there a modelling technique like LOESS that allows for zero, one, or more discontinuities, where the timing of the discontinuities are not known apriori?
• If a technique exists, is there an existing implementation in R?
• discontinuities at known x-values, or at unknown x-values? (known x is easy enough) – Glen_b Sep 7 '10 at 7:04
• @glen I updated the question: I'm interested in situations where timing of discontinuities is not known apriori. – Jeromy Anglim Sep 7 '10 at 7:08
• This may be a moot/foolish question, but you say "timing": is this for use with time series? I believe most of the answers below assume this ("changepoint, etc"), though LOESS can be applied in non-time-series situations, with discontinuities. I think. – Wayne Dec 17 '10 at 16:52

It sounds like you want to perform multiple changepoint detection followed by independent smoothing within each segment. (Detection can be online or not, but your application is not likely to be online.) There's a lot of literature on this; Internet searches are fruitful.

• DA Stephens wrote a useful introduction to Bayesian changepoint detection in 1994 (App. Stat. 43 #1 pp 159-178: JSTOR).
• More recently Paul Fearnhead has been doing nice work (e.g., Exact and efficient Bayesian inference for multiple changepoint problems, Stat Comput (2006) 16: 203-213: Free PDF).
• A recursive algorithm exists, based on a beautiful analysis by D Barry & JA Hartigan
• Product Partition Models for Change Point Models, Ann. Stat. 20:260-279: JSTOR;
• A Bayesian Analysis for Change Point Problems, JASA 88:309-319: JSTOR.
• One implementation of the Barry & Hartigan algorithm is documented in O. Seidou & TBMJ Ourda, Recursion-based Multiple Changepoint Detection in Multivariate Linear Regression and Application to River Streamflows, Water Res. Res., 2006: Free PDF.

I haven't looked hard for any R implementations (I had coded one in Mathematica a while ago) but would appreciate a reference if you do find one.

• I found the bcp R package jstatsoft.org/v23/i03/paper which implements the Barry & Hartigan algorithm – Jeromy Anglim Sep 8 '10 at 2:39
• @Jeromy: Thank you for the R package and for inserting the links to the references. – whuber Sep 8 '10 at 14:32

do it with koencker's broken line regression, see page 18 of this vignette

http://cran.r-project.org/web/packages/quantreg/vignettes/rq.pdf

In response to Whuber last comment:

This estimator is defined like this.

$x\in\mathbb{R}$, $x_{(i)}\geq x_{(i-1)}\;\forall i$,

$e_i:=y_{i}-\beta_{i}x_{(i)}-\beta_0$,

$z^+=\max(z,0)$, $z^-=\max(-z,0)$,

$\tau \in (0,1)$, $\lambda\geq 0$

$\underset{\beta\in\mathbb{R}^n|\tau, \lambda}{\min.} \sum_{i=1}^{n} \tau e_i^++\sum_{i=1}^{n}(1-\tau)e_i^-+\lambda\sum_{i=2}^{n}|\beta_{i}-\beta_{i-1}|$

$\tau$ gives the desired quantile (i.e. in the example, $\tau=0.9$). $\lambda$ directs the number of breakpoint: for $\lambda$ large this estimator shrinks to no break point (corresponding to the classicla linear quantile regression estimator).

Quantile Smoothing Splines Roger Koenker, Pin Ng, Stephen Portnoy Biometrika, Vol. 81, No. 4 (Dec., 1994), pp. 673-680

PS: there is a open acess working paper with the same name by the same others but it's not the same thing.

• That's a neat idea: thanks for the reference. However, the residuals to that particular fit look pretty bad, which makes me wonder how well it identifies potential changepoints. – whuber Sep 7 '10 at 20:40
• whuber: i do not know how much you are familiar with the theory of quantile regression. These lines have a major advantage over splines: they do not assume any error distribution (i.e. they do not assume the residuals to be Gaussian). – user603 Sep 7 '10 at 20:58
• @kwak This looks interesting. Not assuming a normal error distribution would be useful for one of my applications. – Jeromy Anglim Sep 8 '10 at 2:36
• Indeed, what you get out of this estimation are actual conditionnal quantiles: in a nutshell, these are to splines/LOESS-regressions what boxplots are to the couple (mean, s.d.): a much richer view of your data. They also retain there validity in non gaussian context (such as assymetric errors,...). – user603 Sep 8 '10 at 2:51
• @kwak: The residuals are heavily correlated with the x-coordinate. For example, there are long runs of negative or small positive residuals. Whether they have a Gaussian distribution or not, then, is immaterial (as well as irrelevant in any exploratory analysis): this correlation shows that the fit is poor. – whuber Sep 8 '10 at 14:30

Here are some methods and associated R packages to solve this problem

Wavelet thresolding estimation in regression allows for discontonuities. You may use the package wavethresh in R.

A lot of tree based methods (not far from the idea of wavelet) are usefull when you have disconitnuities. Hence package treethresh, package tree !

In the familly of "local maximum likelihood" methods... among others: Work of Pozhel and Spokoiny: Adaptive weights Smoothing (package aws) Work by Catherine Loader: package locfit

I guess any kernel smoother with locally varying bandwidth makes the point but I don't know R package for that.

note: I don't really get what is the difference between LOESS and regression... is it the idea that in LOESS alrgorithms should be "on line" ?

• Re LOESS: Perhaps my terminology is not quite right. By LOESS I'm referring to models that predict Y from X using some form of localised curve fitting. e.g., as seen in most of these graphs: google.com/… – Jeromy Anglim Sep 7 '10 at 7:12

It should be possible to code a solution in R using the non-linear regression function nls, b splines (the bs function in the spline package, for example) and the ifelse function.