Search code examples
rselfnls

Cannot use self-starting models when manually defining maxiter for nls()?


Data:

structure(list(ID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 
24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 
37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 
50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 59L, 60L, 61L, 
62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L), Stage = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 3L, 3L, 5L, 5L, 5L, 1L, 1L, 6L, 6L, 
4L, 4L, 2L, 2L, 7L, 7L), .Label = c("milpa", "robir", "jurup che", 
"pak che kor", "mehen che", "nu kux che", "tam che"), class = "factor"), 
    Time.Since.Burn = c(4, 2, 0.21, 2, 0.42, 4, 0.33, 0.33, 3, 
    6, 2.5, 5, 4, 5, 1.5, 6, 4, 6, 3, 6.5, 6.5, 6, 4, 2.5, 12, 
    10, 8, 18, 5, 10, 8, 16, 28, 22, 22, 21, 20, 18, 30, 27, 
    30, 36, 36, 40, 32, 28, 50, 32, 60, 60, 60, 60, 60, 60, 60, 
    60, 6, 6, 24, 26, 22, 2, 1, 50, 45, 10, 10, 4, 4, 60, 60), 
    meandec = c(0.3625, 0.3025, 0.275, 0.1075, 0.26, 0.395, 0.265, 
    0.4075, 0.9, 0.9275, 0.7075, 0.9625, 0.7725, 0.9325, 0.9875, 
    0.81, 0.575, 0.3075, 0.4675, 0.6975, 0.33, 0.8725, 0.46, 
    0.19, 0.495, 0.3825, 0.58, 0.2275, 0.45, 0.3925, 0.605, 0.515, 
    0.425, 0.34, 0.2475, 0.1375, 0.4225, 0.505, 0.36, 0.4325, 
    0.26, 0.1575, 0.125, 0.3125, 0.1725, 0.3175, 0.43, 0.3475, 
    0.2025, 0.395, 0.12, 0.1625, 0.3175, 0.1975, 0.1525, 0.2775, 
    0.4975, 0.725, 0.04, 0.326666666666667, 0.1425, 0.445, 0.4725, 
    0.3775, 0.27, 0.2225, 0.23, 0.3275, 0.9725, 0.215, 0.2325
    )), row.names = c(NA, -71L), class = c("grouped_df", "tbl_df", 
"tbl", "data.frame"), vars = c("ID", "Stage"), drop = TRUE)

Problem:

I'm trying to run an exponential decay model on these data. I've done it with similar data, but when I try to do it on this particular dataset, it says that the number of max iterations has been exceeded without convergence.

nonlinmod6<-nls(meandec~SSasymp(Time.Since.Burn, Asym,R0,lrc),data=averaged_perherb)

Error in nls(y ~ cbind(1 - exp(-exp(lrc) * x), exp(-exp(lrc) * x)), data = xy,  :   number of iterations exceeded maximum of 50

So, I tried to manually increase the maximum number of iterations using the code below:

nonlinmod6<-nls(meandec~SSasymp(Time.Since.Burn, Asym,R0,lrc),data=averaged_perherb,nls.control(maxiter=500))

but it then gives me an error saying that :

Error in nls(meandec ~ SSasymp(Time.Since.Burn, Asym, R0, lrc), data = 
averaged_perherb,: parameters without starting value in 'data': Asym, R0, lrc

which I don't think should be the case given that I'm using a self-starting function to identify the starting parameters. Is there any way to resolve this?


Solution

  • The problem is that the SSaymp intialization routine itself uses nls and it is that hidden invocation of nls that is the problem.

    You are going to have to hack the intialization routine. Make a new copy of SSasymp called SSasymp2, grab its initialization routine and call it SSasymp2Init, say. Then use trace to insert into the initialization a new version of nls having the required control argument. To do that we use the partial function in the pryr package. Replace the initialization routine with the hacked one and then run nls.

    library(pryr)
    
    SSasymp2 <- SSasymp
    SSasymp2Init <- attr(SSasymp2, "initial")
    trace(SSasymp2Init, 
      quote(nls <- partial(stats::nls, control = nls.control(maxiter = 500))))
    attr(SSasymp2, "initial") <- SSasymp2Init
    
    nls(meandec ~ SSasymp2(Time.Since.Burn, Asym, R0, lrc), data = averaged_perherb)
    

    giving:

    Tracing (attr(object, "initial"))(mCall = mCall, data = data, LHS = LHS) on entry 
    Nonlinear regression model
      model: meandec ~ SSasymp2(Time.Since.Burn, Asym, R0, lrc)
       data: averaged_perherb
       Asym      R0     lrc 
     0.1641  0.5695 -3.4237 
     residual sum-of-squares: 2.977
    
    Number of iterations to convergence: 15 
    Achieved convergence tolerance: 5.875e-06