filter() 在 For 循环中不起作用的解决方案？

Question

我正在尝试将 r 中的 von Bertalanffy 增长函数 (VGBF) 拟合到按序列号分组的数据中。 这是我的数据的 sn-p：

Serial_No<- c(315,315,315,315,315,315,315,316,316,316,316,317,317,317,317,317,317,317,317,317,318,318,318,318,319,319,319,319)

Year<-c(1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945)

tl<-c(19,33,46,55,63,66,70,22,39,55,65,20,40,51,56,60,62,63,64,65,26,43,54,60,28,43,53,61)

age<-c(1,2,3,4,5,6,7,1,2,3,4,1,2,3,4,5,6,7,8,9,1,2,3,4,1,2,3,4))

df<-data.frame(Serial_No, Year, tl, age)

我一直在关注这个例子：https://www.r-bloggers.com/2020/01/von-bertalanffy-growth-plots-ii/ 并将我的代码更改为如下：

vb <- vbFuns()
predict2 <- function(x) predict(x,data.frame(age=ages))

agesum <- group_by(df,Serial_No) %>%
summarize(minage=min(age),maxage=max(age))

Serial_Nos <- unique(df$Serial_No)
nSerial_Nos <- length(Serial_Nos)

cfs <- cis <- preds1 <- preds2 <- NULL

for (i in 1:nSerial_Nos) {
 cat(Serial_Nos[i],"Loop\n")
 tmp1 <- filter(df,Serial_No==Serial_Nos[i])
 sv1 <- vbStarts(tl~age,data=tmp1)
 fit1 <- nls(tl~vb(age,Linf,K,t0),data=tmp1,start=sv1)
 cfs <- rbind(cfs,coef(fit1))
 boot1 <- Boot(fit1)
 tmp2 <-  confint(boot1)
 cis <- rbind(cis,c(tmp2["Linf",],tmp2["K",],tmp2["t0",]))
 ages <- seq(-1,16,0.2)
 boot2 <- Boot(fit1,f=predict2)
 tmp2 <- data.frame(Serial_No=Serial_Nos[i],age=ages,
                 predict(fit1,data.frame(age=ages)),
                 confint(boot2))
 preds1 <- rbind(preds1,tmp2)
 tmp2 <- filter(tmp2,age>=agesum$minage[i],age<=agesum$maxage[i])
 preds2 <- rbind(preds2,tmp2)
}

代码运行，但从 VBGF 返回的结果对于每个序列号都是相同的，这是不可能的。我认为是过滤器功能在上面的代码中不起作用。

我已经搜索了解决方案，但无法让它发挥作用。

如果有人可以帮忙，或者知道解决方案，我将不胜感激。

提前谢谢你

Answer 1

模型适合包growthrates

以下帖子描述了一种没有 for 循环和 filter 的替代方法。类似的无循环解决方案可以使用常见的nls函数和“base”R中的lapply或“tidyverse”中的group_by来实现。

模型定义

growthrates 包不包含 von Bertalanffy 函数，因此它必须作为用户提供的模型提供，如 package vignette 中所述。这里我借用了 FSA 包中的功能并相应地对其进行了调整：

library("growthrates")

grow_von_bert <- function(time, parms) {
  with(as.list(parms), {
    y <- Linf * (1 - exp(-K * (time - t0)))
    as.matrix(data.frame(time = time, y = y))
  })
}

用一个例子测试模型

p <- c(t0=5, Linf=10, K=.1)
time <- seq(5, 100)
plot(grow_von_bert(time, p), type="l")

单个数据示例的拟合

在对所有人进行此操作之前，先拟合一个或多个单个示例总是一个好主意。

df1 <- subset(df, Serial_No == 315)
fit1 <- fit_growthmodel(df1$age, df1$tl,
                        FUN = grow_von_bert,  p=c(t0=0, Linf=70, K=0.1))
summary(fit1)

拟合所有数据集

这可以在循环中或使用适当的tidyverse 函数来完成，wheate 包 growthrates 已经内置了这样的函数，因此所有模型都可以安装一个函数调用。当然，有必要指定良好的起始参数，无论是对所有曲线还是单个参数集都相同，具体取决于数据的质量。这是包含OP数据的完整代码：

library("growthrates")

df <- data.frame(
  Serial_No = factor(c(315,315,315,315,315,315,315,316,316,316,316,317,317,317,317,
                       317,317,317,317,317,318,318,318,318,319,319,319,319)),
  year = c(1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,
          1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945,1945),
  tl = c(19,33,46,55,63,66,70,22,39,55,65,20,40,51,56,60,62,63,64,65,26,43,54,60,28,
         43,53,61),
  age = c(1,2,3,4,5,6,7,1,2,3,4,1,2,3,4,5,6,7,8,9,1,2,3,4,1,2,3,4)
)

grow_von_bert <- function(time, parms) {
  with(as.list(parms), {
    y <- Linf * (1 - exp(-K * (time - t0)))
    as.matrix(data.frame(time = time, y = y))
  })
}


fit <- all_growthmodels(tl ~ age | Serial_No, 
                        data=df, 
                        FUN = grow_von_bert,
                        p=c(t0=0, Linf=70, K=0.1))

results(fit)
par(mfrow=c(2,3))
plot(fit, las=1)