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:
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.
--
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612
http://dssm.unipa.it/vmuggeo
______________________________________________
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.
