[R] Fitting curves with discontinuous derivatives

Handegard, Nils Olav nils.olav.handegard at imr.no
Thu Jan 8 14:29:06 CET 2004


Dear Group Readers,

I have a response variable Y (fish swimming speed) and three explanatory
variables X1 (distance to ship propeller), X2 (distance to fishing gear
wire) and Depth. Y=f(X1) initially looks like a straight line, then an
increase to a higher level at X1=100, and then a decrease after X1=-800.
Y=f(X2) looks similar, but the increase to a higher level at X2=100 is more
abrupt, and the decrease start earlier, around X2=0. I explain this as the
reaction occurs closer and is stronger to X2 as opposed to X1. Most likely
both variables are important. I would like to assess this "importance". 

The idea was to use GAM to resolve the dependence, but since X1 - X2 is
small (i.e. plot(X1,X2) is along a narrow band around X1=X2) this did not
work. I must therefore work with the variables separately.

I would like to try a continuous model defined as:
       { k_1                  , x<a
f(x) = { \sum_i^4 k_{i+1} x^i , a<=x<b
       { k_6                  , x>=b,
with continuity as an additional requirement. 

1. The first problem is that a and b are not given, and the derivatives with
respect to these are not defined. This can be overcome by searching the
definition space using optimise(), as outlined in:
http://maths.newcastle.edu.au/~rking/R/help/01c/3175.html. However,
suggestions for a neater solution are welcomed.

2. How can I specify this model. If the model were lines, I guess something
like
lm(y~bs(1:20, knots=c(5,11), degree=1))
would do, and if the lines were of similar degree, splines could be used.
But is it possible to define splines with different order between the knots,
and continuity as only requirement?

3. The idea is to bootstrap the estimates of a and b to find the start point
of the reaction. However, using this parameterised curve may be a bit ad
hoc. Other suggestions are warmly appreciated. I have also tried nls with a
logistic function:
est <- try(nls(y~cbind(exp((x-beta)*alpha)/
	(1+exp((x-beta)*alpha)) ), start=c(beta=Beta,alpha=Alpha),
      algorithm="plinear",contr=ctrl), silent=T)
but the fit were bad.

Sincerely,
Nils Olav Handegard
PhD Student

Please send any responses directly to me.


lm




More information about the R-help mailing list