Optimizing Non-Linear Langmuir Parameter Estimation in R

Question

I am interested in using nls to aid in fitting the Langmuir Equation Y =(Qmax*k*X)/(1+(k*X)) similar to what was done in this post Fitting Non-linear Langmuir Isotherm in R. The parameter of the equation I am interested in is Qmax which corresponds to the horizontal asymptote (green line) of the plotted sorption data below. Is there a more robust approach other than nls or a way to improve my use of nls that I could employ to get a Qmax value as close as possible to the visual asymptote (green line) around Qmax=3200?

Lang <- nls(formula = Y ~ (Qmax*k*X)/(1+(k*X)),  data = data, start = list(Qmax = 3600, k = 0.015), algorith = "port") 

Using the following data:

     X        Y
1    3.08   84.735
2    5.13  182.832
3    6.67  251.579
4    9.75  460.077
5   16.30  779.350
6   25.10  996.540
7   40.80 1314.739
8   68.90 1929.422
9  111.00 2407.668
10 171.00 3105.850
11 245.00 3129.240
12 300.00 3235.000

I'm getting a Qmax = 4253.63 (red line) - approximately 1000 units away. Using upper and lower limits only results in a Qmax of what I set the upper limit to and changing initial values doesn't appear to change the outcome. Is this a challenge that can be solved with a different approach to non-linear regression than I've taken in base R or is this a statistical/mathematical problem first and foremost?

Plot of Non-linear Langmuir Isotherm

 summary(Lang)

    Formula: Y ~ (Qmax * k * X)/(1 + (k * X))

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
Qmax 4.254e+03  1.554e+02   27.37 9.80e-11 ***
k    1.209e-02  1.148e-03   10.53 9.87e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 99.14 on 10 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

My attempt at the linearization of the model was less successful:

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

Qmax <- 1/coef(M)[1]
#4319.22

k <- coef(M)[1]/coef(M)[2]
#0.00695

Disclaimer: This is my first post so please bear with me, and I'm relatively new to R. With that being said any technical advice that might help me improve my technique above would be greatly appreciated.

Answer 1

Not sure why you expect Qmax to be that low

I rewrote your dependency in a simplest form, removing multiplication and replacing it with addition (a => 1/k) by dividing both nominator and denominator by k. Result looks perfect to my eye.

library(ggplot2)
library(data.table)

dt <- fread("R/Langmuir.dat", sep = " ")

Lang <- nls(formula = Y ~ (Qmax*X)/(a+X),  data = dt, start = list(Qmax = 3600, a = 100.0), algorithm = "port")
q <- summary(Lang)

Qmax <- q$coefficients[1,1]
a    <- q$coefficients[2,1]

f <- function(x, Qmax, a) {
    (Qmax*x)/(a+x)
}

p <- ggplot(data = dt, aes(x = X, y = Y))
p <- p + geom_point()
p <- p + xlab("T") + ylab("Q") + ggtitle("Langmuir Fit")
p <- p + stat_function(fun = function(x) f(x, Qmax=Qmax, a=a))
print(p)

print(Qmax)
print(a)

Output

4253.631
82.68501

Graph

UPDATE

Basically, too many points at low X, hard to get curve bending for lower Qmax. Designed way to make curve bend is to add weights. For example, if I add weights column after reading data table:

dt[, W := (as.numeric(N)/12.0)^3]

and run nls with weights

Lang <- nls(formula = Y ~ (Qmax*X)/(a+X),  data = dt, start = list(Qmax = 3600, a = 100.0), weights = dt$W, algorithm = "port")

I'll get Qmax and a

[1] 4121.114
[1] 74.89386

with the following graph