Nawa Raj
Nawa Raj

Reputation: 111

Non linear fit in r

My data consist of two columns - time and cumulative number like below:

time <- c(1:14)
cum.num <- c(20, 45, 99, 195, 301, 407, 501, 582, 679, 753, 790, 861, 1011, 1441)

My non linear function is:

B/(B*C*exp(-A*B*time) + 1)

My objective is to model my data using non linear method. I tried the following:

m1 < -nls(cum.num ~ B/((B*C)*exp(-A*B*time) + 1)

I have tried several initial values but got the following error:

Error in nls(cum.vul ~ B/((B * C) * exp(-A * B * time) + 1), 
start = list(B = 60,  : singular gradient

Upvotes: 0

Views: 156

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226881

When this happens you generally have to do a bit of detective work to understand the parameters of your function and get a rough estimate of the values by looking at plots.

time <- c(1:14)
cum.num <- c(20, 45, 99, 195, 301, 407, 501, 582, 679,
           753, 790, 861, 1011, 1441)
  • The first thing to notice is that the top-level structure of the function is hyperbolic (i.e., of the form 1/(1+x)). It will be easier to visualize the data and estimate parameters if we invert the y value, so that we have 1/cum.num ~ C*exp(-A*B*time) + 1/B.
plot(time,1/cum.num,log="y",xlim=c(0,14),ylim=c(0.0005,0.5))

(plotting on a log-y scale, and extending the y-axis limits, based on what I discovered below ...)

enter image description here

  • From the equation above, we know that the asymptote (value for large y) should be 1/B, so we have 1/B ~ 0.001 or B ~ 1000.
  • at time 0, the value should be C + 1/B = C + 0.001. Looking at the graph, we have C ~ 0.5
  • Finally, 1/(A*B) is the characteristic scale of decrease (i.e. the time for an e-fold decrease). It looks like the e-folding time ~ 1 (from t=1 to t=2) so 1/(A*B) ~ 1 so A ~ 0.001

Use these starting values to fit:

m1 <- nls(cum.num ~ B/((B*C)*exp(-A*B*time) + 1),
                   start=list(A=0.001,B=1000,C=0.5))

Seems to work. Plot the predicted values:

tpred <- seq(0,14,length=101)
cpred <- predict(m1,newdata=data.frame(time=tpred))
par(las=1,bty="l")
plot(time,cum.num)
lines(tpred,cpred)

enter image description here

Upvotes: 1

Related Questions