0

I'm trying to fit parabola to my data in R. I already created the fit in MATLAB (so I know the parameter values), but when I run this code:

library("nls2")
fo <- y ~ A * (x ^ 2) + B * x + C
sA <- seq(-.000001,0,len=10)
sB <- seq(0,.001,len=10)
sC <- seq(1,20,len=10)

st1 <- expand.grid(A=sA,B=sB,C=sC)
mod1 <- nls2(fo,start=st1,algorithm="brute-force")
z <- nls(fo,start=coef(mod1))

x <- c(50072,50072,52536,52536,53768,53768,53768,54384,54384,54384, 54384,54692,54692,54846,54846,54846,54846,54846,54923,54923, 54923,54923,54923,54923,54961.5,54961.5,54961.5,54961.5,54961.5,54961.5,55000,55‌​000,55000) 
y <- c(0.46007,0.47185,0.24377,0.28506,0.16560,0.13587,0.16844, 0.08156,0.06917,0.09357,0.08704,0.02968,0.03293,0.01233,0.01509, 0.01095,0.005,0.00704,0.00861,0.00287,0.00456,0.01043,0.00373,0.00594, 0.00628,0.00089,0.00085,0.00055,0.00142,0.0022,-0.00091,-0.00209,-0.00293)

I get the following: z

           A             B             C 
-7.728222e-09  7.177715e-04 -1.610196e+01 
> z
Nonlinear regression model
  model: y ~ A * (x^2) + B * x + C
   data: parent.frame()
         A          B          C 
-7.728e-09  7.178e-04 -1.610e+01 
 residual sum-of-squares: 0.0106

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.457e-07

But confint(z) gives me the following error message:

Waiting for profiling to be done...
Error in approx(sp$y, sp$x, xout = cutoff) : 
  need at least two non-NA values to interpolate

What went wrong ? Thanks !

IRTFM
  • 258,963
  • 21
  • 364
  • 487
Dusan Kojic
  • 339
  • 4
  • 11
  • 1
    Why do you use `nls` for a linear model? – Roland Jan 04 '15 at 16:20
  • Also, which package is `confine` from? Please provide a [reproducible example](http://stackoverflow.com/a/5963610/1412059). – Roland Jan 04 '15 at 16:22
  • sorry, my "confint(z)" got autocorrected into "confine" – Dusan Kojic Jan 04 '15 at 16:25
  • How can I do this as a linear model ? – Dusan Kojic Jan 04 '15 at 16:30
  • here's a sample of my data: x<-c(50072,50072,52536,52536,53768,53768,53768,54384,54384,54384, 54384,54692,54692,54846,54846,54846,54846,54846,54923,54923, 54923,54923,54923,54923,54961.5,54961.5,54961.5,54961.5,54961.5,54961.5,55000,55000,55000) y <- c(0.46007,0.47185,0.24377,0.28506,0.16560,0.13587,0.16844, 0.08156,0.06917,0.09357,0.08704,0.02968,0.03293,0.01233,0.01509, 0.01095,0.005,0.00704,0.00861,0.00287,0.00456,0.01043,0.00373,0.00594, 0.00628,0.00089,0.00085,0.00055,0.00142,0.0022,-0.00091,-0.00209,-0.00293) – Dusan Kojic Jan 04 '15 at 16:47

1 Answers1

2

Since your model is linear, you should use lm:

x <- 1:100
set.seed(42)
y <- 20 * (x - 30)^2 + 4 + rnorm(100)

fit <- lm(y ~ poly(x, degree = 2, raw = TRUE))
confint(fit)
#                                       2.5 %     97.5 %
#(Intercept)                      18003.63406 18004.9075
#poly(x, degree = 2, raw = TRUE)1 -1200.04169 -1199.9835
#poly(x, degree = 2, raw = TRUE)2    19.99984    20.0004

In contrast to nls this doesn't use a numeric optimizer, but calculates an "exact" solution.

FYI: A model is considered non-linear, if the nonlinearity is in the parameters. y ~ a * x^2 + b * x + c is linear (in the parameters), but, e.g., y ~ sin(a - x) is not.

Roland
  • 127,288
  • 10
  • 191
  • 288
  • Thanks! This solved my problem. Any ideas on the error message (for future reference...) ? – Dusan Kojic Jan 04 '15 at 17:01
  • Sorry, no idea. If `confint` fails that's often a sign of a bad fit, but from plotting it I can't see a problem. – Roland Jan 04 '15 at 17:25
  • Side note: if you don't use `poly` then be sure to enclose your formula in `I()` to avoid having `*` and `+` interpreted as combining independent factors. e.g. `lm(y ~ I(a*x^2+b*x+c))` – Carl Witthoft Jan 04 '15 at 17:41