对于此数据集：
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
最佳答案
这里有太多问题，我怀疑它是否可以在SO帖子中充分介绍，但这应该可以帮助您入门。
首先，看起来像您想要Tmax < max(dat$x)，例如Tmax - x < 0对于x的某些值以及当您尝试将负数提高为幂时（在第二项中，您的公式），则得到NA。这是错误消息的原因。
其次，非线性模型的收敛性取决于模型公式以及数据，因此该过程与一组数据收敛而与另一组数据收敛这一事实是完全不相关的。
第三，非线性建模根据参数迭代地最小化残差平方和。如果RSS表面具有局部最小值，并且您的start接近于最小值，则算法会找到它。但是，只有全局最小值是真正的解决方案。您的问题有很多很多局部最小值。
第四，nls(...)默认情况下使用高斯牛顿法。众所周知，高斯牛顿具有不稳定的参数（参数要添加到预测变量中或从预测变量中减去，因此在您的情况下为Tmin和Tmax）是不稳定的。幸运的是，minpak.lm包实现了Levenberg Marquardt方法，该方法在这些条件下更加稳定。该包中的nlsLM(...)函数使用与nls(...)相同的调用顺序，并返回类型为nls的对象，因此该类对象的所有方法也都可以使用。用那个
第五，非线性回归（实际上是所有最小二乘回归）的基本假设是残差呈正态分布。因此，您必须使用Q-Q图来验证任何解决方案。
第六，您的模型具有一组不正确的特征。当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))
这看起来很合适，但事实并非如此：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))
从数学上讲，这是一个非常优越的拟合：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的估计要好得多（并且与早期模型函数的估计有很大不同）。残差仍然不正常，但是在这种情况下，我想检查数据是否存在异常值。