Solved – How to calculate Theil–Sen estimator in R

estimationrregression

I found a useful method named "Theil–Sen estimator" on wikipedia: here

This method seems to be covered with mblm function on the mblm library.

I just downloaded it to do some tests.

I have one doub about it, if I do:

x <- c(1,2,3)
y <- c(2,6,7)
mod <- mblm(x~y)

It work perfectly, but I don't understand how to remove the intercept (set it to zero).

With lm() I do lm(x~y+0) but in this case I still see the intercept value.

So the first question is: How to set the intercept at 0 ?

The second question is about using this function with a matrix, take a look at the code below:

> y
[1] 6 4 1 8 7
> 
> t
     [,1] [,2]
[1,]    2   32
[2,]    5  762
[3,]    6    5
[4,]   23    9
[5,]   34   16
> mblm(y~t)

Call:
mblm(formula = y ~ t)

Coefficients:
(Intercept)            t  
         NA           NA

Why those NA / NA ?

If I do the same thing with lm() I don't have problems.

> lm(y~t)

Call:
lm(formula = y ~ t)

Coefficients:
(Intercept)           t1           t2  
  3.453e+00    1.247e-01    5.435e-06

So the second question is: Can I use a matrix on the right side of the formula?

Thank you

Best Answer

No, you can't use an intercept of zero. It's not clear what that would mean for this test. As the wikipedia page describes, this test computes the slope for all pairs of points and then reports the median of those slopes. If you know you want the fitted line to start at zero, this doesn't quite make sense. Maybe you'd want to compute the slopes of all points compared with that zero point only, but that wouldn't be the Theil-Sen estimate anymore; it would just be the sign test.

And no, you can't use a matrix on the right side; this test is for one predictor only. The slope between two pairs of points then has both a magnitude and a direction, and it's not clear how one would take the median of them.

Finally, the default for this function is to use Siegel repeated medians; for the Thiel-Sen test, you need to specify repeated=FALSE. See the help page for more information.

Related Question