[R] piece-wise linear regression nls function

[Rolf Turner]

Instead of reinventing the wheel, why not use the "segmented" package?
Having  stored your data in a data frame "X" I did:

require(segmented)
m1 <- lm(FM ~ BMIJS,data=X)
m2 <- segmented(m1,seg.Z=~BMIJS,psi=list(BMIJS=35))

which worked instantaneously.  The break point is estimated as 41.580, with
a standard error of 2.094  I then did:

with(X, plot(BMIJS,FM))
plot(m2,add=TRUE)

which seems to look as good as one can expect.

I must say however that the plot of your data does not look to me
as though a broken-stick model is appropriate.  Why not just a straight 
line?

     cheers,

         Rolf Turner

On 01/10/2013 02:33 PM, John Sorkin wrote:
> windows 7, R 2.12
>   
> I am trying to run a piecewise linear regression with a single knot, i.e. a regression composed of two straight lines where the two lines intersect at an x value given by the variable knot. I wish to estimate the slope of both lines, the value of knot, the x value where the two lines intersect, and an intercept. I am using the nls code below, and get the following error message:
> Error in nls(FM ~ blow * BMIJS + bhi * sapply(BMIJS - knot, max, 0), start = list(knot = 25,  :
>    singular gradient
> nls code:
>   
> test <- nls(FM~blow*BMIJS+bhi*sapply(BMIJS-knot,max,0),start=list(knot=25,blow=1,bhi=1),data=FeMaleData)
> summary(test)
>   
> greatly shortened version of my data (the full data set has 450 records)
>   
>   FM    BMIJS
> 2   55.878 40.57273
> 4   34.270 27.76939
> 5   20.123 21.73818
> 6   19.320 19.71203
> 9   49.701 43.55356
> 10  51.188 37.84742
> 11  46.753 37.71003
> 13  65.079 37.23438
> 14  37.097 36.81806
> 15  30.625 29.92783
> 17  50.617 42.42754
> 18  63.954 48.78709
> 20  29.790 26.97648
> 21  36.558 34.79373
> 22  41.275 33.03063
> 24  27.682 27.24508
> 26  37.968 35.41399
> 28  24.878 27.20250
> 30  47.513 35.77961
> 31  51.315 37.46032
> 33  41.944 36.40212
> 34  38.150 32.83818
> 35  60.719 42.48594
> 36  42.643 34.29355
> 38  40.728 32.42817
> 42  34.814 30.57573
> 43  32.896 29.32912
> 44  30.430 25.44183
> 46  48.986 37.90910
> 49  47.485 36.34642
> 52  46.312 38.64647
> 54  45.228 33.08783
> 55  45.391 35.86965
> 59  37.256 32.66507
> 60  27.367 28.49880
> 63  38.663 34.34131
> 64  34.527 29.57858
> 67  58.368 38.97266
> 68  13.473 17.35397
> 69  22.456 20.80958
> 71  28.829 25.50056
> 73  15.487 20.22202
> 76  18.313 21.38991
> 77  41.535 36.85707
> 78  56.124 40.51978
> 80  52.587 40.77256
> 81  24.991 25.48543
> 83  56.327 39.97214
> 84  70.836 36.52915
> 85  62.294 42.45244
> 86  39.689 35.18527
> 87  35.006 35.15136
> 88  47.378 37.54779
> 89  18.149 23.99236
> 90  33.041 28.10476
> 91  28.884 26.74443
> 92  37.670 32.25230
> 94  55.410 43.72364
> 99  34.461 35.05930
> 101 59.727 42.83035
> 102 41.913 35.64677
> 104 66.644 41.01642
> 105 55.250 43.86426
> 107 45.196 31.78370
> 108 36.476 33.45537
> 109 34.386 29.08402
> 110 39.277 36.98500
> 111 53.789 45.54654
> 112 33.077 29.09559
> 116 57.246 39.98031
> 120 52.546 40.12191
> 122 34.409 29.70977
> 123 31.188 28.75295
> 126 54.567 38.15226
> 129 19.193 22.71878
> 133 39.322 33.45712
> 134 41.415 31.28980
> 136 57.616 36.94016
> 140 28.162 24.40219
> 142 37.524 29.92673
> 143 29.611 29.15452
> 144 26.780 26.53462
> 146 47.219 35.14919
> 147 35.341 28.68955
> 148 44.827 37.68317
> 149 54.180 41.12226
> 150 41.636 30.00930
> 151 33.626 28.00164
> 156 34.334 29.64970
> 160 36.317 30.12031
> 161 46.823 35.64603
> 163 39.506 34.27740
> 164 61.619 39.20019
> 169 48.984 35.77558
> 171 66.467 41.59008
> 172 70.144 42.79996
> 173 37.324 31.56521
> 174 66.882 46.04938
> 182 54.239 38.21065
> 184 48.800 32.01630