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 ______________________________________________ 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.