Thanks, that is what I was trying to do. ________________________________ From: Peter Ehlers [via R] [mailto:ml-node+3350009-943873321-216...@n4.nabble.com] Sent: Saturday, March 12, 2011 4:27 AM To: Anthony Seeman Subject: Re: Kendall Theil line as fit?
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:[hidden > email]</user/SendEmail.jtp?type=node&node=3350009&i=0&by-user=t>] > 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 > ______________________________________________ [hidden email]</user/SendEmail.jtp?type=node&node=3350009&i=1&by-user=t> 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. ________________________________ If you reply to this email, your message will be added to the discussion below: http://r.789695.n4.nabble.com/Kendall-Theil-line-as-fit-tp3344617p3350009.html To unsubscribe from Kendall Theil line as fit?, click here<http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3344617&code=YXNlZW1hbkBpYXNveWJlYW5zLmNvbXwzMzQ0NjE3fDI0NTMyMDgxOQ==>. -- View this message in context: http://r.789695.n4.nabble.com/Kendall-Theil-line-as-fit-tp3344617p3354235.html Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]] ______________________________________________ R-help@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.