[R] Piecewise

[Ben Bolker]

Joe Waddell <joecwaddell <at> gmail.com> writes:

> 
> Hi,
> I am a biologist (relatively new to R) analyzing data which we predict 
> to fit a power function.  I was wondering if anyone knew a way to model 
> piecewise functions in R, where across a range of values (0-x) the data 
> is modeled as a power function, and across another range (x-inf) it is a 
> linear function.  This would be predicted by one of our hypotheses, and 
> we would like to find the AICs and weights for a piecewise function as 
> described, compared with a power function across the entire range.
> 
> I have looked into the polySpline function, however it appears to use 
> only polynomials, instead of the nls models I have been using.  Thanks 
> in advance for any help you might be able to offer.

 You should be able to do it in nls as follows ...
There are some potentially subtle issues if you get into
the details of likelihood profile confidence intervals
for the breakpoint (since the likelihood profile is flat 
between/discontinuous across the locations of data points),
but hopefully that won't bite you in your applications.

## function to generate piecewise power-law/linear "data"
f <- function(x,brk,alpha,beta,b,sd) {
  mu <- ifelse(x<brk,alpha*x^beta,(alpha*brk^beta)-b*(x-brk))
  rnorm(length(x),mean=mu,sd=sd)
}

## generate "data"
set.seed(1001)
x <- runif(1000,max=100)
y <- f(x,brk=50,alpha=100,beta=-0.5,b=1,sd=5)

## take a quick look, vs. theoretical curve
plot(y~x)
curve(f(x,50,100,-0.5,1,0),col=2,add=TRUE,n=1000,lwd=2)

## fit, using the "I()" command to do the piecewise part
dat <- data.frame(x,y)
fit1 <- nls(y~I(x<brk)*alpha*x^beta+I(x>brk)*((alpha*brk^beta)-b*(x-brk)),
            start=list(brk=60,alpha=110,beta=-0.75,b=2),data=dat)

## plot the fit
xvec <- seq(0,100,length=200)
lines(xvec,predict(fit1,newdata=data.frame(x=xvec)),col=4,lwd=2)
## testing ...
with(as.list(coef(fit1)),
     lines(xvec,f(xvec,brk,alpha,beta,b,sd=0),col=5,lwd=2))


  Ben Bolker