[R] nls for piecewise linear regression not converging to least square
Thomas Lumley
tlumley at u.washington.edu
Mon Apr 19 23:40:45 CEST 2010
On Mon, 19 Apr 2010, Karen Chang Liu wrote:
> Hi R experts,
>
> I'm trying to use nls() for a piecewise linear regression with the first
> slope constrained to 0. There are 10 data points and when it does converge
> the second slope is almost always over estimated for some reason. I have
> many sets of these 10-point datasets that I need to do. The following
> segment of code is an example, and sorry for the overly precise numbers,
> they are just copied from real data.
>
> y1<-c(2.37700445, 1.76209775, 0.09795576, 2.21834963, 6.62262243,
> 15.70471269, 21.92956392, 36.39401717, 32.43620195, 44.77442277)
> x1<-c(24.6, 28.9, 33.2, 37.6, 42.0, 46.4, 50.9, 55.3, 59.8, 64.3)
>
> dat <- data.frame(x1,y1)
> nlmod <- nls(y1 ~ ifelse(x1 < xint+(yint/slp), yint, yint +
> (x1-(xint+(yint/slp)))*slp),
> data=dat, control=list(minFactor=1e-5,maxiter=500,warnOnly=T),
> start=list(xint=39.27464924, yint=0.09795576, slp=2.15061064),
> na.action=na.omit, trace=T)
>
> ##plotting the function
> plot(dat$x1,dat$y1)
> segments(x0=0, x1=coef(nlmod)[1]+coef(nlmod)[2]*coef(nlmod)[3],
> y0=coef(nlmod)[2], y1=coef(nlmod)[2])
> segments(x0=coef(nlmod)[1]+coef(nlmod)[2]*coef(nlmod)[3],x1=80,
> y0=coef(nlmod)[2], y1=80*coef(nlmod)[3]+coef(nlmod)[2])
>
> As you can see from the plot, the line is above all data points on the
> second segment. This seems to be the case for different datasets. I'm
> wondering if anyone can help me understand why this happens. Is this because
> there are too few data points or is it because the likelihood function is
> just not smooth enough?
>
I think there's something wrong with your graph. If I do
points(x1,fitted(nlmod),col="red")
I get points that are on the horizontal line segment, but then go through the data nicely on the right.
-thomas
Thomas Lumley Assoc. Professor, Biostatistics
tlumley at u.washington.edu University of Washington, Seattle
More information about the R-help
mailing list