R how to use bounds and the "port" algorithm in nlslist?

Question

I would like to fit curves on a number of data sets, grouped by treatment. This was working very well with nlslist, but now I would like to introduce upper bounds to my parameters.

Introducing bounds does work out very nicely when I fit each group seperately with nls, but apparently not when I want to speed up my work (I have many more data sets) with nlslist.

Could anybody help me out here with any idea how to solve this?

An example of my dataset:

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

Trying to fit this dataset with only nls works well:

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

But as soon as I try to combine things in nlslist, it doesnt work:

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 message:

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

I just encountered the same problem - this issue I believe really should be fixed at source code level!

As a workaround you could perhaps try to construct the nlsList object yourself, something along the lines of

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

This almost gets me there, except that I don't quite know how to correctly fill in the "call" slot (in nlsList it is constructed using match.call() - anybody that knows how to do this by any chance?

If you would like to check the correct structure this can be done by looking e.g. at

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