R help - Jul 2007 - Piecewise Regression with a known slope

Hey, all.  I'm working on a data set with a broken stick linear  
regression where I know one of the two slopes.  It is a negative  
linear function until the line intersects with the x-axis, at which  
point it becomes 0.  It is not a nonlinear asymptotic function, and,  
indeed, using negative exponential or logistic types of fits as an  
approximation has tended to lead to an under or overestimation of  
values.  I am also very interested to know just what the breakpoint  
in the data is.

Essentially

if x<psi y = a + bx + error, where b is negative
else y=0+error

and I want to know the value of psi, as well as a and b.  The error  
is also most likely gamma distributed, as values<0 are not possible  
(nutrient concentrations).

I have attempted to use the segmented package, specifying something  
close to the visually estimated breakpoint as the value of psi, but  
it continues to return fits with a breakpoint that occurs somewhere  
in the middle of the linear portion of the line, and the slope and  
intercept of the second half of the regression is not 0.

Is there either a package that exists that will allow me to estimate  
such a model, or a function that I can use for optim or nlm (I admit,  
I am a rank novice at coding such functions).

Thanks so much!

-Jarrett

Dear Jarrett,

If I understand correctly your post, your constraints may be achieved
straightforwardly in segmented. See the code below.
> The error is also most likely gamma distributed..[SNIP]
The 'error' component can be specified in the 'initial' model by
means of
the family argument in the glm() function
> I have attempted to use the segmented package, specifying something
> close to the visually estimated breakpoint as the value of psi, but
> it continues to return fits with a breakpoint that occurs somewhere
> in the middle of the linear portion of the line, ..[SNIP]
Mhmm..you can specify different starting values and to assess the
differences in the 'final' estimates. However if you says "..it
continues
to return fits with a breakpoint that occurs somewhere in the middle.."
probably the ML estimate of the breakpoint is really in the middle.

Hope this helps you,

best,
vito



####simulate data
n<-50
x<-1:n/n
y<- 0-pmin(x-.5,0)+rnorm(50)*.03
plot(x,y) #This should be your scatterplot..
abline(0,0,lty=2)

o<-lm(y~x) #or glm(y~x,family=..)
o1<-segmented(o,seg.Z=~x,psi=list(x=.3))
slope(o1) #get the slope
points(x,fitted(o1),col=2)

#a parsimonious modelling: constrain right slope=0
o<-lm(y~1)
xx<- -x
o2<-segmented(o,seg.Z=~xx,psi=list(xx=-.3))
slope(o2)
points(x,fitted(o2),col=2)

#now constrain \hat{\mu}(x)=0 for x>psi

o<-lm(y~0)
o3<-segmented(o,seg.Z=~xx,psi=list(xx=-.3))
slope(o3)
points(x,fitted(o3),col=3)
> Hey, all.  I'm working on a data set with a broken stick linear
> regression where I know one of the two slopes.  It is a negative
> linear function until the line intersects with the x-axis, at which
> point it becomes 0.  It is not a nonlinear asymptotic function, and,
> indeed, using negative exponential or logistic types of fits as an
> approximation has tended to lead to an under or overestimation of
> values.  I am also very interested to know just what the breakpoint
> in the data is.
>
> Essentially
>
> if x<psi y = a + bx + error, where b is negative
> else y=0+error
>
> and I want to know the value of psi, as well as a and b.  The error
> is also most likely gamma distributed, as values<0 are not possible
> (nutrient concentrations).
>
> I have attempted to use the segmented package, specifying something
> close to the visually estimated breakpoint as the value of psi, but
> it continues to return fits with a breakpoint that occurs somewhere
> in the middle of the linear portion of the line, and the slope and
> intercept of the second half of the regression is not 0.
>
> Is there either a package that exists that will allow me to estimate
> such a model, or a function that I can use for optim or nlm (I admit,
> I am a rank novice at coding such functions).
>
> Thanks so much!
>
> -Jarrett
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>
>