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