From: Peter Ehlers <ehlers_at_ucalgary.ca>

Date: Sat, 12 Mar 2011 02:25:44 -0800

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 12 Mar 2011 - 10:31:38 GMT

Date: Sat, 12 Mar 2011 02:25:44 -0800

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@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
**>
*

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 12 Mar 2011 - 10:31:38 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Mon 14 Mar 2011 - 16:30:21 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*