Warning: lmdif: info = 0. Improper input parameters for nls.lm() function

Question

I have the following solutions.

$$\frac{{dx}}{{dt}} = k_+(2d11 + d12 - 2x - z) - k_{-}x$$

$$\frac{{dy}}{{dt}} = k_+(2d22 + d12 - 2y - z) - k_{-}y$$

$$\frac{{dz}}{{dt}} = 2k_+(2d11 + d12 - 2x - z)(2d22 + d12 - 2y - z) - k_{-}z$$

I want to write a program in R to solve the model equations and then fit the data in the table from an article to the solution of equations. I want to find the best fit k+ and k− according to a Levenberg–Marquardt least-squares error algorithm.

The following code is my attempt to do this.

library(deSolve)
library(minpack.lm)

model = function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dxdt = k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x
    dydt = k_plus*(2*d22 + d12 - 2*y - z)-k_minus * y
    dzdt = 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z
    return(list(c(dxdt, dydt, dzdt)))
  })
}

# Defining an objective function
objective = function(par){
  k_plus = par[1]
  k_minus = par[2]
  x0 = 1.71
  y0 = 1.75
  z0 = 0.62
  d11 = 1.71
  d12 = 1.75
  d22 = 0.62
  
# Vector of times
times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)

initial_state = c(x = x0, y = y0, z = z0)
  
  # Solving the ODE
  out = lsode(y = initial_state, 
              times = times_vec,
              func = model,
              parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
              maxsteps = 100000)
  
  error = sum((out[,"x"] - xdata)^2 + (out[,"y"] - ydata)^2 + (out[,"z"] - zdata)^2)
  return(error)
}



# The table from the article
xdata = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05)
ydata = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08)
zdata = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95)


# 
fit = nls.lm(par = c(k_plus = 200,
                     k_minus = 2),
             fn = objective)

When I run the code everything works fine except when using the nls.lm function. I get the following error

Warning: lmdif: info = 0. Improper input parameters.

I have tried different parameters. I looked up the function to make sure my parameters have the correct format but I cant seem to find what the problem is.

Answer 1

To debug the problem, I suggest to compare model and data before attempting to fit the model. We can see clearly, that the model results with supplied numerical values are far away from the data.

library(deSolve)
library(minpack.lm)

model = function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dxdt = k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x
    dydt = k_plus*(2*d22 + d12 - 2*y - z) - k_minus*y
    dzdt = 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z
    return(list(c(dxdt, dydt, dzdt)))
  })
}

times_vec = c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)
k_plus = 200; k_minus = 2
x0 = 1.71; y0 = 1.75; z0 = 0.62
d11 = 1.71; d12 = 1.75; d22 = 0.62
initial_state = c(x = x0, y = y0, z = z0)

out = lsode(y = initial_state, 
            times = times_vec,
            func = model,
            parms = c(k_plus = k_plus, k_minus = k_minus, d11 = d11, d12 = d12, d22 = d22),
            maxsteps = 100000)


# The table from the article
data = data.frame(
  time = times_vec,
  x = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05),
  y = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08),
  z = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95)
)

plot(out, obs=data, mfrow=c(1, 3))

Therefore, I would recommend to check:

• the model equations and
• the numerical values, especially as the initial states are obviously the same as the d..-parameters.

Furthermore, I would recommend to use the lsoda solver instead of lsode. One may also consider another optimization package, e.g. FME, see https://doi.org/10.18637/jss.v033.i03 This package is a wrapper for several optimization algorithms (including nls.lm), but has more support for fitting ODE models.

Solution

The model did not fit well, because the start parameter for k_plus was to extreme. After reducing it to, say 2.0, it was technically possible to fit the model, but the curves were too far away from the data.

Then, after including all parameters, it was possible to get a reasonable fit.

library(deSolve)
library(FME)

model <- function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dxdt = k_plus*(2*d11 + d12 - 2*x - z) - k_minus*x
    dydt = k_plus*(2*d22 + d12 - 2*y - z) - k_minus*y
    dzdt = 2*k_plus*(2*d11 + d12 - 2*x - z)*(2*d22 + d12 - 2*y - z) - k_minus*z
    return(list(c(dxdt, dydt, dzdt)))
  })
}

times_vec <- c(1/30, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 16, 18, 20)
initial_state <- c(x = 1.71, y = 1.75, z = 0.62)

# note different start value for k_plus
pstart <- c(k_plus = 1, k_minus = 1, d11 = 1.71, d12 = 1.75, d22 = 0.62)

# The table from the article
data <- data.frame(
  time = times_vec,
  x = c(1.71, 1.56, 1.40, 1.28, 1.14, 1.17, 1.14, 1.06, 1.10, 1.10, 1.07, 1.09, 1.06, 1.03, 1.06, 1.06, 1.05),
  y = c(1.75, 1.59, 1.43, 1.31, 1.17, 1.20, 1.17, 1.09, 1.13, 1.13, 1.10, 1.13, 1.09, 1.06, 1.09, 1.09, 1.08),
  z = c(0.62, 0.93, 1.26, 1.50, 1.78, 1.71, 1.78, 1.93, 1.86, 1.86, 1.91, 1.87, 1.93, 1.99, 1.94, 1.94, 1.95)
)

## test without fit
out <- ode(y = initial_state, times = times_vec, func = model, parms = pstart)
plot(out, obs=data, mfrow=c(1, 3))

objective <- function(par) {
  # cat(par, ", ") # uncomment for debugging
  out <- ode(y = initial_state,  times = times_vec, func = model, parms = par)
  cost <- modCost(out, data)
  # cat(cost$model, "\n") # uncomment for debugging
 cost
}

# method = "Marq" is essentially from nls.lm
# the model produces some warnings, but finally converges
# fit <- modFit(f = objective, p = pstart, method="Marq") 

# method = "Nelder" runs more smoothly in this case
fit <- modFit(f = objective, p = pstart, method="Nelder") 

# another way would be to use box constraints and then, 
# for example, L-BFGS-B or Port
#fit <- modFit(f = objective, p = pstart, 
#              lower = c(k_plus = .1, k_minus = .1, d11 = .1, d12 = .1, d22 = .1), 
#              upper = c(k_plus = 10, k_minus = 10, d11 = 10, d12 = 10, d22 = 10),
#              method="Port") 

# show detailed results
summary(fit)

# simulate model with fitted parameters
out2 <- ode(y = initial_state, times = times_vec, func = model, parms = fit$par)
plot(out2, obs=data, mfrow=c(1,3))