Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am dealing with linear data with outliers, some of which are at more the 5 standard deviations away from the estimated regression line. I'm looking for a linear regression technique that reduces the influence of these points.

So far what I did is to estimate the regression line with all the data, then discard the data point with very large squared residuals (say the top 10%) and repeated the regression without those points.

In the literature there are lots of possible approaches: least trimmed squares, quantile regression , m-estimators, etc. I really don't know which approach I should try, so I'm looking for suggestions. The important for me is that the chosen method should be fast because the robust regression will be computed at each step of an optimization routine. Thanks a lot!

share|improve this question
A method that you did not mention is the use of Student-$t$ errors with unknown degrees of freedom. However, this may not be as fast as you need. – user10525 Dec 19 '12 at 11:37
@Procrastinator: (It's easy to imagine a configuration of outliers where) this will not work. – user603 Dec 19 '12 at 12:06
@user603 That is true for any method, there is no Panacea ;). I was simply pointing out another method. +1 to your answer. – user10525 Dec 19 '12 at 12:18
@Procrastinator: I agree that all methods will fail for some rate of contamination. And 'failure' in this context can be defined quantitatively and empirically. But the idea is to still favour those methods that will fail only at higher rates of contamination. – user603 Dec 19 '12 at 12:24
Since this is being done repeatedly during an optimization routine, perhaps the data in the regression are (eventually) changing slowly. This suggests an algorithm adapted to your situation: start with some form of robust regression, but when taking small steps during the optimization, simply assume in the next step that any previous outlier will remain an outlier. Use OLS on the data, then check whether the presumptive outliers are still outlying. If not, restart with the robust procedure, but if so--which might happen often--you will have saved a lot of computation. – whuber Dec 19 '12 at 16:38
up vote 41 down vote accepted

If your data contains a single outlier, then it can be found reliably using the approach you suggest (without the iterations though). A formal approach to this is

Cook, R. Dennis (1979). Influential Observations in Linear Regression. Journal of the American Statistical Association (American Statistical Association) 74 (365): 169–174.

For finding more than one outlier, for many years, the leading method was the so-called $M$-estimation family of approach. This is a rather broad family of estimators that includes Huber's $M$ estimator of regression, Koenker's L1 regression as well as the approach proposed by Procastinator in his comment to your question. The $M$ estimators with convex $\rho$ functions have the advantage that they have about the same numerical complexity as a regular regression estimation. The big disadvantage is that they can only reliably find the outliers if:

  • the contamination rate of your sample is smaller than $\frac{1}{1+p}$ where $p$ is the number of design variables,
  • or if the outliers are not outlying in the design space (Ellis and Morgenthaler (1992)).

You can find good implementation of $M$ ($l_1$) estimates of regression in the robustbase (quantreg) R package.

If your data contains more than $\lfloor\frac{n}{p+1}\rfloor$ outlier potentially also outlying on the design space, then, finding them amounts to solving a combinatorial problem (equivalently the solution to an $M$ estimator with re-decending/non-convex $\rho$ function).

In the last 20 years (and specially last 10) a large body of fast and reliable outlier detection algorithms have been designed to approximately solve this combinatorial problem. These are now widely implemented in the most popular statistical packages (R, Matlab, SAS, STATA,...).

Nonetheless, the numerical complexity of finding outliers with these approaches is typically of order $O(2^p)$. Most algorithm can be used in practice for values of $p$ in the mid teens. Typically these algorithms are linear in $n$ (the number of observations) so the number of observation isn't an issue. A big advantage is that most of these algorithms are embarrassingly parallel. More recently, many approaches specifically designed for higher dimensional data have been proposed.

Given that you did not specify $p$ in your question, I will list some references for the case $p<20$. Here are some papers that explain this in greater details in these series of review articles:

Rousseeuw, P. J. and van Zomeren B.C. (1990). Unmasking Multivariate Outliers and Leverage Points. Journal of the American Statistical Association, Vol. 85, No. 411, pp. 633-639.

Rousseeuw, P.J. and Van Driessen, K. (2006). Computing LTS Regression for Large Data Sets. Data Mining and Knowledge Discovery archive Volume 12 Issue 1, Pages 29 - 45.

Hubert, M., Rousseeuw, P.J. and Van Aelst, S. (2008). High-Breakdown Robust Multivariate Methods. Statistical Science, Vol. 23, No. 1, 92–119

Ellis S. P. and Morgenthaler S. (1992). Leverage and Breakdown in L1 Regression. Journal of the American Statistical Association, Vol. 87, No. 417, pp. 143-148

