带有 nls 的循环函数答案

【问题标题】：loop function with nls带有 nls 的循环函数
【发布时间】：2021-02-04 16:02:40
【问题描述】：

我正在为循环 nls 函数而苦苦挣扎。所以这里是单个样本的示例数据集

dat<-read.table(text="time y
1 4.62
2 13.55
3 30.82
6 93.97
12 145.93
24 179.93", header = TRUE)
plot(data);lines(data)
model <- nls(y ~ Max * (1-exp(-k * (time - Lag))),data=dat,start=list(Max = 200, k = 0.1, Lag = 0.5))

但是如果我想将model 应用于多列样本怎么办？例如

dat<-read.table(text="time gluc starch solka
+ 1 6.32 7.51 1.95
+ 2 20.11 25.49 6.43
+ 3 36.03 47.53 10.39
+ 6 107.52 166.31 27.01
+ 12 259.28 305.19 113.72
+ 24 283.40 342.56 251.14
+ 48 297.55 353.66 314.22", header = TRUE)

我如何让 R 为每个样品（gluc、淀粉、solka）求解 Max、k 和 Lag？

【问题讨论】：

melt data.frame 并使用包 nlme 中的nlsList。

标签： r loops curve-fitting nls

【解决方案1】：

在以下所有备选方案中，我们使用这些值：

long <- tidyr::pivot_longer(dat, -1, values_to = "y")
long$name <- factor(long$name)
st0 <- list(Max = 200, k = 0.1, Lag = 0.5)

1) nls 分组数据 将 dat 转换为长格式，然后使用nls 的分组数据功能这个解决方案是这里介绍的最适合测试某些参数在三个名称，因为如果要在名称中通用，很容易删除参数上的下标。拟合本身不使用任何包，但我们显示 ggplot2 和 lattice 包图形进行绘图。

# get better starting values
model0 <- nls(y ~ Max * (1-exp(-k * (time - Lag))), long, start = st0)  
st <- with(as.list(coef(model0)), 
  list(Max = rep(Max, 3), k = rep(k, 3), Lag = rep(Lag, 3)))

model <- nls(y ~ Max[name] * (1-exp(-k[name] * (time - Lag[name]))), 
  long, start = st)
model

给予：

Nonlinear regression model
  model: y ~ Max[name] * (1 - exp(-k[name] * (time - Lag[name])))
   data: long
     Max1      Max2      Max3        k1        k2        k3      Lag1      Lag2 
306.48737 389.84657 361.82290   0.12214   0.03857   0.13747   1.38072   2.02205 
     Lag3 
  1.31770 
 residual sum-of-squares: 7167

Number of iterations to convergence: 8 
Achieved convergence tolerance: 9.186e-06

ggplot2 图形可以这样完成。

library(ggplot2)
fitdf <- transform(long, fit = fitted(model))
ggplot(fitdf, aes(x = time, y = y, color = name)) +
    geom_point() +
    geom_line(aes(y = fit))

可以使用 R 附带的格子图形生成外观略有不同的图，因此不必安装包。代码特别紧凑。

library(lattice)
xyplot(fit + y ~ time | name, fitdf, type = c("l", "p"), auto.key = TRUE)

2) nlsList 如果您不需要调查名称中参数的常用设置，那么另一种可能性是在 nlme 包中使用 nlsList（它与 R 一起提供，所以您不需要）不必安装它）。 long 和 st0 来自上方。

library(nlme)
fit <- nlsList(y ~ Max * (1-exp(-k * (time - Lag))) | name, long, start = st0)

给出一个nlsList 对象，其3 个组件是通过为每个name 运行nls 获得的三个nls 对象。

> fit
Call:
  Model: y ~ Max * (1 - exp(-k * (time - Lag))) | name 
   Data: long 

Coefficients:
            Max          k      Lag
gluc   306.4875 0.12214330 1.380713
solka  389.8449 0.03856544 2.022057
starch 361.8231 0.13747402 1.317698

Degrees of freedom: 21 total; 12 residual
Residual standard error: 24.43858

我们可以绘制数据并进行拟合：

levs <- levels(long$name)
col <- setNames(rainbow(length(levs)), levs) 
plot(y ~ time, long, col = col[name], pch = 20, cex = 1.5)
for(lv in levs) lines(fitted(fit[[lv]]) ~ time, dat, col = col[lv])
legend("bottomright", leg = levs, col = col, pch = 20, cex = 1.5)

3) 子集 与（2）类似的方法是使用subset= 执行三个nls 运行以选择数据。这将返回nls 对象的命名列表。 st0 和 long 来自上方。没有使用任何包。

fit <- Map(function(nm) nls(y ~ Max * (1-exp(-k * (time - Lag))), data = long, 
    start = st0, subset = name == nm), levels(long$name))

(2) 中的图形代码也适用于此。

【讨论】：

【解决方案2】：

构建要用作字符串的公式：

outcomes = c("gluc", "starch", "solka")
my_formulas = paste(outcomes, "~ Max * (1-exp(-k * (time - Lag)))")
model_list = list()

for(i in seq_along(outcomes)) {
  model_list[[outcomes[i]]] = nls(
    as.formula(my_formulas[i], 
    data = dat,
    start = list(Max = 200, k = 0.1, Lag = 0.5)
  )
}

这将创建一个模型列表，您可以使用例如summary(model_list[[1]]) 或summary(model_list[["solka"]]) 进行访问

【讨论】：