[R] piece-wise linear regression nls function

Thu Jan 10 03:17:21 CET 2013

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}}