[R] piecewise linear approximation

[Rolf Turner]

On 30/08/2007, at 12:45 PM, Naxerova, Kamila wrote:

> Dear list,
>
> I have a series of data points which I want to approximate with  
> exactly two
> linear functions. I would like to choose the intervals so that the  
> total
> deviation from my fitted lines is minimal. How do I best do this?
>
> Thanks!
> Kamila

If you want to minimize the total ***squared*** deviations, then the  
code that
I enclose below should do the trick for you.  If you really want to  
minimize
the total (presumably absolute) deviation, then you'll have to adapt the
code a bit.

The argument ``b'' to the bsr() function is a vector of candidate  
values for
the ``join point'', e.g. b <- seq(lo,hi,length=100), where lo and hi are
suitably chosen lower and upper bounds on the location of the join  
point.

				cheers,

					Rolf Turner


bsr <- function(x,y,b) {
#
# Broken stick regression.
#
if(length(b)==1) {
         x1 <- ifelse(x<=b,0,x-b)
         temp <- lm(y~x+x1)
         ss <- sum(resid(temp)^2)
         rslt <- list(fit=temp,ss=ss,b=b)
         class(rslt) <- "bsr"
         return(rslt)
}
temp <- list()
for(v in b) {
         x1 <- ifelse(x<=v,0,x-v)
         fff <- lm(y~x+x1)
         temp[[paste(v)]] <- sum(resid(fff)^2)
}
ss <- unname(unlist(temp))
v <- b[which.min(ss)]
x1 <- ifelse(x<=v,0,x-v)
rslt <- list(fit=lm(y~x+x1),ss=unname(unlist(temp)),
              b=b,ss.min=min(ss),b.min=v)
class(rslt) <- "bsr"
rslt
}


===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+===+=== 
+===+===

plot.bsr <- function(x,cont=TRUE,add=FALSE,rr=NULL,npts=101,...) {
#
# Plot method for broken stick regressions.
#
fit <- x
if(cont) {
         if(is.null(rr)) {
                 x <- fit$fit$model$x
                 rr <- range(x)
         }
         xx <- seq(rr[1],rr[2],length=npts)
         b  <- fit$b.min
         if(is.null(b)) b <- fit$b
         x1 <- ifelse(xx<=b,0,xx-b)
         yh <- predict(fit$fit,newdata=data.frame(x=xx,x1=x1))
} else {
         xx <- fit$fit$model[,2]
         yh <- predict(fit$fit)
}
if(add) lines(x=xx,y=yh,...) else plot(x=xx,y=yh,...)
}


######################################################################
Attention:\ This e-mail message is privileged and confidenti...{{dropped}}