Warning: lmdif: info = 0. Improper input parameters for nls.lm() function
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.
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))
- Converting the row of a data.table to a vector
- ggplot multi-line log-log labels
- dplyr filter vector by position in pipe %>%
- Removing white spaces in a histogram?
- How to assign and match land cover type from Corine to a dataframe with a set of lon lat coordinates?
- Error in list2env(list, envir = .GlobalEnv) : names(x) must be a character vector of the same length as x
- How to traverse hclust internal node in R
- How are the abbreviations from the par function like mgp, mar, mai, deciphered?
- How to map over function with multiple arguments (and changing data)?
- Join/fill in data based on both row and column values
- How to plot a Log Pearson III distribution using `ggplot2` in R?
- Variable names and Easystats reports
- R plotly(): Adding regression line to a correlation scatter plot
- How to Use ifelse Statements in R
- Creating subgraphs from neighbors with max weight in R using igraph
- merge based on id and dates with missing data
- Automatically adjust the size of the text inside each bar on a bar plot using ggplot2 in R
- gghighlight and facet_wrap one always highlighted
- Mask and Crop functions in Terra R
- How to use a variable in dplyr::filter?
- Changing R Temporary Files
- How do you add a general label to facets in ggplot2?
- change colour of controls in R DT datatable
- Click Button Within Datatable to Display Popup Twice in a row
- How to display values in filter with ascending order
- Select Data From Column in Data Table
- Warning: Error in dataTableOutput: unused argument (height = "auto")
- How to have cells with rounded corners in {gt}?
- Two colours in the same box
- Change color of datatable filter dropdown box