将nls拟合为分组数据R.

Question

我试图将非线性模型拟合到整个季节中几个地块上收集的一系列测量数据中.以下是较大数据集的子样本.数据:

dput(nee.example)结构(列表(julian = 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),.标签= 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 ="factor"),plot = 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),.标签="a",类="因子"),布= 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)),. Name = c("julian","blk","type","plot","trt","cloth","plotID", "flux","ChT","par"),class ="data.frame",row.names = c(NA,-29L))"flux","ChT","par"),class ="data.frame",row.names = c(NA,-29L))"flux","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 <0则增加).出于好奇,我绘制了一个导致nls失败的数据集:

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

并发现它不是单调的,我甚至会说它根本没有任何趋势!我想即使你强行为nls这个案子找到一些解决方案,也很可能是一个虚假的解决方案,所以你可能只想让它不合适(即NA).

我还建议对输入数据和拟合模型质量进行目视检查.随着R和像包reshape2和ggplot2,你可以很容易地绘制数百人,甚至采取了快速看一下他们会帮助你远离麻烦.