R help - Jan 2013 - piece-wise linear regression nls function

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  Thanks,John 
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the sole use of the
intended recipient(s) and may contain confidential and privileged information. 
Any unauthorized use, disclosure or distribution is prohibited.  If you are not
the intended recipient, please contact the sender by reply email and destroy all
copies of the original message.

On Jan 9, 2013, at 5: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)
I was surprised to see `sapply` inside a formula expression. I instead imagined
that this might have been what was meant:
> test <- nls( FM ~ blow*BMIJS + bhi*pmax(BMIJS-knot,0) ,
             
start=list(knot=25,blow=1,bhi=1),data=FeMaleData)
> summary(test)
Formula: FM ~ blow * BMIJS + bhi * pmax(BMIJS - knot, 0)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
knot  21.4960     3.2095   6.698 1.39e-09 ***
blow   0.8983     0.1264   7.106 2.02e-10 ***
bhi    0.9551     0.1610   5.931 4.63e-08 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

Residual standard error: 5.638 on 97 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 8.684e-09 

I offer not particular opinion on whether this is sensible, only htat it does
not break the interpreter's understanding of function application. and the
know seems within the range of the values, albeit to the left hand edge:
> with(FeMaleData, plot(FM~BMIJS) )
> lines(seq(15, 50), predict(test, newdata=list(BMIJS=seq(15, 50)) ) )

-- 
David.
> 
> 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  Thanks,John 
> John David Sorkin M.D., Ph.D.
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
> Confidentiality Statement:
> This email message, including any attachments, is for ...{{dropped:15}}

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