5

我想在许多数据集上拟合曲线，按治疗分组。这与 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))

适合矿物

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))
4

2 回答 2

3

我刚刚遇到了同样的问题——我相信这个问题真的应该在源代码级别解决！

作为一种解决方法，您也许可以尝试自己构建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
于 2014-12-23T20:29:08.150 回答
1

调整汤姆提出的解决方法的建议，我现在使用以下代码：

DF2 <- lapply(unique(DF1$treatment),function(i) {datasubs=DF1[DF1$treatment==i,];
                                              coef(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)))
})

DF2 <- data.frame(DF2)
names(DF2) <- levels(DF1$treatment)
DF2 <- t(DF2)
DF2

这给了我拟合的系数：

                a         b           c
mineral  13.78825 -5.868506 -0.01260393
residues 12.76157 -4.201832 -0.01006236

目前对我来说足够好，也适用于更大的数据集。

于 2015-01-07T13:01:55.963 回答