在 R 中拟合非线性朗缪尔等温线
Fitting non-linear Langmuir Isotherm in R
我想为 R 中的以下数据拟合等温模型。最简单的等温模型是此处给出的 Langmuir 模型 model is given in the bottom of the page。下面给出了我的 MWE，它抛出了错误。不知道有没有等温模型的R包。
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)
已编辑
一些非线性模型可以转换为线性模型。我的理解是非线性模型的估计与其线性模型形式之间可能存在一对一的关系，但它们对应的标准误差彼此无关。这个说法是真的吗？通过转换为线性来拟合非线性模型是否有任何陷阱？
我不知道有这样的包，我个人认为你不需要一个，因为这个问题可以使用基本 R 来解决。
nls 对起始参数很敏感，所以你应该从一个好的起始猜测开始。您可以轻松计算 Q，因为它对应于 x-->Inf 处等温线的渐近极限，因此从 Q=323 开始是合理的（这是 Y 的最后一个值在您的示例数据集中）。
接下来，您可以 plot(Data) 并添加一行，其中包含与您的起始参数 Q 和 b 相对应的等温线，并调整 b 以得出一个合理的猜测。
下图显示了您的数据集（点）和由 with(Data,lines(X,323*0.5*X/(1+0.5*X),col='red')) 生成的 Q = 323 和 b = 0.5 的探针等温线（红线）。对我来说这似乎是一个合理的开始猜测，我试了一下 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)
并绘制预测线以确保 nls 找到正确的解决方案：
lines(Data$X,predict(LangIMfm2),col='green')
话虽如此，我建议使用更有效的策略，基于模型的线性化，通过重写倒坐标中的等温线方程：
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
如您所见，这两种方法产生的结果基本相同，但线性模型更稳健并且不需要起始参数（据我所知，这是等温线分析的标准方法在实验物理化学中）。
您可以使用 R 的 nlme 包中的 SSmicmen 自启动函数（参见 Ritz 和 Streibig，2008，Nonlinear Regression with R），它根据 Michaelis-Menten 线性化形式的拟合计算初始参数(MM) 方程。幸运的是，MM 方程具有适用于 Langmuir 方程的形式，S = Smax*x/(KL + x)。我发现 nlshelper 和 tidyverse 包对于建模和将 nls 命令的结果导出到表格和图中非常有用，尤其是在对样本组建模时。这是我对一组吸附数据建模的代码：
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)
为简单起见，此处将朗缪尔亲和常数建模为 1/KL。应用此代码，我得到与上面给出的@Marat 相同的参数估计值。
下面的简单代码允许整理数据以创建一个 ggplot 对象，其中包含原始点和拟合线（即，geom_point 代表原始 X 和 Y 数据，geom_line 将代表原始 X 加 YHat）。
FitY <- tibble(predict(lang.fit))
YHat <- FitY[,1]
Data2 <- cbind(Data, YHat)
如果您想对多组数据建模（例如，基于“Sample_name”列，那么 lang.fit 变量将按如下方式计算，这次使用 nlsList 命令：
lang.fit <- nlsList(Y ~ SSmicmen(X,Smax,InvKL) | Sample_name, data=Data)
问题出在起始值上。我们展示了两种方法以及一种甚至使用问题中的起始值收敛的替代方法。
1) plinear 右手边在Q*b中是线性的所以最好将b吸收到Q中然后我们有一个线性进入的参数所以它是更容易解决。此外，对于线性算法，线性参数不需要起始值，因此只需指定 b 的起始值。对于 plinear，nls 公式的右侧应指定为与线性参数相乘的向量。 运行ning nls 在下面给出 fm0 的结果将是名为 b 和 .lin 的系数，其中 Q = .lin / b.
我们已经从 fm0 那里得到了答案，但是如果我们想要一个干净的 运行 b 和 Q 而不是 b 和 .lin 我们可以 运行 问题中的原始公式，使用 fm0 返回的系数隐含的起始值，如图所示。
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
给予
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
我们可以显示结果。点是数据，红线是拟合曲线。
plot(Data)
lines(fitted(fm) ~ X, Data, col = "red")
（剧情后续）
2) mean 或者，对 Q 使用 mean(Data$Y) 的起始值似乎效果很好。
nls(Y ~ Q*b*X/(1+b*X), Data, start = list(b = 0.5, Q = mean(Data$Y)))
给予：
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
这个问题已经有了我们使用的 b 的合理起始值，但如果需要的话，可以将 Y 设置为 Q*b 以便它们取消并 X 表示 (Data$X) 并求解 b 以给出 b = 1 - 1/mean(Data$X) 作为可能的起始值。尽管未显示使用 b 的起始值和 mean(Data$Y) 作为 Q 的起始值也导致收敛。
3) optim 如果我们使用 optim 算法即使在问题中使用初始值也会收敛。我们形成残差平方和并将其最小化：
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)
给予：
$par
[1] 366.27028219 0.07213613
$value
[1] 920.62
$counts
function gradient
249 NA
$convergence
[1] 0
$message
NULL
我想为 R 中的以下数据拟合等温模型。最简单的等温模型是此处给出的 Langmuir 模型 model is given in the bottom of the page。下面给出了我的 MWE，它抛出了错误。不知道有没有等温模型的R包。
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)
已编辑
一些非线性模型可以转换为线性模型。我的理解是非线性模型的估计与其线性模型形式之间可能存在一对一的关系，但它们对应的标准误差彼此无关。这个说法是真的吗？通过转换为线性来拟合非线性模型是否有任何陷阱？
我不知道有这样的包，我个人认为你不需要一个，因为这个问题可以使用基本 R 来解决。
nls 对起始参数很敏感，所以你应该从一个好的起始猜测开始。您可以轻松计算 Q，因为它对应于 x-->Inf 处等温线的渐近极限，因此从 Q=323 开始是合理的（这是 Y 的最后一个值在您的示例数据集中）。
接下来，您可以 plot(Data) 并添加一行，其中包含与您的起始参数 Q 和 b 相对应的等温线，并调整 b 以得出一个合理的猜测。
下图显示了您的数据集（点）和由 with(Data,lines(X,323*0.5*X/(1+0.5*X),col='red')) 生成的 Q = 323 和 b = 0.5 的探针等温线（红线）。对我来说这似乎是一个合理的开始猜测，我试了一下 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)
并绘制预测线以确保 nls 找到正确的解决方案：
lines(Data$X,predict(LangIMfm2),col='green')
话虽如此，我建议使用更有效的策略，基于模型的线性化，通过重写倒坐标中的等温线方程：
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
如您所见，这两种方法产生的结果基本相同，但线性模型更稳健并且不需要起始参数（据我所知，这是等温线分析的标准方法在实验物理化学中）。
您可以使用 R 的 nlme 包中的 SSmicmen 自启动函数（参见 Ritz 和 Streibig，2008，Nonlinear Regression with R），它根据 Michaelis-Menten 线性化形式的拟合计算初始参数(MM) 方程。幸运的是，MM 方程具有适用于 Langmuir 方程的形式，S = Smax*x/(KL + x)。我发现 nlshelper 和 tidyverse 包对于建模和将 nls 命令的结果导出到表格和图中非常有用，尤其是在对样本组建模时。这是我对一组吸附数据建模的代码：
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)
为简单起见，此处将朗缪尔亲和常数建模为 1/KL。应用此代码，我得到与上面给出的@Marat 相同的参数估计值。
下面的简单代码允许整理数据以创建一个 ggplot 对象，其中包含原始点和拟合线（即，geom_point 代表原始 X 和 Y 数据，geom_line 将代表原始 X 加 YHat）。
FitY <- tibble(predict(lang.fit))
YHat <- FitY[,1]
Data2 <- cbind(Data, YHat)
如果您想对多组数据建模（例如，基于“Sample_name”列，那么 lang.fit 变量将按如下方式计算，这次使用 nlsList 命令：
lang.fit <- nlsList(Y ~ SSmicmen(X,Smax,InvKL) | Sample_name, data=Data)
问题出在起始值上。我们展示了两种方法以及一种甚至使用问题中的起始值收敛的替代方法。
1) plinear 右手边在Q*b中是线性的所以最好将b吸收到Q中然后我们有一个线性进入的参数所以它是更容易解决。此外，对于线性算法，线性参数不需要起始值，因此只需指定 b 的起始值。对于 plinear，nls 公式的右侧应指定为与线性参数相乘的向量。 运行ning nls 在下面给出 fm0 的结果将是名为 b 和 .lin 的系数，其中 Q = .lin / b.
我们已经从 fm0 那里得到了答案，但是如果我们想要一个干净的 运行 b 和 Q 而不是 b 和 .lin 我们可以 运行 问题中的原始公式，使用 fm0 返回的系数隐含的起始值，如图所示。
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
给予
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
我们可以显示结果。点是数据，红线是拟合曲线。
plot(Data)
lines(fitted(fm) ~ X, Data, col = "red")
（剧情后续）
2) mean 或者，对 Q 使用 mean(Data$Y) 的起始值似乎效果很好。
nls(Y ~ Q*b*X/(1+b*X), Data, start = list(b = 0.5, Q = mean(Data$Y)))
给予：
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
这个问题已经有了我们使用的 b 的合理起始值，但如果需要的话，可以将 Y 设置为 Q*b 以便它们取消并 X 表示 (Data$X) 并求解 b 以给出 b = 1 - 1/mean(Data$X) 作为可能的起始值。尽管未显示使用 b 的起始值和 mean(Data$Y) 作为 Q 的起始值也导致收敛。
3) optim 如果我们使用 optim 算法即使在问题中使用初始值也会收敛。我们形成残差平方和并将其最小化：
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)
给予：
$par
[1] 366.27028219 0.07213613
$value
[1] 920.62
$counts
function gradient
249 NA
$convergence
[1] 0
$message
NULL