R如何在nlslist中使用边界和 "port"算法？

Question

我想在许多数据集上拟合曲线，并按治疗分组。这与 nlslist 配合得很好，但现在我想引入参数的上限。

当我用 nls 分别拟合每个组时，引入边界确实效果很好，但当我想用 nlslist 加速我的工作(我有更多的数据集)时，显然效果不佳。

有人可以帮我解决这个问题吗？

我的数据集示例:

DF1<-data.frame(treatment = rep(c("mineral","residues"),4),
            N_level = c(0,0,100,100,200,200,300,300),
            yield = c(8,8.5,10,10.5,11,9.8,9.5,9.7))

DF1

    treatment N_level yield
1   mineral       0   8.0
2  residues       0   8.5
3   mineral     100  10.0
4  residues     100  10.5
5   mineral     200  11.0
6  residues     200   9.8
7   mineral     300   9.5
8  residues     300   9.7

尝试仅使用 nls 来拟合此数据集效果很好:

fit_mineral <- nls(formula = yield ~ a + b*0.99^N_level +c*N_level, 
                data=subset(DF1, subset = treatment == "mineral"), 
                algorithm = "port", start = list(a = 12, b = -8, c= -0.01), 
                upper = list(a=1000, b=-0.000001, c=-0.000001))

fit_mineral

Nonlinear regression model
model: yield ~ a + b * 0.99^N_level + c * N_level
data: subset(DF1, subset = treatment == "mineral")
    a       b       c 
13.7882 -5.8685 -0.0126 
residual sum-of-squares: 0.4679

但是当我尝试在 nlslist 中组合内容时，它就不起作用:

fit_mineral_and_residues <- nlsList(model = yield ~ a + b*0.99^N_level +c*N_level 
                      | treatment, data=DF1, 
                      algorithm = "port", start = list(a = 12, b = -8, c= -0.01), 
                      upper = list(a=1000, b=-0.000001, c=-0.000001))

错误信息:

Error in nlsList(model = yield ~ a + b * 0.99^N_level + c * N_level |  : 
unused arguments (algorithm = "port", upper = list(a = 1000, b = -1e-06, c = -1e-06))

Answer 1

我刚刚遇到了同样的问题 - 我相信这个问题确实应该在源代码级别解决!

作为一种解决方法，您也许可以尝试自己构造 nlsList 对象，类似于

library(nlme)
DF1=data.frame(treatment = rep(c("mineral","residues"),4),
               N_level = c(0,0,100,100,200,200,300,300),
               yield = c(8,8.5,10,10.5,11,9.8,9.5,9.7))
nlslist=lapply(unique(DF1$treatment),function(i) {datasubs=DF1[DF1$treatment==i,];
                                                             nls(yield ~ a + b*0.99^N_level +c*N_level, 
                                                                 data=datasubs, 
                                                                 start = list(a = 12, b = -8, c= -0.01), 
                                                                 upper = list(a=1000, b=-0.000001, c=-0.000001), 
                                                                 algorithm="port", 
                                                                 control=list(maxit=100000,tol=1e-10,warnOnly=T,minFactor=1e-10) )
                                                  })
names(nlslist)=unique(DF1$treatment)
attr(nlslist, "dims")=list(N = nrow(DF1), M = length(nlslist))
attr(nlslist, "call")=NA # this line is not correct - should be fixed
attr(nlslist,"groups")=names(nlslist)
attr(nlslist, "origOrder")=1:length(unique(DF1$treatment))
attr(nlslist, "pool")=TRUE
attr(nlslist, "groupsForm")=formula(~treatment)
class(nlslist)=c("nlsList", "lmList")

这几乎让我到达那里，除了我不太知道如何正确填写 "call" 槽(在 nlsList 中，它是使用 构造的>match.call() - 有人知道如何做到这一点吗？

如果您想检查正确的结构，可以通过查看例如在

test=nlsList(uptake ~ SSasympOff(conc, Asym, lrc, c0),
               data = CO2, start = c(Asym = 30, lrc = -4.5, c0 = 52))
class(test)=NULL
test