15
  • 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?
amoeba
  • 93,463
  • 28
  • 275
  • 317
Jeromy Anglim
  • 42,044
  • 23
  • 146
  • 250
  • 1
    discontinuities at known x-values, or at unknown x-values? (known x is easy enough) – Glen_b Sep 07 '10 at 07:04
  • @glen I updated the question: I'm interested in situations where timing of discontinuities is not known apriori. – Jeromy Anglim Sep 07 '10 at 07: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

4 Answers4

15

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.

whuber
  • 281,159
  • 54
  • 637
  • 1,101
7

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.

user603
  • 21,225
  • 3
  • 71
  • 135
  • 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 07 '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 07 '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 08 '10 at 02: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 08 '10 at 02: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 08 '10 at 14:30
  • Whuber: the line you see on the plot is an estimate of the 90% conditionnal quantile (i.e. the function $f(x):p(y_i – user603 Sep 08 '10 at 17:39
  • @Kwak: Thank you for the clarification. But if the estimate is true only asymptotically and globally, why should we expect it to perform at all well in identifying changepoints? Is there some additional property the broken line regression has, of which I am unaware, that provides some comfort in this respect? – whuber Sep 10 '10 at 02:36
  • @Whuber:> i appended my answer with some mathematical details and a much needed scientific citation. Best. – user603 Sep 10 '10 at 11:24
  • @Kwak: Thank you! I see how the case \tau = 0.5 might be effective. (BTW, you probably meant to write \tau instead of \alpha in the formula.) It would be interesting to compare its performance to methods specifically designed to identify changepoints. – whuber Sep 10 '10 at 14:02
  • @Whuber:> with the caveat that this is not an 'online' algorithm. If you add observations, the locations of the change points will vary (although this behavior is fairly easy to 'fix', for some reasons, the author have never explored that direction). – user603 Sep 10 '10 at 17:59
6

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

robin girard
  • 6,335
  • 6
  • 46
  • 60
  • 1
    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: http://www.google.com/images?um=1&hl=en&biw=1076&bih=569&tbs=isch:1&sa=1&q=lowess+graph&aq=f&aqi=&aql=&oq=&gs_rfai= – Jeromy Anglim Sep 07 '10 at 07:12
2

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.

Andrew Robinson
  • 739
  • 5
  • 9