1

I am trying to write a basic function to add some lines of best fit to plots using nls. This works fine unless the data just happens to be defined exactly by the formula passed to nls. I'm aware of the issues and that this is documented behaviour as reported here.

My question though is how can I get around this and force a line of best fit to be plotted regardless of the data exactly being described by the model? Is there a way to detect the data matches exactly and plot the perfectly fitting curve? My current dodgy solution is:

#test data
x <- 1:10
y <- x^2
plot(x, y, pch=20)

# polynomial line of best fit
f <- function(x,a,b,d) {(a*x^2) + (b*x) + d}
fit <- nls(y ~ f(x,a,b,d), start = c(a=1, b=1, d=1)) 
co <- coef(fit)
curve(f(x, a=co[1], b=co[2], d=co[3]), add = TRUE, col="red", lwd=2)

Which fails with the error:

Error in nls(y ~ f(x, a, b, d), start = c(a = 1, b = 1, d = 1)) : 
  singular gradient

The easy fix I apply is to jitter the data slightly, but this seems a bit destructive and hackish.

# the above code works after doing...
y <- jitter(x^2)

Is there a better way?

Community
  • 1
  • 1
thelatemail
  • 91,185
  • 12
  • 128
  • 188
  • All the same, in the real world, this situation will never happen. There's always measurement error. Unless you're a teacher who's giving an R exam and devilishly supplying your students with perfect datasets, that is :-) . – Carl Witthoft Dec 19 '12 at 12:38
  • @CarlWitthoft I've had this kind of problems with data exported from Excel as CSV (and thereby rounded to visible digits) if n was small. – Roland Dec 19 '12 at 14:27
  • @Roland, I suppose if you round your test data inappropriately (i.e. losing valid sig figs), you get what you deserve :-) – Carl Witthoft Dec 19 '12 at 15:13
  • @CarlWitthoft Well, if you do rounding appropriate for measurement precision and then do regression with n=4 (which shouldn't be done, but such is life) ... – Roland Dec 19 '12 at 15:18

1 Answers1

6

Use Levenberg-Marquardt.

x <- 1:10
y <- x^2

f <- function(x,a,b,d) {(a*x^2) + (b*x) + d}
fit <- nls(y ~ f(x,a,b,d), start = c(a=1, b=0, d=0)) 

Error in nls(y ~ f(x, a, b, d), start = c(a = 1, b = 0, d = 0)) : 
  number of iterations exceeded maximum of 50

library(minpack.lm)
fit <- nlsLM(y ~ f(x,a,b,d), start = c(a=1, b=0, d=0))
summary(fit)

Formula: y ~ f(x, a, b, d)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a        1          0     Inf   <2e-16 ***
b        0          0      NA       NA    
d        0          0      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.49e-08

Note that I had to adjust the starting values and the result is sensitive to starting values.

fit <- nlsLM(y ~ f(x,a,b,d), start = c(a=1, b=0.1, d=0.1))

Parameters:
    Estimate Std. Error    t value Pr(>|t|)    
a  1.000e+00  2.083e-09  4.800e+08  < 2e-16 ***
b -7.693e-08  1.491e-08 -5.160e+00  0.00131 ** 
d  1.450e-07  1.412e-08  1.027e+01  1.8e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.191e-08 on 7 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 1.49e-08
Roland
  • 127,288
  • 10
  • 191
  • 288