7
$\begingroup$

I've just written some C code for Theil-Sen, after some Googling (I don't have any definitive documentation on it). My understanding of the intercept calculation is that I first calculate the median slope, and then construct a line through every data point with this slope, find the intercept of every line, and then take the median intercept.

The only way I can find to test the code is to compare the results against the Kendall-Theil Robust Line program, from the USGS. On a dataset of 237 points (healthcare data, with a Pearson correlation of ~0.55), we agree exactly on the median slope, but disagree on the intercept (by 1.4%). According to my figures, the KTRL intercept isn't the median intercept, but is instead 46% of the way through the range.

After some digging around in the KTRL code, it appears that they calculate the intercept by creating a single 'median line', rather than the median of all intercepts. Their intercept is medianY - medianX * median slope.

Any feedback on which is the "right" way to do this, if there is one, or how this is handled in R/etc?

Thanks.

$\endgroup$
  • $\begingroup$ Isn't the intercept usually the median of the residuals $y_i-\beta x_j$? $\endgroup$ – Dimitriy V. Masterov Feb 22 '13 at 0:34
  • $\begingroup$ @DimitriyV.Masterov: no! that would not be regression equivariant! $\endgroup$ – user603 Feb 22 '13 at 2:59
  • 3
    $\begingroup$ @user603 There are papers (possibly Sen's paper) which discuss what's it's equivariant to and what it's not. This is the most common form I have seen. $\endgroup$ – Glen_b -Reinstate Monica Feb 22 '13 at 3:05
8
$\begingroup$

The Theil-Sen estimator is essentially an estimator for the slope alone; the line has been constructed in a host of different ways.

In fact there are a large variety of ways to calculate the intercept.

You said:

My understanding of the intercept calculation is that I first calculate the median slope, and then construct a line through every data point with this slope, find the intercept of every line, and then take the median intercept.

That's one way that it's done.

There are a couple of approaches that compute kinds of weighted-medians but based off that same approach (putting more weight on the points further apart in x-space).

Another common one (probably the most common) is to compute median($y-bx$). This is what Sen looked at, for example.

Another is to try to get an estimator with higher efficiency at the normal (akin to that of the slope estimator in typical situations) and similar breakdown point to the slope estimate (there's probably little point in having better breakdown at the expense of efficiency), such as using the Hodges-Lehmann estimator (median of pairwise averages) on $y-bx$. This has a kind of symmetry in the way the slopes and intercepts are defined ... and generally gives something very close to the LS line when the normal assumptions nearly hold, whereas the Sen-intercept can be - relatively speaking - quite different.

And some people just compute the mean.

There are still other suggestions that have been looked at. There's really no 'one' intercept to go with the slope estimate.

Dietz lists several possibilities, possibly even including all the ones I mentioned, but that's by no means exhaustive.

$\endgroup$
  • $\begingroup$ +1. @OP: you really gotta wonder first if the ST estimator is really what you want. Google search for FastLTS or FastS too $\endgroup$ – user603 Feb 22 '13 at 3:02
  • 1
    $\begingroup$ @Glen_b - thanks - on a quick scan, the Dietz paper looks very good. She looks at 10 intercept methods, and mine (+Wikipedia/etc) is ${\hat{\alpha }}_{M}$. The KTRL code appears to use the Conover estimator, ${\hat{\alpha }}_{C}$ (it's in Visual Basic, largely uncommented, and uncited, so it's difficult to see what's going on). On the error distributions that Dietz uses, Conover always performs better than ${\hat{\alpha }}_{M}$. $\endgroup$ – EML Feb 22 '13 at 10:15
  • $\begingroup$ @user603: probably not, but this is an unquantifiable performance metric on groups of humans, with unexplainable outliers, and the major estimator selection criteria is that they should be able to find it easily on the web, and get a warm feeling... :) $\endgroup$ – EML Feb 22 '13 at 10:21
  • $\begingroup$ @Glen_b - what are what you consider the current viable alternatives for TS estimation? $\endgroup$ – EngrStudent - Reinstate Monica Aug 27 '14 at 15:56
  • 1
    $\begingroup$ @EngrStudent The kinds of robust lines offered by typical robust regression routines are generally quite good but aren't quite as conceptually simple. They're easy to fit in (say) R - see ltsreg in MASS or the packages robust and robustbase for example. Which kind of estimator you need depends on what you need robustness against, how much of it you might need to deal with and how much you worry much about efficiency at some ideal model. [For example, if x-outliers are not an issue, quantile regression or M-estimation may do quite well, but if they are an issue, I wouldn't use those.] $\endgroup$ – Glen_b -Reinstate Monica Aug 27 '14 at 20:39
1
$\begingroup$

The Kendall–Theil Robust Line program from the USGS has a companion PDF.

On page 8 (PDF page 15) it states the method used and formula as you found but gives the reference as Conover.

Intercept

The estimate of the intercept is calculated by use of the Conover (1980) equation

$$b = Y_\text{median} -m\times X_\text{median}\ \ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6)$$

where
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ b\ \ \ \ \ $ is the estimated intercept,
$\ \ \ \ \ \ \ Y_\text{median}\ \ \ \ $ is the median of the response variables,
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ m\ \ \ \ $ is the estimated slope,
and
$\ \ \ \ \ \ X_\text{median}\ \ \ \ $ is the median of the explanatory variables.

I confirm this does produce the same result as the program. Whether there are superior methods and so on, as always, a matter of opinion and your particular circumstances.

The M-estimation algorithm is arguably erroneous.

for i = 1,   # dat-1 do
for j = i+1, # dat do

Change there is j indexing from i+1 and don't process instances of i == j.

Then either rank (sort) the result choosing index as described, or arguably take the median, which will give a slightly different result. If you plot, the data will look like a CDF plot.

$\endgroup$
  • $\begingroup$ The intercept is calculated so that the line will run through the median of input data, would be the self-consistent non-parametric solution. $\endgroup$ – Carl Aug 24 '16 at 0:12

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.