问题描述
我正在为循环 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
?
解决方法
构建要用作字符串的公式:
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"]])
在下面的所有替代方案中,我们使用这些值:
long <- tidyr::pivot_longer(dat,-1,values_to = "y")
long$name <- factor(long$name)
st0 <- list(Max = 200,Lag = 0.5)
1) nls 分组数据 将 dat 转换为长格式,然后使用 nls
的分组数据特性 此方案最适合这里介绍的用于测试某些参数是否通用的方案在三个名称之间,因为如果要在名称中通用,很容易简单地删除参数上的下标。拟合本身不使用任何包,但我们显示 ggplot2 和格子包图形用于绘图。
# 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,Lag = rep(Lag,3)))
model <- nls(y ~ Max[name] * (1-exp(-k[name] * (time - Lag[name]))),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
不必安装)。 long
和 st0
来自上方。
library(nlme)
fit <- nlsList(y ~ Max * (1-exp(-k * (time - Lag))) | name,start = st0)
给出一个 nlsList
对象,其 3 个组件是通过为每个 nls
运行 nls
获得的三个 name
对象。
> 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,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,cex = 1.5)
3) 子集 与(2) 类似的方法是使用nls
执行三个subset=
运行来选择数据。这将返回 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) 中的图形代码也适用于此。