Nonlinear model fitting issue: NaNs produced in minpack.lm::nlsLM

Machine Learning and Modeling
1

I am trying to fit a nonlinear model to empirical data using the minpack.lm::nlsLM function. The initial conditions for the model work, but I am getting this warning message:
In log((dayOfYear_num - tau_2_A)/delta_2_L) : NaNs produced.

I am unable to resolve this problem. I tried to impose conditions (for example, parameters > 0) like this, but it doesn't work:

 lower = c(tau_1_A = 0.0001, tau_2_A = 0.0001, h_1_L = 0.0001, h_2_L = 0.0001, delta_1_L = 0.0001, delta_2_L = 0.0001,
                                    alpha_1_L = 0.0001, alpha_2_L = 0.0001)

Here is the data:

z <- structure(list(dayOfYear = structure(1:33, levels = c("128", 
                                                           "136", "137", "138", "142", "143", "144", "157", "162", "163", 
                                                           "164", "166", "167", "173", "184", "185", "186", "187", "189", 
                                                           "197", "198", "199", "200", "201", "232", "233", "235", "236", 
                                                           "253", "255", "256", "299", "308"), class = "factor"), density = c(0, 
                                                                                                                              0, 0, 0, 0.0153846153846154, 0.0974358974358974, 0.0745192307692308, 
                                                                                                                              0.130769230769231, 0.0307692307692308, 0.0259615384615385, 0.158974358974359, 
                                                                                                                              0.207692307692308, 0.0375, 0.0510989010989011, 0.00769230769230769, 
                                                                                                                              0, 0, 0.00384615384615385, 0, 0.00192307692307692, 0.0201923076923077, 
                                                                                                                              0.00341880341880342, 0.00538461538461539, 0, 1.11923076923077, 
                                                                                                                              0.39957264957265, 1.32735042735043, 0.564102564102564, 0.562820512820513, 
                                                                                                                              0.486410256410256, 0.343956043956044, 0, 0), dayOfYear_num = c(128, 
                                                                                                                                                                                             136, 137, 138, 142, 143, 144, 157, 162, 163, 164, 166, 167, 173, 
                                                                                                                                                                                             184, 185, 186, 187, 189, 197, 198, 199, 200, 201, 232, 233, 235, 
                                                                                                                                                                                             236, 253, 255, 256, 299, 308)), row.names = 54:86, class = "data.frame")

Here is the code:

mod <- minpack.lm::nlsLM(density ~ ifelse(dayOfYear_num < tau_1_A, 0, 
                                          ifelse(dayOfYear_num < tau_2_A, 
                                                 h_1_L*exp(-1/2*((((log((dayOfYear_num - tau_1_A)/delta_1_L))^2)/alpha_1_L)^2)),
                                                 h_2_L*exp(-1/2*((((log((dayOfYear_num - tau_2_A)/delta_2_L))^2)/alpha_2_L)^2)))),
                         data=z,
                         start=list(tau_1_A = 100,
                                    tau_2_A = 150,
                                    h_1_L = 0.5,
                                    h_2_L = 1,
                                    delta_1_L = 50,
                                    delta_2_L = 50,
                                    alpha_1_L = 0.1,
                                    alpha_2_L = 0.1), algorithm="port")
2

Hi @peresophie11 ,

this is a problem with the way ifelse works. Look at this example from the documentation of ifelse itself :

x <- c(6:-4)
#...
## Note: the following also gives the warning !
ifelse(x >= 0, sqrt(x), NA)

you would think that it should be fine - but it still throws a warning message.

So refactoring this part:

ifelse(dayOfYear_num < tau_1_A, 0, 
                                          ifelse(dayOfYear_num < tau_2_A, 
                                                 h_1_L*exp(-1/2*((((log((dayOfYear_num - tau_1_A)/delta_1_L))^2)/alpha_1_L)^2)),
                                                 h_2_L*exp(-1/2*((((log((dayOfYear_num - tau_2_A)/delta_2_L))^2)/alpha_2_L)^2)))),
                         data=z,

