R nlsModel 奇异梯度矩阵在初始参数估计误差初始参数估计问题答案

【问题标题】：R nlsModel singular gradient matrix at initial parameter estimates error initial parameter estimation problemR nlsModel 奇异梯度矩阵在初始参数估计误差初始参数估计问题
【发布时间】：2020-11-24 01:25:24
【问题描述】：

我需要针对作物/年份数据进行 nls 拟合（请参阅 Question on nls fit in R - why is this such a strange fit?）。

我从这些数据开始，并使用来自常规线性模型拟合的系数作为初始参数估计值。

x <- c(1979L, 1980L, 1981L, 1982L, 1983L, 1984L, 1985L, 1986L, 1987L, 
1988L, 1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 1997L, 1998L, 
2000L, 2001L, 2002L, 2003L, 2004L, 2012L, 2014L, 2018L)

y <- c(94.4, 100, 120.6, 97.3, 110, 110, 80, 110, 117, 115, 50, 120, 
68.4, 137, 83, 106, 124, 113, 95.8, 115.7, 60, 105, 60, 74, 95.7, 
100)

mod_lm <- lm(y~x)

我现在在使用 minpack.lm 包时遇到此错误，

library(minpack.lm)

mod_nls <- nlsLM(y ~ a+x^b, start=list(a = mod_lm$coefficients[1],b=mod_lm$coefficients[2]))    

Error in nlsModel(formula, mf, start, wts) : 
 singular gradient matrix at initial parameter estimates

或者这个，如果我尝试使用 nls2 包。

library(nls2) 

mod_nls <- nls2(y ~ a+x^b, start=list(a = mod_lm$coefficients[1],b=mod_lm$coefficients[2]))

Error in numericDeriv(form[[3L]], names(ind), env) : 
 Missing value or an infinity produced when evaluating the model

我假设问题是初始参数估计不佳。最初如何获得更好的参数估计？谢谢。

【问题讨论】：

我得到了同样的错误，但在使用x <- 1:26 时没有。在不使用起始参数时，它似乎也能正常工作。
您的解决方案运行良好（只需一次调整）。无论如何，谢谢。

标签： r non-linear-regression nls

【解决方案1】：

部分问题在于您的数据中没有强关系：

cor(x, y)
# [1] -0.219275
cor(x, y, method="spearman")
# [1] -0.2145792
cor(x, y, method="kendall")
# [1] -0.117833
plot(y~x)

从相关性和图中可以看出，x 和 y 之间存在弱负相关。基本上没有强有力的证据表明非线性拟合会优于简单的线性回归。

情节是一条证据线，等级相关性是另一条证据线。如果该关系是一个简单的多项式，则 Spearman 和 Kendall 相关性将大于 Pearson 相关性。正如@Erik Krantz 指出的那样，让函数猜测初始值 (1, 1) 比使用线性回归系数效果更好。

这是您的初始条件描述的线与数据的对比图：

plot(x, coef(mod_lm)[1] + x^coef(mod_lm)[2], ylim=c(0, 1000), type="l")
points(x, y)

看图就知道第一个值太大了，100 左右的起始值会更好：

mod_nls <- nlsLM(y ~ a+x^b, start=list(a = 100, b = coef(mod_lm)[2]))
mod_nls
# Nonlinear regression model
#   model: y ~ a + x^b
#    data: parent.frame()
#       a     b.x 
# 98.1638 -0.1306 
#  residual sum-of-squares: 12142
# 
# Number of iterations to convergence: 13 
# Achieved convergence tolerance: 1.49e-08

【讨论】：

是的，我意识到在这种情况下，简单的线性拟合就可以了，但我无法查看 1000 个这些图来确定这一点，所以我需要一种方法来准确地自动化它。
第一个检查可能是将 Spearman/Kendall 相关性与 Pearson 进行比较。

【解决方案2】：

我们可以使用 nls2 来获得近似拟合，我们看到它几乎是一条水平线。使用下面显示的 b 系数得到大约 x^b = 1，因此第二项是常数，可以吸收到 a 中，我们也可以使用模型 y ~ a 和在这个新的简化模型 a估计为mean(y)。关键是b 的任何足够大的负值都会给出这个结果，因此模型被过度参数化了，我们可以将其减少到仅a 项。

library(nls2)

fm <- nls2(y ~ a + x^b, start = list(a = range(y), b = c(-10, 10)), alg = "brute")
fm
## Nonlinear regression model
##   model: y ~ a + x^b
##    data: parent.frame()
##      a      b 
## 99.714 -4.286 
##  residual sum-of-squares: 12179
##
## Number of iterations to convergence: 64 
## Achieved convergence tolerance: NA

plot(y ~ x)
lines(fitted(fm) ~ x, col = "red")

【讨论】：