Search code examples
rnls

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

Solution

  • 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