将nls拟合成组数据R

Question

我试图将非线性模型与整个季节几个地块上收集的一系列测量数据进行拟合。下面是来自较大数据集的子示例。数据：

dput(nee.example)结构(list=c(159 L、159 L、169 L、169 L、159 L、169 L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L、2L)、.Label = c("e“、”w“、”类“)、类型=结构(c( 1L、1L、1L、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、2L)、.Label = c("b”、"g")，( 1L、1L、1L、2L、2L、2L、2L、2L、3L、3L、3L、3L、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、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、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、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、1L、1L、1L、1L、1L、1L、1L、1L1L、1L、1L、1L、1L、1L、.Label = "a“、”因子“)、布匹= c(25L、50L、75L、100 L、0L、0L、25L、50L、75L、0L、25L、25L、75L、100L、25L、50L、100L、25L、50L、75L、100L )、plotID = c(1L、1L、1L、1L、2L、2L、25L、3L、25L、50L、75L、100L)，plotID=c(1L、1L、1L、1L、2L、2L、2L、3L、3L，3L，3L，3L，13L，13L，13L，13L，14L，14L，14L，14L，15L，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、24.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，( .Names =c(“朱利安”、"blk“、”类型“、”情节“、"trt”、“布”、"plotID“、”通量“、"ChT”、"par")、类= "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中允许的分组结构。换句话说，我可以在朱利安中嵌套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 <0，则会增加)。出于好奇，我绘制了一个导致nls失败的数据集：

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

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

​

​

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