nls - 收敛错误

Question

对于此数据集:

dat = structure(list(x = c(5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 15L, 
15L, 15L, 15L, 17L, 17L, 17L, 17L, 20L, 20L, 20L, 20L, 20L, 20L, 
20L, 20L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 24L, 25L, 25L, 25L, 
25L, 27L, 27L, 27L, 27L, 30L, 30L, 30L, 30L, 35L, 35L, 35L, 35L), 
y = c(2.2, 2.2, 1.95, 1.9, 4.1, 3.95, 3.75, 3.4, 5.15, 4.6, 
4.75, 5.15, 3.7, 4.1, 3.9, 3.5, 7, 6.7, 6.7, 6.95, 4.95, 6, 6.45, 
6.4, 7, 4.45, 6.15, 6.4, 7, 6.6, 6.7, 7, 4.5, 4.7, 5.75, 4.35, 
5.4, 5.15, 5.7, 5.7, 0, 0, 0.5, 0, 0, 0, 0, 0)), .Names = c("x", "y"), 
row.names = c(6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
15L, 16L, 17L, 34L, 35L, 36L, 37L, 18L, 19L, 20L, 21L, 38L, 39L, 
40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 22L, 23L, 24L, 
25L, 50L, 51L, 52L, 53L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L), 
class = "data.frame")

其中"x"是温度,"y"是生物过程的响应变量

我正在尝试适应这个功能

beta.reg<-function(x, Yopt,Tmin,Topt,Tmax, b1) {
Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
}

mod <- nls(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat,
       start=c(Yopt=6, Tmin=0.1, Topt=24, Tmax=30, b1=1),
       control=nls.control(maxiter=800))

但是,我有这个消息错误:

错误en numericDeriv(form [[3L]],names(ind),env):评估模型时产生的缺失值或无穷大

我已尝试使用其他类似数据集的相同功能并正确拟合...

 rnorm<-(10)
 y <- c(20,60,70,49,10)
 rnorm<-(10)
 y <- c(20,60,70,49,10)
 dat<-data.frame(x = rep(c(15,20,25,30,35), times=5),
              rep = as.factor(rep(1:5, each=5)),
              y = c(y+rnorm(5), y+rnorm(5),y+rnorm(5),y+rnorm(5),y+rnorm(5)))

有人可以帮我吗？

会话信息:

R version 3.1.1 (2014-07-10)
Platform: x86_64-pc-linux-gnu (64-bit)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] nlme_3.1-118        latticeExtra_0.6-26 RColorBrewer_1.0-5  lattice_0.20-29    

loaded via a namespace (and not attached):
[1] grid_3.1.1  tools_3.1.1

Answer 1

这里有很多问题,我怀疑在SO帖子中可以充分涵盖它,但这应该让你开始.

首先,它看起来像你想要的Tmax < max(dat$x),例如,<35.这会导致一个问题,因为那时Tmax - x < 0某些值x以及当你试图将一个负数加到一个幂时(在公式的第二项中),你就得到了NA.这是错误消息的原因.

其次,非线性模型的收敛取决于模型公式以及数据,因此该过程与一组数据而非另一组数据收敛的事实完全无关紧要.

第三,非线性建模迭代地最小化作为参数的函数的残差平方和.如果RSS表面具有 局部最小值,并且您start接近于1,则算法将找到它.但只有全局最小值才是真正的解决方案.你的问题有很多很多的本地最小值.

第四,nls(...)默认使用Gauss Newton方法.高斯牛顿是与移位参数(其被添加到或从预测器中减去,所以参数熟知的不稳定Tmin和Tmax你的情况).幸运的是,该minpak.lm软件包实现了Levenberg Marquardt方法,该方法在这些条件下更加稳定.该nlsLM(...)包中的函数使用nls(...)与返回和类型对象相同的调用序列nls,因此该类对象的所有方法也可以工作.用那个.

第五,非线性回归(实际上是所有最小二乘回归)的基本假设是残差是正态分布的.因此,您必须使用QQ图验证任何解决方案.

