do a nonlinear least square fit in r

Question

I have two vectors:

y <- c(0.044924, 0.00564, 0.003848, 0.002385, 0.001448, 0.001138, 
0.001025, 0.000983, 0.00079, 0.000765, 0.000721, 0.00061, 0.000606, 
0.000699, 0.000883, 0.001069, 0.001226, 0.001433, 0.00162, 0.001685, 
0.001604, 0.001674, 0.001706, 0.001683, 0.001505, 0.001497, 0.001416, 
0.001449, 0.001494, 0.001544, 0.00142, 0.001458, 0.001544, 0.001279, 
0.00159, 0.001756, 0.001749, 0.001909, 0.001885, 0.002063, 0.002265, 
0.002137, 0.002391, 0.002619, 0.002733, 0.002957, 0.003244, 0.003407, 
0.003563, 0.003889, 0.004312, 0.004459, 0.004946, 0.005248, 0.005302, 
0.00574, 0.006141, 0.006977, 0.007386, 0.007843, 0.008473, 0.008949, 
0.010164, 0.010625, 0.011279, 0.01191, 0.012762, 0.014539, 0.01477)

x <- 0:68

I am trying to use the non-linear least squares function to fit the data but I keep getting the error:

Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates

My code is:

a=0.00012
b=0.08436
k=0.21108
fit = nls(y ~ (a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), start=list(a=a,b=b,k=k))

The parameters I have entered are parameters that I know are close to the expected values. Does anyone know what am I doing wrong here?

I have tried various initial values for the parameters a, b and k, but I always get some kind of error.

When I compare your starting curve to the points, for example with `plot(x, y); lines((a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), col = "red")`, it doesn't look like it's that similar- in particular the bump near 0 is not represented. Why are you confident that the data can be fit by this curve, and that the parameter values are close? — David Robinson, Aug 05 '15 at 13:23
I'm with @DavidRobinson. It doesn't look like the equation you are using is capable of fitting an upturned parabolic curve. I used your equation with a range of starting values and the only curves with upturns in x less than 10 and x greater than 60 were strange looking, with sharp angles in them. — Jean V. Adams, Aug 05 '15 at 15:53
Thank you very much @DavidRobinson and Jean V. Adams. I am agree both of you. In this case, y=death rate, x=age, and the model use is mixing distribution (frailty models). `(a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1)))` it's gamma distribution with Gompertz be the baseline hazard. How can I find the fit of this? Need your helps. Thank you very much — Desak Ristia Kartika, Aug 05 '15 at 16:34

Dag Hjermann · Answer 1 · 2015-08-06T10:54:39.417

Use optim() instead. You have to make a function which takes a,b and k as input (collected as a vector), and which returns the squared error as a result:

func <-function(pars) {
  a <- pars["a"]
  b <- pars["b"]
  k <- pars["k"]
  fitted <- (a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1)))
  sum((y-fitted)^2)  
  }

Then we run optim() using the initial values:

result <- optim(c(a=0.00012, b=0.08436, k=0.21108), func)

To test the resulting fit:

plot(x, y)
a <- result$par["a"]
b <- result$par["b"]
k <- result$par["k"]
lines((a*b*exp(b*x)*k)/((k*b)+(a*(exp(b*x)-1))), col = "blue")

do a nonlinear least square fit in r

1 Answers1