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

`nlme::nlsList` doesn't support different algorithms. – Roland Dec 19 '14 at 17:43 — Roland, Dec 19 '14 at 17:43

Tom Wenseleers · Answer 1 · 2014-12-23T20:57:31.900

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

Hi Tom, Sorry to hear you find yourself in the same situation, but good to hear that I'm not alone in this. Would indeed be good if this could be solved at source code level. In the mean time: Did you get your workaround working? — renske hijbeek, Jan 06 '15 at 12:47
Hi Tom, I have adjusted your code to use it as a loop and store the coefficients in a dataframe (so without the trouble of making a nlslist function - which would be a nicer solution I think, but can't figure it out). Works good enough for me on the moment, although not a long term solution. Thanks for your help! I can check your answer as the solution, but maybe you want to adjust to a workable one? — renske hijbeek, Jan 07 '15 at 13:04
Hi Renske - Glad I could help a little - well maybe leave my post as it is and wait a little with accepting an answer, maybe someone sooner or later will post a solutions to show how to properly construct a working `nlsList` object, or fix the `nlsList` source code, don't know... — Tom Wenseleers, Jan 07 '15 at 13:24
Ok, I'll leave it like this then and wait to see if anybody comes up with a suggestion to fix nlslist — renske hijbeek, Jan 07 '15 at 13:44

score 1 · Answer 2 · answered Jan 07 '15 at 13:01

Adjusting Tom's suggestion for making a workaround, I am now using the following code:

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

Which gives me the coefficients of the fits:

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

Works good enough for me on the moment, also for larger datasets.

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

2 Answers2