第六,你的模型具有一组不正常的特征.当Tmin -> -Inf模型中的第一项接近时1.事实证明,这产生的RSS低于任何其他Tmin小于的值min(dat$x),因此算法都倾向于驱动Tmin大的负值.你可以很容易地看到这一点:

library(minpack.lm)
mod <- nlsLM(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat,
             start=c(Yopt=6,Tmin=0,Topt=24,Tmax=50, b1=1),
             control=nls.lm.control(maxiter=1024,maxfev=1024))
coef(summary(mod))
#         Estimate   Std. Error     t value     Pr(>|t|)
# Yopt    6.347019    0.2919686 21.73870235 8.055342e-25
# Tmin -155.530098 2204.0011003 -0.07056716 9.440694e-01
# Topt   21.157545    0.6702713 31.56564484 2.240134e-31
# Tmax   35.000000   11.4838614  3.04775537 3.933164e-03
# b1      3.321326    9.1844548  0.36162468 7.194035e-01
sum(residuals(mod)^2)
# [1] 50.24696

par(mfrow=c(1,2))
plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))

这看起来非常不错,但事实并非如此:QQ图表显示残差不是很正常.这两个事实Tmin,并b1都很差估计,和值Tmin不是物理意义,与数据,不配合的问题.

第七,事实证明上面的拟合实际上是局部最小值.我们可以做一个网格搜索看到Tmin,Tmax和b1(离开了Yopt,并Topt节省时间,因为这些参数都很好,无论出发点的估计).

init <- c(Yopt=6, Topt=24)
grid <- expand.grid(Tmin= seq(0,4,len=100),
                    Tmax= seq(35,100,len=10),
                    b1  = seq(1,10,len=10))
mod.lst <- apply(grid,1,function(gr){
  nlsLM(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat,
        start=c(init,gr),control=nls.control(maxiter=800)) })
rss <- sapply(mod.lst,function(m)sum(residuals(m)^2))
mod <- mod.lst[[which.min(rss)]]   # fit with lowest RSS
coef(summary(mod))
#        Estimate   Std. Error      t value     Pr(>|t|)
# Yopt   6.389238    0.2534551 25.208557840 2.177168e-27
# Topt  22.636505    0.5605621 40.381798589 7.918438e-36
# Tmin  35.000002  104.6221159  0.334537316 7.396005e-01
# Tmax  36.234602  133.4987344  0.271422809 7.873647e-01
# b1   -41.512912 7552.0298633 -0.005496921 9.956395e-01
sum(residuals(mod)^2)
# [1] 34.24019

plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))

在数学上,这是一个明显优越的拟合:RSS较低,残差更接近正态分布.同样,参数未被很好地估计并且在物理上没有意义的事实是数据(可能是模型公式)的问题,而不是拟合过程.

所有上述内容都表明您的模型存在问题.从数学上讲,它的一个问题是函数在x外部是未定义的(Tmin,Tmax).由于您有数据输出x=35,拟合算法永远不会产生Tmax < 35(如果它收敛).处理此问题的方法会稍微改变您的模型函数,以便在该范围之外剪切为0.(我不知道这是否合法,基于你的问题的物理学,虽然......).

beta.reg<-function(x, Yopt,Tmin,Topt,Tmax, b1) {
  ifelse(x>Tmax,0,
    ifelse(x<Tmin,0,
      Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
  ))
}

使用此函数运行上面的代码会产生:

coef(summary(mod))
#         Estimate   Std. Error     t value     Pr(>|t|)
# Yopt   6.1470413   0.21976766   27.970636 3.202940e-29
# Tmin -52.8172658 184.16899439   -0.286787 7.756528e-01
# Topt  23.0777898   0.63750721   36.200045 7.638121e-34
# Tmax  30.0039413   0.02529877 1185.984187 1.038918e-98
# b1     0.5966129   0.32439982    1.839128 7.280793e-02

sum(residuals(mod)^2)
# [1] 28.10144

par(mfrow=c(1,2))
plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))
qqline(residuals(mod))

实际上,网格搜索产生与起点无关的完全相同的结果.请注意,RSS低于早期模型的任何结果,并且b1估计得更好(与早期模型函数的估计值非常不同).残差仍然不正常,但在这种情况下,我想检查数据是否有异常值.