[R] Piecewise

[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:
> 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 at 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