將nls與分組數據擬合R

Question

我試圖將一個非線性模型擬合到整個賽季中幾個地塊上收集的一系列測量數據中。以下是來自較大數據集的子樣本。 數據：

dput（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，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，2 L），標籤= 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 = structure（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），.Label =「a」，class =「factor」）， cloth = c（25L，50L，75L，100L，0L，25L，50L，75L，100L，0L， 25L， ，75L，100L，0L，25L，50L，100L，0L，25L，50L，75L， 100L，0L，25L，50L，75L，75L，100L），繪圖ID = 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， ），fl ux = c（0.76,0.6,0.67,0.7,1.72,1.63，-7.8,0.89,0.）0.51,0.76,0.48,0.62,0.18,0.21,3.87,2.44,1.26，-1.39,0.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， （129.9,210.2,305.4,796.6），其中所述標準化合物的特徵在於，所述化合物的特徵在於： ，1.3,6.2.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， ），「N」= c（「julian」，「blk」，「type」， 「plot」，「trt」，「cloth」，「plotID」，「flux」，「ChT」 「），class =」data.frame「， row.names = C（NA， -29L））

我需要適合下面的模型（rec.hyp，下文），以在每個日期的每個情節和檢索每個儒略-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等包，您可以輕鬆地繪製數百個包，甚至快速查看它們將幫助您避免麻煩。