NLS регрессия и хранение выходных коэффициентов и участков

Question

я пытался сделать воспроизводимые данные, как было предложено How to make a great R reproducible example?:

structure(list(ID = structure(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, 2L), .Label = c("ANK.26.1", 
"ANK.35.10"), class = "factor"), DAY = c(2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 
18L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
15L, 16L, 17L, 18L), carbon = c(1684.351094778, 3514.451339358, 
6635.877888654, 10301.700591252, 11361.360992769, 11891.934331254, 
12772.885869486, 13545.127224369, 14022.00520767, 14255.045990397, 
14479.813468278, 14611.749542181, 14746.382638335, 14942.733363567, 
14961.338739162, 15049.433738817, 15047.197961499, 1705.361701104, 
3293.593040601, 4788.872254899, 6025.622715999, 6670.80499518, 
7150.526272512, 7268.955557607, 7513.61998338, 7896.202773246, 
8017.953574608, 8146.09464786, 8286.148260324, 8251.229520243, 
8384.244997158, 8413.034235219, 8461.066691601, 8269.360979031 
), g.rate.perc = c(NA, 1.08653133557123, 0.888168948119852,0.55242467750436, 
0.102862667394628, 0.0466998046116733, 0.0740797513417739, 0.060459426536321, 
0.0352066079115925, 0.0166196474238596, 0.0157675729725753, 0.00911172469120847, 
0.00921402983026387, 0.0133151790542558, 0.00124511193115184, 
0.00588817626489591, -0.000148562222127446, NA, 0.931316411333049, 
0.45399634862756, 0.258255053647507, 0.107073129133681, 0.0719135513148148, 
0.0165623173150578, 0.0336588143694119, 0.0509185706373581,0.0154189051191185, 
0.0159817679236518, 0.0171927308137518, -0.00421410998016991, 
0.0161206856006937, 0.00343373053515927, 0.00570929049366353, 
-0.0226573929218994), max.carb = c(15049.433738817, 15049.433738817, 
15049.433738817, 15049.433738817, 15049.433738817, 15049.433738817, 
15049.433738817, 15049.433738817, 15049.433738817, 15049.433738817, 
15049.433738817, 15049.433738817, 15049.433738817, 15049.433738817, 
15049.433738817, 15049.433738817, 15049.433738817, 8461.066691601, 
8461.066691601, 8461.066691601, 8461.066691601, 8461.066691601, 
8461.066691601, 8461.066691601, 8461.066691601, 8461.066691601, 
8461.066691601, 8461.066691601, 8461.066691601, 8461.066691601, 
8461.066691601, 8461.066691601, 8461.066691601, 8461.066691601 
)), .Names = c("ID", "DAY", "carbon", "g.rate.perc", "max.carb" 
), row.names = c(NA, 34L), class = "data.frame") 

'data.frame': 34 obs. of 5 variables: 
$ ID   : Factor w/ 150 levels "ANK.26.1","ANK.35.10",..: 1 1 1 1 1 1 1 1 1 1 ... 
$ DAY  : int 2 3 4 5 6 7 8 9 10 11 ... 
$ carbon  : num 1684 3514 6636 10302 11361 ... 
$ g.rate.perc: num NA 1.087 0.888 0.552 0.103 ... 
$ max.carb : num 15049 15049 15049 15049 15049 ... 

В ID выборке данных только имеет два уровня не указанный 150.

Моего NLS выглядит как то:

res.sample <- ddply (
    d, .(ID), 
    function(x){ 
    mod <- nls(carbon~phi1/(1+exp(-(phi2 + phi3 * DAY))), 
     start=list(
      phi1 = x$max.carb, 
      phi2 = int[1], 
      phi3 = mean(x$g.rate.perc)), 
     data=x,trace=TRUE) 
return(coef(mod)) 
} 
) 

phi2 на самом деле результат для перехвата из

int <- coef(lm(DAY~carbon,data=sample)) 

К сожалению, он больше не работает, так как я пытался обернуть его в окружение ddply, но я не могу пройти через все оригинальные 150 уровней идентификатора вручную.

Кроме того, я хотел бы сохранить все три выходных значения для phi1-phi3 в кадре данных/списке с соответствующим идентификатором. Я имел в виду, чтобы сделать это

return(coef(mod)) 

Вишневый на вершине будет участки фактических данных и подогнанной кривой на вершине. Вручную с подмножеством я тоже могу это сделать, но это слишком много времени. Мой уменьшенный код ggplot для этого является

ggplot(data=n, aes(x = DAY, y = carbon))+ 
geom_point(stat="identity", size=2) + 
geom_line(aes(DAY,predict(logMod)))+ 
ggtitle("ID") 

, если каким-то образом идентификатор, который содержит триплет информации, менее полезна, вот как вы можете вернуть его в другой версии

sep_sample <- sample %>% separate(ID, c("algae", "id", "nutrient")) 

я чувствую, это так много, чтобы спросить, но я действительно пробовал, и я могу потратить столько дней на это.

Вот краткое изложение:

мне нужно запустить модель на каждом уровне ID/каждая комбинация водорослей & питательных веществ, если разделить его.

Выходные данные phi должны храниться в виде кадра/списка/таблицы с соответствующей идентификацией того, где они принадлежат.

В идеале есть способ включить ggplots во все это, которые автоматически генерируются и сохраняются.

, как я сказал, модель сама по себе уже работал, но когда я положил в в ddply структуре я получаю следующее сообщение об ошибке:

Error in numericDeriv(form[[3L]], names(ind), env) : 
    Missing value or an infinity produced when evaluating the model 

Я надеюсь, что это то, что вы можете работать с каким-то образом, и это представляется разумным вопросом. Если есть страницы, которые уже предоставляют аналогичное решение, которое я нашел, я с удовольствием посмотрю.

Большое спасибо!

Answer 1

Су я пришел к этому решению, которое не все я хотел, но я думаю, что на один шаг ближе, так как он работает

coef_list <- list() 
curve_list <- list() 
for(i in levels(d$ALGAE)) { 
for(j in levels(d$NUTRIENT)) { 
dat = d[d$ALGAE == i & d$NUTRIENT == j,] 

#int <- coef(lm(DAY~carbon,data=dat)) 

mod <- nls(carbonlog~phi1/(1+exp(-(phi2+phi3*DAY))), 
      start=list(
      phi1=9.364, 
      phi2=0, 
      phi3= 0.135113), 
      data=dat,trace=TRUE) 
coef_list[[paste(i, j, sep = "_")]] = coef(mod) 

plt <- ggplot(data = dat, aes(x = DAY, y = carbonlog)) + geom_point()+ 
    geom_line(aes(DAY,predict(mod)))+ 
    ggtitle(paste(i,"RATIO",j,sep=" ")) + 
    theme.plot 
curve_list[[paste(i, j, sep = "_")]] = plt 
    } 
} 

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

если я применить

curve_list[["ANK_1"]] 

я получаю сообщение об ошибке, хотя:

Error: Aesthetics must be either length 1 or the same as the data (17): x, y 

я только получает сообщение, когда я использовать логарифмические значения углерода. когда я использую углерод в его исходном формате, он отображает все.