Fitting non-linear Langmuir Isotherm in R

Question

I want to fit Isotherm models for the following data in R. The simplest isotherm model is Langmuir model given here model is given in the bottom of the page. My MWE is given below which throw the error. I wonder if there is any R package for Isotherm models.

X <- c(10, 30, 50, 70, 100, 125)
Y <- c(155, 250, 270, 330, 320, 323)
Data <- data.frame(X, Y)
LangIMfm2 <- nls(formula = Y ~ Q*b*X/(1+b*X),  data = Data, start = list(Q = 1, b = 0.5), algorith = "port")

Error in nls(formula = Y ~ Q * b * X/(1 + b * X), data = Data, start = list(Q = 1,  : 
  Convergence failure: singular convergence (7)

Edited

Some nonlinear models can be transform to linear models. My understanding is that there might be one-to-one relationship between the estimates of nonlinear model and its linear model form but their corresponding standard errors are not related to each other. Is this assertion true? Are there any pitfalls in fitting Nonlinear Models by transforming to linearity?

Answer 1

I am not aware of such packages and personally I don't think that you need one as the problem can be solved using a base R.

nls is sensitive to the starting parameters, so you should begin with a good starting guess. You can easily evaluate Q because it corresponds to the asymptotic limit of the isotherm at x-->Inf, so it is reasonable to begin with Q=323 (which is the last value of Y in your sample data set).

Next, you could do plot(Data) and add a line with an isotherm that corresponds to your starting parameters Q and b and tweak b to come up with a reasonable guess.

The plot below shows your data set (points) and a probe isotherm with Q = 323 and b = 0.5, generated by with(Data,lines(X,323*0.5*X/(1+0.5*X),col='red')) (red line). It seemed a reasonable starting guess to me, and I gave it a try with nls:

LangIMfm2 <- nls(formula = Y ~ Q*b*X/(1+b*X),  data = Data, start = list(Q = 300, b = 1), algorith = "port")
# Nonlinear regression model
#   model: Y ~ Q * b * X/(1 + b * X)
#    data: Data
#        Q        b 
# 366.2778   0.0721 
#  residual sum-of-squares: 920.6
# 
# Algorithm "port", convergence message: relative convergence (4)

and plotted predicted line to make sure that nls found the right solution:

lines(Data$X,predict(LangIMfm2),col='green')

Having said that, I would suggest to use a more effective strategy, based on the linearization of the model by rewriting the isotherm equation in reciprocal coordinates:

z <- 1/Data
plot(Y~X,z)
abline(lm(Y~X,z))
M <- lm(Y~X,z)

Q <- 1/coef(M)[1]
# 363.2488 

b <- coef(M)[1]/coef(M)[2]
# 0.0741759 

As you could see, both approaches produce essentially the same result, but the linear model is more robust and doesn't require starting parameters (and, as far as I remember, it is the standard way of the isotherm analysis in the experimental physical chemistry).

Answer 2

You can use the SSmicmen self-starter function (see Ritz and Streibig, 2008, Nonlinear Regression with R) in the nlme package for R, which calculates initial parameters from the fit of the linearized form of the Michaelis-Menten (MM) equation. Fortunately, the MM equation possesses a form that can be adapted for the Langmuir equation, S = Smax*x/(KL + x). I've found the nlshelper and tidyverse packages useful for modeling and exporting the results of the nls command into tables and plots, particularly when modeling sample groups. Here's my code for modeling a single set of sorption data:

library(tidyverse)
library(nlme)
library(nlshelper)

lang.fit <- nls(Y ~ SSmicmen(X,Smax,InvKL), data=Data)
fit.summary <- tidy(lang.fit)
fit.coefs <- coef(lang.fit)

For simplicity, the Langmuir affinity constant is modeled here as 1/KL. Applying this code, I get the same parameter estimates as @Marat given above.

The simple code below allows for wrangling the data in order to create a ggplot object, containing the original points and fitted line (i.e., geom_point would represent the original X and Y data, geom_line would represent the original X plus YHat).

FitY <- tibble(predict(lang.fit))
YHat <- FitY[,1]
Data2 <- cbind(Data, YHat)

If you want to model multiple groups of data (say, based on a "Sample_name" column, then the lang.fit variable would be calculated as below, this time using the nlsList command:

lang.fit <- nlsList(Y ~ SSmicmen(X,Smax,InvKL) | Sample_name, data=Data)

Answer 3

The problem is the starting values. We show two approaches to this as well as an alternative that converges even using the starting values in the question.

1) plinear The right hand side is linear in Q*b so it would be better to absorb b into Q and then we have a parameter that enters linearly so it is easier to solve. Also with the plinear algorithm no starting values are needed for the linear parameter so only the starting value for b need be specified. With plinear the right hand side of the nls formula should be specified as the vector that multiplies the linear parameter. The result of running nls giving fm0 below will be coefficients named b and .lin where Q = .lin / b.

We already have our answer from fm0 but if we want a clean run in terms of b and Q rather than b and .lin we can run the original formula in the question using the starting values implied by the coefficients returned by fm0 as shown.

fm0 <- nls(Y ~ X/(1+b*X), Data, start = list(b = 0.5), alg = "plinear")

st <- with(as.list(coef(fm0)), list(b = b, Q = .lin/b))
fm <- nls(Y ~ Q*b*X/(1+b*X), Data, start = st)
fm

giving

Nonlinear regression model
  model: Y ~ Q * b * X/(1 + b * X)
   data: Data
       b        Q 
  0.0721 366.2778 
 residual sum-of-squares: 920.6

Number of iterations to convergence: 0 
Achieved convergence tolerance: 9.611e-07

We can display the result. The points are the data and the red line is the fitted curve.

plot(Data)
lines(fitted(fm) ~ X, Data, col = "red")

(contineud after plot)

2) mean Alternately, using a starting value of mean(Data$Y) for Q seems to work well.

nls(Y ~ Q*b*X/(1+b*X), Data, start = list(b = 0.5, Q = mean(Data$Y)))

giving:

Nonlinear regression model
  model: Y ~ Q * b * X/(1 + b * X)
   data: Data
       b        Q 
  0.0721 366.2779 
 residual sum-of-squares: 920.6

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.818e-06

The question already had a reasonable starting value for b which we used but if one were needed one could set Y to Q*b so that they cancel and X to mean(Data$X) and solve for b to give b = 1 - 1/mean(Data$X) as a possible starting value. Although not shown using this starting value for b with mean(Data$Y) as the starting value for Q also resulted in convergence.

3) optim If we use optim the algorithm converges even with the initial values used in the question. We form the residual sum of squares and minimize that:

rss <- function(p) {
  Q <- p[1]
  b <- p[2]
  with(Data, sum((Y - b*Q*X/(1+b*X))^2))
}
optim(c(1, 0.5), rss)

giving:

$par
[1] 366.27028219   0.07213613

$value
[1] 920.62

$counts
function gradient 
     249       NA 

$convergence
[1] 0

$message
NULL