“评估模型时产生的缺失值或无穷大”答案

【问题标题】："Missing value or an infinity produced when evaluating the model"“评估模型时产生的缺失值或无穷大”
【发布时间】：2018-11-01 01:52:05
【问题描述】：

我正在尝试使用 nls 函数解决 R 中的两分量衰减模型，但遇到了错误。等式是：

其中t是时间，Ctot是C1+C2，p1和p2是Ctot的已知比例。

我的数据 (dd) 是：

> head(dd,n=15)
      t    Ctot
1  0.00 6.62
2  0.33 6.45
3  0.50 6.38
4  0.67 6.44
5  0.83 6.38
6  1.00 6.39
7  1.17 6.35
8  1.33 6.33
9  1.50 6.33
10 1.67 6.28
11 1.83 6.17
12 2.00 6.11
13 2.17 6.07
14 2.33 5.89
15 2.50 5.86

使用我尝试过的 nls：

p1 <- 0.3
p2 <- 0.7   
z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01))

但是我得到了：

z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01))
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

如果有人有建议，我将不胜感激！

【问题讨论】：

可能k1和k2的最优值比较接近，所以只有一个term的模型会更好。
@G.Grothendieck 是的，但它是一个双分量衰减模型。
如果 k1 = k2 则不是。
似乎 p 和 C 可能无法单独估计。如果Ctot == C1+C2，那么p1和p2不是必须为1或者是固定比例吗？
我猜Ctot==C1+C2 仅适用于t==t0。否则C1 和C2 不能是常量，因为Ctot 不会随着时间的推移而保持不变。

标签： r nls

【解决方案1】：

数据似乎相当有限并且显然不完整，因为它只是头部。如果我们为测试方法编造一些数据……而忽略令人困惑的 p1 和 p2：

 t=seq(0, 20, by=.3)
 Ctot = 3 * exp( -1 * t) + 4 * exp(-5*t)
 # following hte example on gnm::gnm's help page:

saved.fits <- list(); library(gnm)
for (i in 1:10) {
      saved.fits[[i]] <- suppressWarnings( gnm(Ctot ~ Exp(1 + t, inst = 1) +
                                                      Exp(1 + t, inst = 2),
                    verbose=FALSE))}
plot(Ctot~t)
lines(saved.fits[[3]]$fitted~t)
lines(saved.fits[[3]]$fitted~t,col="red")

我不熟悉 gnm 包，因此最终阅读了前几节，然后阅读了其小插图中的工作 2 组件数据拟合示例：https://cran.r-project.org/web/packages/gnm/vignettes/gnmOverview.pdf。大多数拟合都符合预期，但有些人会发现局部最大值的可能性不是全局最大值：

> saved.fits[[1]]$coefficients
                     (Intercept) Exp(. + t, inst = 1).(Intercept) 
                    1.479909e-12                     1.098612e+00 
          Exp(1 + ., inst = 1).t Exp(. + t, inst = 2).(Intercept) 
                   -1.000000e+00                     1.386294e+00 
          Exp(1 + ., inst = 2).t 
                   -5.000000e+00 
attr(,"eliminated")
[1] 0
> exp( saved.fits[[1]]$coefficients[4] )
Exp(. + t, inst = 2).(Intercept) 
                               4 
> exp( saved.fits[[1]]$coefficients[2] )
Exp(. + t, inst = 1).(Intercept) 
                               3

【讨论】：

【解决方案2】：

问题中显示的数据似乎效果不佳，但如果您对其他参数模型持开放态度，那么这个 3 参数模型似乎是合理的。

fm <- nls(Ctot ~ 1 / (a + b * t^c), dd, st = list(a = 1, b = 1, c = 1))
plot(dd)
lines(fitted(fm) ~ t, dd, col = "red")

【讨论】：