如何使用 nlsLM 避免“错误：初始参数估计时的奇异梯度矩阵”？

Question

这是我的第一篇文章，所以我希望我的所有内容都正确无误，并感谢您提供的任何帮助。

我正在尝试将曲线拟合到我收集的一组原始数据。数据是 sigmoid，应该用 Shae-Ackers 模型（生物学中常用）来解释。在查看了此站点上的其他示例后，我最初尝试使用“nls”功能来执行此操作。但是我收到以下错误：

“错误：初始参数估计处的奇异梯度矩阵”。

随后，我尝试使用“nlsLM”功能。正如我从本网站上的其他问题中了解到的那样，这应该使用替代算法来适应我的模型，从而避免错误。但是，我仍然得到同样的错误：

“错误：初始参数估计处的奇异梯度矩阵”。

这是我目前正在尝试工作的代码：

library("minpack.lm") 
library("nls2")

L = c(0, 0.0001, 0.0005, 0.001, 0.005, 
0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50) # X- values, inducer concentration.
Fluo = c(23.10263117,
     21.5111054,
     25.01080225,
     32.63906667,
     82.6671047,
     287.1788694,
     812.2339928,
     1308.71973,
     2822.260637,
     3675.085354,
     4634,
     4399.131349,
     4096.759224) # Y-values, fluorescence measurement.

model <- nlsLM(Fluo ~ Fmax*((k1 + k3C0*((L^n)/(kd^n + L^n)))/(1 + k1 + k2C0*((L^n)/(kd^n + L^n)) + k3C0*((L^n)/(kd^n + L^n)))),
           start = c(k1 = 0.009,k2C0 = 37.5,k3C0 = 3.4,kd = 0.09,n = 2.8,Fmax = 7650), algorithm = "LM", trace = TRUE)

plot(log10(L), log10(Fluo), main = "data")
lines(log10(L), log10(fitted(model)), col = "red", lty = 2)  

我提供的初始估计值应该与我正在寻找的值非常接近，因为这些值之前已在文献中针对我试图拟合的系统进行了报道。

所以我的总体问题是这样的：
1. 有没有办法可以找到可以避免此错误的替代起始估计值？
2. 除了“nls”或“nlsLM”之外，我还可以使用哪些其他函数来拟合我的模型？ （我也看过'nls.lm'，但不确定如何正确实施）。

非常感谢您的帮助，如果需要，我很乐意提供任何其他信息！

Answer 1

问题不在于使用哪个优化器。问题是模型参数无法识别。

设右手边是：

rhs <- function(k1, k2C0, k3C0, kd, n, Fmax, L) {
  r <- L^n / (kd^n + L^n)
  Fmax*((k1 + k3C0 * r)/(1 + k1 + k2C0*r + k3C0*r))
}

然后注意这两个对 rhs 的调用给出了相同的值（除了 0 是退化的）：

st <- c(k1 = 0.009, k2C0 = 37.5, k3C0 = 3.4, 
  kd = 0.09, n = 2.8, Fmax = 7650)

rhs1 <- with(as.list(st), rhs(k1, k2C0, k3C0, kd, n, Fmax, L))

r <- with(as.list(st), L^n/(kd^n+L^n))
rhs2 <- with(as.list(st), rhs(2*k1, (k2C0 * r - k1 - r * k3C0)/r, 2*k3C0, 
 kd, n, Fmax/2, L))

all.equal(rhs1[-1], rhs2[-1])
## [1] TRUE

奇异值

我们可以查看奇异值，我们看到最小的两个奇异值非常小，甚至第三个最小的也比最大的小几个数量级，所以我们将不得不修复 2 或 3 个参数为了应用任何优化算法。

library(nls2)
fm2 <- nls2(Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L), st = st, alg = "brute")
svd(fm2$m$Rmat())[-2]

给予：

$`d`
[1] 1.815654e+04 2.063566e+03 4.100935e+02 4.992959e+01 1.234634e-06
[6] 1.526028e-09

$v
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,] -9.993013e-01  0.0373287823  0.0014810192  1.127961e-03 -2.744777e-09
[2,]  9.184235e-05  0.0058555553 -0.0879560839  3.069373e-03 -9.960997e-01
[3,] -1.388630e-03 -0.0755201221  0.9926417526 -3.431623e-02 -8.820121e-02
[4,]  3.726825e-02  0.9958422212  0.0744652781 -3.692667e-02 -8.316237e-04
[5,]  2.458570e-03  0.0341651105  0.0371291870  9.987232e-01 -2.370720e-09
[6,] -1.845720e-06 -0.0000361664  0.0004803576 -1.661631e-05  2.309941e-03
              [,6]
[1,] -1.187056e-06
[2,]  2.343452e-03
[3,] -2.763882e-04
[4,]  1.622277e-06
[5,] -1.799465e-11
[6,]  9.999972e-01

修正参数

例如，如果我们将参数 2、5 和 6 (k2C0, n, Fmax) 固定在问题中给出的起始值，那么我们会得到以下结果（也显示为最后图中的蓝色拟合）。

fm <- with(as.list(st), nls(Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L), 
  start = st[-c(2, 5:6)]))

给予：

Nonlinear regression model
  model: Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L)
   data: parent.frame()
      k1     k3C0       kd 
 0.04214 49.21873  2.02506 
 residual sum-of-squares: 1737050

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.947e-06

注意

如果您愿意使用不同的公式，那么这个公式只有 4 个参数（与问题中的 6 个参数相比）并提供了合理的拟合。它在最后的图表中显示为红色拟合。

顺便说一句，在问题中，log10(Fluo) 是针对 log10(L) 绘制的，但是最小化 log10(Flou) - log10(rhs(...)) 的平方和与最小化Flou - rhs(...) 的平方和，尽管答案可能很接近。

fm2 <- nls(Fluo ~ L^d / (a + b * L^c), 
 start = c(a = 10e-5, b = 10e-5, c = 1, d = 1))

plot(Fluo ~ L)
lines(fitted(fm) ~ L, col = "blue")
lines(fitted(fm2) ~ L, col = "red")
legend("bottomright", legend = c("fm", "fm2"),
  col = c("blue", "red"), lty = 1, lwd = 2)