# Robust monotonic regression in R

I have the following table in R

df <- structure(list(x = structure(c(12458, 12633, 12692, 12830, 13369,
13455, 13458, 13515), class = "Date"), y = c(6080, 6949, 7076,
7818, 0, 0, 10765, 11153)), .Names = c("x", "y"), row.names = c("1",
"2", "3", "4", "5", "6", "8", "9"), class = "data.frame")

> df
x     y
1 2004-02-10  6080
2 2004-08-03  6949
3 2004-10-01  7076
4 2005-02-16  7818
5 2006-08-09     0
6 2006-11-03     0
8 2006-11-06 10765
9 2007-01-02 11153


I can plot the points and a Tukey's linear fitting (line function in R) via

plot(data=df,  y ~ x)
lines(df$x, line(df$x, df$y)$fitted.values)


which produces:

All fine. The above plot shows energy consumption values, expected only to increase, so I'm happy with the fit not passing through those two points (which will be subsequently flagged as outliers).

However, "just" removing the last point and replot again

df <- df[-nrow(df),]
plot(data=df,  y ~ x)
lines(df$x, line(df$x, df$)$fitted.values)


The result is completely different.

My need is to have ideally the same result in both scenarios above. R doesn't seem to have ready to use function for monotonic regression, besides isoreg which however is piecewise constant.

EDIT:

As @Glen_b pointed out the outliers-to-sample size ratio is too big (~28%) for the regression technique used above. However, I believe there might be something else to consider. If I add the points at the beginning of the table:

df <- rbind(data.frame(x=c(as.Date("2003-10-01"), as.Date("2003-12-01")), y=c(5253,5853)), df)


and recalculate again like above plot(data=df, y ~ x); lines(df$x, line(df$x,df$y)$fitted.values) I get the same result, with a ration of ~22%

• Could you explain to us what you mean by "Tukey's line"? (He used various resistant line fitting methods.) – whuber Jan 15 '14 at 15:27
• @whuber oh I see, sorry. Is the method implemented in the R function line. You can have more details by typing ?line in the r console – Michele Jan 15 '14 at 15:30
• Thanks, but I'm afraid that's no good at all: the help merely refers to Tukey's 1977 EDA book--with which I am quite familiar and in which I can identify many line-fitting methods--and the code simply invokes a C program. Maybe we could make progress if you could explain more clearly what you are trying to achieve. How would you characterize (in general) the difference between your two "scenarios"? Why do you prefer the first solution? – whuber Jan 15 '14 at 15:33
• (+1) "Should only increase" is key: you are asking about how to perform (robust) monotonic regression. It would help to emphasize that point more in your question: you will get better answers. – whuber Jan 15 '14 at 15:51
• @Michele Maybe you could have a look at the nnls package (non-negative least squares). That should help you with the positivity constraints, but not with the outliers. – Matteo Fasiolo Jan 16 '14 at 11:14

I note that after you delete the last point, you only have seven values of which two (28.6%!) are outliers. Many robust methods don't have a breakdown point quite that high (e.g., Theil regression breaks down right at that point for n=7, though at large $$n$$ it goes to 29.3%), but if you must have such a high breakdown that it can manage so many outliers, you need to choose some approach that actually has that higher breakdown point.

There are some available in R; the rlm function in MASS (M-estimation) should deal with this particular case (it has high breakdown against y-outliers), but it won't have robustness to influential outliers.

Function lqs in the same package should deal with influential outliers, or there are a number of good packages for robust regression on CRAN.

You may find Fox and Weisberg's Robust Regression in R (pdf) a useful resource on several robust regression concepts.

All this is just dealing with robust linear regression and is ignoring the monotonicity constraint, but I imagine that will be less of a problem if you get the breakdown issue sorted. If you're still getting negative slope after performing high-breakdown robust regression, but want a nondecreasing line, you would set the line to have slope zero - i.e. choose a robust location estimate and set the line to be constant there. (If you want robust nonlinear-but-monotonic regression, you should mention that specifically.)

In response to the edit:

You seem to have interpreted my example of the Theil regression as a comment about the breakdown point of line. It was not; it was simply the first example of a robust line that came to me which broke down at a smaller proportion of contamination.

As whuber already explained, we can't easily tell which of several lines is being used by line. The reason why line breaks down as it does depends on which of several possible robust estimators that Tukey mentions and line might use.

For example, if it's the line that goes 'split data into three groups and for the slope use the slope of the line joining the medians of the outer two thirds' (sometimes called the three-group resistant line, or the median-median line), then its breakdown point is asymptotically 1/6, and its behavior in small samples depends on exactly how the points are allocated to the groups when $$n$$ is not a multiple of 3.

Please note that I am not saying it is the three group resistant line that is implemented in line - in fact I think it isn't - but simply that whatever they have implemented in line may well have a breakdown point such that the resulting line cannot deal with 2 odd points out of 8, if they're in the 'right' positions.

In fact the line that is implemented in line has some bizarre behavior - so odd that I wonder if it might have a bug - if you do this:

 x = y = 1:9 #all points lie on a line with slope 1
plot(x,y)
abline(line(x,y),col=2)


Then the line line has slope 1.2:

Off the top of my head, I don't recall any of Tukey's lines having that behavior.

• Hi thanks for all of that. At the moment my need is linear monotonic regression. The breakdown point is very interesting and I definitely should have considered. However, can you please comment the third chart I have just added? I add two point at the beginning which should bring the ratio to 22% – Michele Jan 16 '14 at 10:46
• by the way lqs does the job! So I accept your answer :-) thanks a lot. If you could still help me out in understanding the third chart would be great! Cheers – Michele Jan 16 '14 at 11:24
• I have made an edit that I hope clarifies things somewhat. If I've missed something, let me know. – Glen_b -Reinstate Monica Jan 16 '14 at 22:51
• Many thanks. I'm not a statistician (easy to tell), and you've been very helpful! I think I'll stick with lqs for now. Thanks again – Michele Jan 17 '14 at 10:22