Skip to content

Piecewise

3 messages · Joe Waddell, Ben Bolker, vito muggeo

Joe Waddell
# 
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.

Cheers,
-Joe Waddell
Ben Bolker
# 
Joe Waddell <joecwaddell <at> gmail.com> writes:
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
vito muggeo
# 
Hi,
In addition to useful Ben's suggestion, you have a, possibly simpler, 
alternative.

If you are willing to assume to know the power of you piecewise 
polynomial (beta parameter according to code of Ben) you can use the 
package segmented. Using the data generated by Ben in his previous 
email, you get

olm<-lm(y~0+I(x^(-.5))) #note that 0 as the true model has not interc
os1<-segmented(olm,seg.Z=~x,psi=40)

Results are comparable with the ones returned by the nls() (even as 
regard the St.Err of the breakpoint), although segmented returns a 
somewhat smaller residual sum of squares :-)

***segmented() vs. nls(): some possibles differences***

1)In segmented, you are assuming to know the power of polynomial. In 
practice you can try several different values such as {-1,-.5,1,2,3}, 
say, and to assess the fit.. However the St.Err from the other estimates 
might be underestimated..

2)In segmented, you need to specify starting value only for the breakpoint.

3)In segmented, you can easy generalize the model: multiple breakpoints 
and/or additional "linear" covariates and/or non-continuous response 
and/or non-identity link functions (e.g. gamma with log-link)...


Hope this helps you,

vito



Ben Bolker ha scritto: