Установка nls для сгруппированных данных R

Question

Я пытаюсь подгонять нелинейную модель к серии измерений, собранных на нескольких участках в течение всего сезона. Ниже приведена подвыборка из более крупного набора данных. данные:

dput (nee.example) структура (список (Julian = с (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, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2 L), .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), .Label = "a", class = "factor"), cloth = 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, 3L, 3L, 3L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L ), fl (0,76, 0,6, 0,67, 0,7, 1,72, 1,63, -7,8, 0,92, -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", "plot", "trt", "cloth", "plotID", "flux", "ChT", "par "), class =" data.frame ", row.names = с (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 внутри julian? Может ли это быть причиной моей ошибки.

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, вы можете легко построить сотни из них, и даже если вы быстро взглянете на них, это поможет вам избежать неприятностей.