r 中 nls 的奇异梯度误差;我怎么知道它是代码还是数据?

singular gradient error with nls in r; how do I know if its the code or the data?

我正在尝试 运行 对此数据进行非线性回归:

Flux<-c(192.09536, 199.47616, 137.63245, 133.60358, -89.28360, -23.17639, -27.14659, 107.25287,  52.72565, NA, 167.43277, 113.59047)
Par<-c(4.166667e-01, 4.347826e-02, 4.583333e-01, 1.845833e+02, 1.122688e+03, 1.059048e+03, 6.384000e+02, 3.326087e+02, 7.094762e+02, 4.180000e+02, 3.953333e+02, 3.998636e+02)
Obs<-c(1,2,3,4,5,6,7,8,9,10,11,12)
curve1<-data.frame(Flux, Par, Obs)
curve1<-do.call("cbind", curve1)

这是我尝试的第一个模型,它已经在其他一些类似的数据集上运行:

model1 <- nls(Flux~b*Par/(c+Par)-a, data = curve1, start=list(a=180,
                                                          b=-200, c=800))

但是对于此数据模型 1 给出:

Error in nls(Flux ~ b * Par/(c + Par) - a, data = curve1, start =
    list(a = 180, : singular gradient

我想这可能是因为我的启动参数不对所以我试着把它变成一个自启动模型(我也尝试了很多不同的启动参数):

model2<-with(curve1, nls(Flux~SSasymp(Par, a, b, c)))

这给出了同样的错误。 但是我认为我在这种情况下错误地使用了 SSasymp,因为它使不正确的曲线适合我能够适合 model1 的数据。 我认为这是因为我混淆了 R 关于 a、b 和 c(?)。我在使用 SSasymp 时读过: b 是 'the horizontal asymptote (a) -the response when x is 0' 而 c 是速率常数。

在我的模型1中的原始方程中,b是水平渐近线,c是速率常数,a是x为0时的响应。

如果我尝试编写一个自启动模型来反映这一点:

model3<-with(curve1, nls(Flux~SSasymp(Par, b, (b-a), c)))

我收到这个错误: 另外: 警告信息: 在 nls(Flux ~ SSalymp(Par, b, (b - a), c)) 中: 没有为某些参数指定起始值。 将“a”初始化为“1”。 考虑指定 'start' 或使用 selfStart 模型

我正在寻求建议 1) model1 不工作是因为我的 code/incorrect 起始参数错误 还是因为模型不适合数据 ?

如果是后者,有没有办法强制R尽最大努力去拟合一个非线性模型呢?从生态学上讲,这应该是一条饱和曲线。

2) I/How 我可以将方程拟合到自启动模型中吗?我是否从根本上误解了如何使用 SSASymp?

非常感谢任何帮助。抱歉,如果我没有很好地解释或格式化它,这是我的第一个 post,我不是经验丰富的 R 用户或统计学家!

是这样的吗?

model1<-nls(Flux~b*Par/(c+Par)-a, data = curve1, start=list(a=180, b=-200, c=-2000))
plot(Flux~Par,curve1)
curve(predict(model1,newdata=data.frame(Par=x)),add=TRUE)
summary(model1)
# Formula: Flux ~ b * Par/(c + Par) - a
#
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# a  -179.17      22.86  -7.837 5.06e-05 ***
# b  1009.36    2556.44   0.395    0.703    
# c -5651.20   11542.41  -0.490    0.638    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 42.43 on 8 degrees of freedom
# ...

您的数据有些病态。形式的函数

y = b * x / (c+x)

b < 0c > 0时是上凹的;当 b > 0c < 0 时,它们向下凹,提供 |c| > max(x)(否则会有垂直渐近线,如评论之一所示)。因为您的数据是 "nearly" 线性的,并且有很大的分散性,所以最佳拟合(例如,最小化残差平方和的参数集 a、b 和 c)是下凹的(c < 0).因此,如果您从估计 c < -max(x) 开始,您就会收敛。

现在,我从你的问题中了解到 c 有一些物理意义,要求它 > 0。这里的问题是你的模型被过度指定(参数太多)。在饱和曲线中,速率常数由曲率确定。但在你的情况下,没有曲率(或者,如果有的话,曲率是负的),所以你无法确定速率常数。在数学上,对于 x << c

b * x / (c + x) ~ (b/c) * x

在您的情况下,斜率约为 -0.25,因此 b/c ~ -0.25。但是 bc 有无数个值可以产生这个比率。因此,虽然您对比率 b/c 了解很多,但您对 bc 分别一无所知。这就是为什么这些参数的标准误差在上面的拟合中如此之大(并且 p 值如此之高)。

底线是,在这种特定情况下,您没有足够的数据来分别准确地确定 a、b 和 c。

[两个小问题]

  1. 您的数据中 NA 的存在与此无关 - nls(...) 默认删除包含 NA 的行。
  2. 您不需要以下行:curve1<-do.call("cbind", curve1)