NLS - 误差收敛问题

Question

NLS - 误差收敛问题

3

对于这个数据集：

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))

但是，我遇到了以下错误信息：

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

我已经尝试了另一个类似的数据集，同样的函数可以正确地适配...

 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

- Juanchi

这是在R语言中吗？如果是的话，你应该添加[R]标签。 - John Saunders

2个回答

1

为@jlhoward的解决方案添加另一个可能的解决方案...

我发现这个nls2包：

library("nls2")

从原始数据集中排除x〜17,35：

newdat <- subset(dat, x!=17 & x!=35 )

将该函数应用于精简数据集：

beta.reg<-with(newdat,  
           y ~ Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / Tmax-Topt))^b1
           )

创建一组启动程序：

st1 <- expand.grid(Yopt = seq(4, 8, len = 4),
                   Tmin = seq(0, 4, len = 4), 
                   Topt = seq(15, 25, len = 4),
                   Tmax= seq(28, 38, len = 4),
                   b1 = seq(0, 4, len = 4))

适配模型：

mod <- nls2(beta.reg, start = st1, algorithm = "brute-force")

提取系数：

round(coef(summary(mod)),3)

#     Estimate Std. Error t value Pr(>|t|)
# Yopt    6.667      0.394  16.925    0.000
# Tmin    0.000     12.023   0.000    1.000
# Topt   21.667      0.746  29.032    0.000
# Tmax   31.333      1.924  16.289    0.000
# b1      1.333      1.010   1.320    0.197

诊断：

sum(residuals(mod)^2)

# [1] 50.18246

最后，调整后的函数和QQ正态图：

par(mfrow=c(1,2))
with(newdat,plot(y~x,xlim=c(0,35))) 
points(fitted(mod)~I(newdat$x), pch=19)
with(as.list(coef(mod)),
 curve(
  Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1,
   add=TRUE, col="red"))

qqnorm(residuals(mod))
qqline(residuals(mod))

- Juanchi

记录一下，nls2(...)（您正在使用它）不会最小化RSS，它会在1024个网格点中计算每个点的RSS并报告具有最低RSS的点。这就是为什么您得到Tmin=0的原因。较低的Tmin值将给出较低的RSS，但这是您网格中的最低值。 - jlhoward

是的。通过这种方式，我试图将“Tmin”的估计值限制在一些生物感觉值范围内，牺牲了残差平方和（RSS）。这与您上一个模型的限制是否相同？beta.reg <- function(x, Yopt, Tmin, Topt, Tmax, b1) { 如果（x> Tmax，0，如果（x <Tmin，0，

Yopt * ((x-Tmin) /（Topt-Tmin）) ^（b1 *（Topt-Tmin）/（Tmax-Topt））*（（Tmax-x）/（Tmax-Topt））^ b1   )) }

- Juanchi

不，上面的模型仅限制函数在范围（Tmin，Tmax）之外返回0。它根本不限制Tmin或Tmax。你所做的是在所选参数空间中找到最小的RSS（或多或少，这是一个非常粗略的网格）。从RSS的角度来看，这是“最佳拟合”，但你应该知道，当你以这种方式进行时，拟合的统计数据（参数的se值等）完全没有意义。 - jlhoward

此外，所有的技术都讲述了关于 Tmin 的同样故事 - 在给定数据和模型的情况下，无法准确估计。在另一个答案中的最后一个模型中，Tmin ~ -52 +/- 360，因此可能是任何值。然而，请注意，不同的方法给出了非常不同的 b1 估计值。 - jlhoward

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jlhoward · Accepted Answer

这里有很多问题，我怀疑在一个SO帖子中无法充分涵盖，但这应该可以让你开始着手解决。

首先，看起来你想要Tmax < max(dat$x)，例如，< 35。这会导致问题，因为对于一些值x，Tmax - x < 0，当你试图将负数提高到幂（在公式的第二项中）时，你会得到NA。这就是错误消息的原因。

其次，非线性模型的收敛取决于模型公式和数据，所以过程在一个数据集上收敛而在另一个数据集上不收敛是完全无关紧要的。

第三，非线性建模迭代地将残差平方和作为参数的函数最小化。如果RSS曲面具有局部极小值，并且你的start接近其中一个局部极小值，算法将找到它。但只有全局最小值才是真正的解。你的问题有很多，很多局部极小值。

第四，nls（...）默认使用Gauss-Newton方法。Gauss-Newton在移动参数（添加到或从预测器中减去的参数，因此在你的情况下是Tmin和Tmax）时是出了名的不稳定。幸运的是，minpak.lm包实现了Levenberg Marquardt方法，在这些条件下更加稳定。该包中的nlsLM（...）函数使用与nls（...）相同的调用序列，并返回一个类型为nls的对象，因此该类对象的所有方法也适用。使用它。

第五，非线性回归（事实上所有最小二乘回归）的一个基本假设是残差服从正态分布。因此，必须使用Q-Q图验证任何解决方案。

第六，你的模型具有一组反常的特征。当Tmin -> -Inf时，模型中的第一项趋近于1。结果，比min(dat$x)小的任何其他Tmin值都会产生较低的RSS，因此所有算法都倾向于将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))

这看起来非常适合，但事实并非如此：Q-Q图表明残差根本不是正态分布的。 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))

从数学上讲，这是一个明显更好的拟合：残差平方和更小且残差分布更加接近正态分布。同样的，参数不太容易估计且没有物理意义，这是数据（或者也可能是模型公式）的问题，而非拟合过程的问题。

所有这些都表明你的模型存在一些问题。一个数学上的问题是当 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的估计要好得多（与早期模型函数的估计值非常不同）。残差仍然不正常，但在这种情况下，我想检查数据是否存在异常值。