[R] Kendall Theil line as fit?

[Peter Ehlers]

On 2011-03-11 17:10, jonbfish wrote:
> Thanks again. I guess I have been fitting in splus gui and in R Commander, not really calling it myself. I was thinking maybe there was a way to define it as a class or something.
>
> I will have to look at segments. Any thoughts if this might be easy to use with years on an axis?

Alright, I see what you're after. (It would have helped
if you had mentioned R Commander sooner.)

All you need is to use lines() with appropriate endpoints.
Here's how:

   x <- 1999:2010
   y <- c(17.942,3.43,14.062,14.814,13.778,13.706,
          9.748,13.088,12.1309728,9.644646671,9.134,8.84)
   mod <- lm(y ~ x)

   kt <- function(x, y){
     yy <- outer(y, y, "-")
     xx <- outer(x, x, "-")
     z  <- yy / xx
     s  <- z[lower.tri(z)]
     slope <- median(s)
     intercept <- median(y) - slope * median(x)
     list(intercept = intercept, slope = slope)
   }

   ind1 <- which.min(x)
   ind2 <- which.max(x)

   w <- kt(x, y)
   a <- w[["intercept"]]
   b <- w[["slope"]]

   x.ends <- c(x[ind1], x[ind2])
   y.ends.kt <- a + b * x.ends

   a <- coef(mod)[1]
   b <- coef(mod)[2]
   y.ends.lm <- a + b * x.ends

   plot(y ~ x)
   lines(x.ends, y.ends.kt, col = "blue")
   lines(x.ends, y.ends.lm, col = "red")

Peter Ehlers

>
> From: Peter Ehlers [via R] [mailto:ml-node+3349445-1733968405-216550 at n4.nabble.com]
> Sent: Friday, March 11, 2011 06:38 PM
> To: Anthony Seeman
> Subject: Re: Kendall Theil line as fit?
>
> On 2011-03-11 14:43, jonbfish wrote:
>
>> Thanks for the response, sorry I didn't post it initially.
>>
>> kt.mat<-
>> function(x,y,z){
>> for(i in 1:length(x)){for(j in
>> 1:length(y)){z[i,j]<-(y[j]-y[i])/(x[j]-x[i])}}
>> return(z)}
>>
>>
>> kt.slope<-
>> function(x,y,z,s){
>> count<-0
>> for(i in 1:length(x)){for(j in 1:length(y)){
>> if(j>= i+1) {
>> count<-count+1
>> s[count]<-z[i,j]}
>> }}
>> print(count)
>> return(s)}
>>
>> #Site23
>>
>> x<- c(1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010)
>> y<-
>> c(17.942,3.43,14.062,14.814,13.778,13.706,9.748,13.088,12.1309728,9.644646671,9.134,8.84)
>>
>> z<-matrix(0:0,length(x),length(y))
>> z<-kt.mat(x,y,z)
>> z
>>
>> s<-c(1:(length(x)*(length(x)-1)/2))
>> s<-kt.slope(x,y,z,s)
>> s
>> slope=median(s)
>> intercept=median(y)-slope*median(x)
>> cbind(slope,intercept)
>> plot(x,y)
>>
>> abline(intercept,slope)
>
> Okay, you're using abline() for the KT line.
> But I still don't know what you're after.
>   From your original post:
>
>     Is there a way to make it appear like a regression fit instead
>     of a line that extends from the edges of the plot? I would like
>     to have the OLS appear as a dotted line and the KT a solid line
>     but as it is the KT line is longer.
>
> So how are plotting your 'regression fit'?
> abline( lm( y ~ x ) ) would also extend across the plot.
> I suppose that you could use segments() with the
> range of x-values.
>
> BTW, here's a shorter version of your code:
>
>    yy<- outer(y, y, "-")
>    xx<- outer(x, x, "-")
>    z<- yy / xx
>    s<- z[lower.tri(z)]
>    slope<- median(s)
>    intercept<- median(y) - slope * median(x)
>
> Peter Ehlers
>