to an if - else if or switch statement will get rid of the warnings.

This has nothing to do with minpack.

3

Hi @vedoa,

Thank you very much for your response. I had tried an if - else if, but I was getting this error message: Error in if (dayOfYear_num < tau_1_A) { : the condition has length > 1.

Here is the modified code with an if - else if:

f <- function(dayOfYear_num, tau_1_A, tau_2_A, h_1_L, h_2_L, delta_1_L, delta_2_L, alpha_2_L){
  if(dayOfYear_num < tau_1_A){
    result <- 0
  } else if(dayOfYear_num < tau_2_A){
    result <- h_1_L*exp(-1/2*((((log((dayOfYear_num - tau_1_A)/delta_1_L))^2)/alpha_1_L)^2))
  } else {
    result <- h_2_L*exp(-1/2*((((log((dayOfYear_num - tau_2_A)/delta_2_L))^2)/alpha_2_L)^2))
  }
  
  return(result)
}

mod <- minpack.lm::nlsLM(density ~ f(dayOfYear_num, tau_1_A, tau_2_A, h_1_L, h_2_L, delta_1_L, delta_2_L, alpha_2_L),
                         data=z,
                         start=list(tau_1_A = 100,
                                    tau_2_A = 150,
                                    h_1_L = 0.5,
                                    h_2_L = 1,
                                    delta_1_L = 50,
                                    delta_2_L = 50,
                                    alpha_1_L = 0.1,
                                    alpha_2_L = 0.1), algorithm="port")

Any advice is greatly appreciated.

4

@peresophie11 a few remarks:

algorithm="port" does nothing since according to the documentation of nlsLM algorithm only method "LM" (Levenberg-Marquardt) is implemented. So you should skip that.

Running your initial code i have no idea what R does with a formula having a ifelse statement

Formula: density ~ ifelse(dayOfYear_num < tau_1_A, 0, ifelse(dayOfYear_num < 
    tau_2_A, h_1_L * exp(-1/2 * ((((log((dayOfYear_num - tau_1_A)/delta_1_L))^2)/alpha_1_L)^2)), 
    h_2_L * exp(-1/2 * ((((log((dayOfYear_num - tau_2_A)/delta_2_L))^2)/alpha_2_L)^2))))

The warnings definitely come from there.

So to give you at least something that i understand i rewrote the function as the following:

f <- function(dayOfYear_num, l){
  vapply(dayOfYear_num, function(x){
    if(x < l$tau_1_A){
      return(0)
    } else if(x < l$tau_2_A){
      return(l$h_1_L*exp(-1/2*((((log((x - l$tau_1_A)/l$delta_1_L))^2)/l$alpha_1_L)^2)))
    } else {
      return(l$h_2_L*exp(-1/2*((((log((x - l$tau_2_A)/l$delta_2_L))^2)/l$alpha_2_L)^2)))
    }
  }, numeric(1))

}

fn <- function(l){
 z$density - f(z$dayOfYear_num, l)
}

vapply is just a nice loop ensuring that on every iteration a numeric value will be returned.

Now calling:

mod_new <- minpack.lm::nls.lm(par = list(tau_1_A = 100,
                                    tau_2_A = 150,
                                    h_1_L = 0.5,
                                    h_2_L = 1,
                                    delta_1_L = 50,
                                    delta_2_L = 50,
                                    alpha_1_L = 0.1,
                                    alpha_2_L = 0.1),
                          fn = fn)

Should in principle do what you wanted:

  • based on the entry of your dayOfYear_num vector decide which function to call (for every entry separately)
  • run the Levenberg-Marquardt algorithm
  • has no side effect calls (or whatever is happening inside the formula in your original call)

The results are the same for your initial mod and the mod_new when running summary(mod) and summary(mod_new).

Btw. this will still throw 2 Warnings (not 50 :D) which i think comes from non exhaustive checking in the if else if statement.

5

Thank you so much, @vedoa. Your answers were invaluable to me. Thanks again for your help !