将 nls 拟合到分组数据 R

Question

我正在尝试将非线性模型拟合到整个季节在多个地块上收集的一系列测量值。下面是来自较大数据集的子样本。 数据：

输入（nee.example） 结构（列表（朱利安 = c（159L，159L，159L，159L，159L，159L， 159L, 159L, 159L, 159L, 159L, 159L, 159L, 159L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L, 169L), blk = 结构(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("e", "w"), class= "factor"), type = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("b", "g"), class= "因子"), 绘图 = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L，2L，3L，3L，3L，3L，3L，1L，1L，1L，1L，2L，2L，2L，2L，2L， 3L, 3L, 3L, 3L, 3L, 3L), trt = 结构(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "a", class= "factor"), 布 = c(25L, 50L, 75L, 100L, 0L, 25L, 50L, 75L, 100L, 0L, 25L, 50L, 75L, 100L, 0L, 25L, 50L, 100L, 0L, 25L, 50L, 75L, 100L, 0L, 25L, 50L, 75L, 75L, 100L), plotID = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L ), 通量 = c(0.76, 0.6, 0.67, 0.7, 1.72, 1.63, -7.8, 0.89, 0.51、0.76、0.48、0.62、0.18、0.21、3.87、2.44、1.26、-1.39、 2.18、1.9、0.81、-0.04、-0.83、1.99、1.55、0.57、-0.02、-0.16、 -2.12), ChT = c(18.6, 19.1, 19.6, 19.1, 16.5, 17.3, 18.3, 19、18.6、17.2、18.4、19、19.2、20.6、22、21.9、22.4、23.8、 20.7、21.5、22.5、23.3、23.8、20.1、20.8、21.2、21.8、21.8、 21.4), par = c(129.9, 210.2, 305.4, 796.6, 1.3, 62.7, 149.9, 171.2、453.3、1.3、129.7、409.3、610、1148.6、1.3、115.2、 237、814.6、1.3、105.4、293.4、472.1、955.9、1.3、100.5、 290, 467, 413.6, 934.2)), .Names = c("julian", "blk", "type", “情节”，“trt”，“布料”，“plotID”，“通量”，“ChT”，“par”），class= “data.frame”，row.names = c（NA， -29L))

我需要将以下模型（rec.hyp，如下）拟合到每个日期的每个图，并检索每个 julian-plotID 组合的参数估计值。经过一番摸索，听起来 nlsList 将是一个理想的函数，因为它具有分组方面：

library(nlme)
rec.hyp <- nlsList(flux ~ Re - ((Amax*par)/(k+par)) | julian/plotID,
             data=nee.example,
             start=c(Re=3, k=300, Amax=5),
             na.action=na.omit)
coef(rec.hyp)

但是我不断收到相同的错误消息：

Error in nls(formula = formula, data = data, start = start, control = control) : 
step factor 0.000488281 reduced below 'minFactor' of 0.000976562

我尝试调整 nls.control 中的控件以增加 maxIter 和 tol，但仍显示相同的错误消息。而且我已经更改了初始起始值，但无济于事。

需要注意的是，为了与之前的工作保持一致，我需要使用最小二乘来拟合模型。

问题：

1. 在 nlsList 中是否允许我的分组结构。换句话说，我可以在 julian 中嵌套 plotID 吗？这可能是我错误的根源吗？

2. 我已经读到不适当的起始参数估计会导致错误消息，但在更改它们后我得到相同的消息。

我觉得我在这里遗漏了一些简单的东西，但我的大脑被炸了。

提前致谢。

Answer 1

Q1 的答案：您的分组结构是正确的。您可以通过在数据子集上运行 nls 来验证它：

rec.hyp.test <- nls(flux ~ Re - ((Amax*par)/(k+par)),
                   data=subset(nee.example,julian==159 & plotID==3),
                   start=c(Re=3, k=300, Amax=5),
                   na.action=na.omit)
coef(rec.hyp.test)
#        Re           k        Amax 
# 0.7208943 792.4412287   0.8972519 

coef(rec.hyp)[3,]
#              Re        k      Amax
# 159/3 0.7208943 792.4412 0.8972519

对 Q2 的回答：某些数据集无法正确拟合给定模型。从flux ~ Re - ((Amax*par)/(k+par)) 公式中，可以预期flux 随着par 单调减少（或增加，如果Amax nls 失败的数据集：

plot(flux~par,subset(nee.example,julian==159 & plotID==1)) 

发现它不是单调的，我什至会说它根本没有任何趋势！我想即使你强制 nls 为这种情况找到一些解决方案，它也很可能是一个虚假的解决方案，所以你可能只想让它不合适（即 NA）。

我还建议对输入数据和拟合模型质量进行目视检查。使用R 和reshape2 和ggplot2 之类的软件包，您可以轻松绘制数百个，甚至快速查看它们也可以帮助您避免麻烦。