A recent reference book on the problem of outlier identification is:

Maronna R. A., Martin R. D. and Yohai V. J. (2006). Robust Statistics: Theory and Methods. Wiley, New York.

These (and many other variations of these) methods are implemented (among other) in the robustbase R package.

share|improve this answer
Now that is a great answer! – Peter Flom Dec 19 '12 at 12:43
Thanks a lot user603! In my problem $p < 10$ and there are no outliers in the design space (because the explanatory variables are simulated from a normal distribution). So maybe I can try with the m-estimator? In any case all the other references you have given me will be very useful once I will start working on more complex applications ($p$ >> 10) of my algorithm. – Matteo Fasiolo Dec 19 '12 at 13:40
@Jugurtha: In that case (no outlier in the design space and $p<10$) $M$ estimators are indeed the preferred solution. Consider the '' function in the robustbase package, the 'rlm' function in the MASS package or the l1 regression in the quantreg package. I would still also run the LTS-regression in a few case and compare the results, since they can withstand more outliers. I would do this just as a check of whether the contamination rate is not higher than you suspect. – user603 Dec 19 '12 at 13:47
"A big advantage is that most of these algorithms are embarrassingly parallel." I like the wording. ;) – Mateen Ulhaq Apr 23 at 0:38
@Mateen, well, it is the term of art after all. :) – J. M. 19 hours ago

For simple regression (single x), there's a lot to be said for the Theil-Sen line in terms of robustness to y-outliers and to influential points as well as generally good efficiency (at the normal) compared to LS for the slope. The breakdown point for the slope is nearly 30%; as long as the intercept (there are a variety of possible intercepts people have used) doesn't have a lower breakdown, the whole procedure copes with a sizable fraction of contamination quite well.

Its speed might sound like it would be bad - median of $\binom{n}{2}$ slopes looks to be $O(n^2)$ even with an $O(n)$ median - but my recollection is that it can be done more quickly if speed is really an issue ($O(n \log n)$, I believe)

Edit: user603 asked for an advantage of Theil regression over L1 regression (they have the same breakdown point on y-outliers). The answer is the other thing I mentioned - influential points:


The red line is the $L_1$ fit (from the function rq in the quantreg package). The green is a fit with a Theil slope. All it takes is a single typo in the x-value - like typing 533 instead of 53 - and this sort of thing happens. So the $L_1$ fit isn't robust to a single typo in the x-space.

share|improve this answer
it can indeed be computed in time $n\log n$. Could you elaborate on what advantage (in the single x case) the T-S estimator has over, say, $l_1$ regression? – user603 Mar 10 '13 at 9:39
@user603 See the edit. – Glen_b Mar 10 '13 at 11:28
(+1) thanks for the edit. It's important to point this feature out. – user603 Mar 11 '13 at 9:27
And what's the advantage over an MM-estimate, such as lmrob() from R package robustbase or even {no need to install anything but 'base R'} rlm(*, ... method="MM") from package MASS? These have full breakdown point (~ 50%) and are probably even more efficient at the normal. – Martin Mächler May 17 '13 at 12:41
@MartinMächler It seems like you're arguing against a claim I haven't made there. If you'd like to put up an answer which also contains a comparison of other high-breakdown robust estimators, especially ones that are roughly as simple to understand for someone at the level of the OP, I'd look forward to reading it. – Glen_b May 17 '13 at 23:05

Have you looked at RANSAC (Wikipedia)?

This should be good at computing a reasonable linear model even when there are a lot of outliers and noise, as it is built on the assumption that only part of the data will actually belong to the mechanism.

share|improve this answer
yea but adding a simple re-weighting step yields an estimator (LTS) that is equally robust and so much more stable and statistically efficient. Why not do? – user603 Apr 29 '13 at 1:13

I found the $l_1$ penalized error regression best. You can also use it iteratively and reweight samples, which are not very consistent with the solution. The basic idea is to augment your model with errors: $$y=Ax+e$$ where $e$ is the unknown error vector. Now you perform the regression on $$\parallel y-Ax-e \parallel_2^2+ \lambda \parallel e \parallel_1$$. Interestingly you can of course use "fused lasso" for this when you can estimate the certainty of your measurements in advance and put this as weighting into $$W=diag(w_i)$$ and to solve the new slighty different task $$\parallel y-Ax-e \parallel_2^2 + \lambda \parallel W e \parallel_1$$

More information can be found here:

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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