Use of "algorithm = 'port' and control lower limit in nlsList

Question

I am using the data from How can I get the coefficients from nlsList into a dataframe?

library(nlme)
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)
long <- tidyr::pivot_longer(dat, -1, values_to = "y")
long$name <- factor(long$name)
st0 <- list(Max = 200, k = 0.1, Lag = 0.5)
nlsList(y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))) | name,
 long, 
algorithm="port",
lower=c(k = 0.1, Max =-Inf, Lag = -Inf), 
start = st0)

What I need differently is to not have k lower than 0.1, so I used algorithm="port", lower=c(k = 0.1, Max =-Inf, Lag = -Inf) as in nls() Prevent a nls-fit from falling below zero. It doesn't look like nlsList is taking those 2 commands.

Error in nlsList(y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))) | name,  : 
  unused arguments (algorithm = "port", lower=c(k = 0.1, Max =-Inf, Lag = -Inf))

How do I work around this issue?

Ben Bolker · Accepted Answer · 2022-10-14T18:59:32.370

It looks like nlsList doesn't take those additional arguments — as far as I can tell that's just an oversight on the part of the authors. (This might be worth a feature/enhancement request on the R bug tracker ...)

In the meantime you can use tidyverse as here to split-apply-modify-combine ...

Note that the values in the lower argument must be specified in the same order as start: names are silently ignored ... (Submitted to the R bug tracker here ...).

models <- (long
    |> group_by(name)
    |> nest()
    |> mutate(fit = map(data,
                        nls,
                        form = y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))),
                        algorithm="port",
                        lower=c(Max =-Inf, k = 0.1,  Lag = -Inf), 
                        start = st0))
)

coefs <- (models
    |> mutate(cc = map(fit, broom::tidy))
    |> select(name, cc)
    |> unnest(cols = cc)
)

  name   term  estimate std.error statistic    p.value
  <fct>  <chr>    <dbl>     <dbl>     <dbl>      <dbl>
1 gluc   Max    300.      16.7        18.0  0.0000561 
2 gluc   k        0.162    0.0382      4.23 0.0134    
3 gluc   Lag      2.43     0.515       4.71 0.00924   
4 starch Max    357.      11.8        30.1  0.00000722
5 starch k        0.161    0.0211      7.64 0.00157   
6 starch Lag      1.80     0.234       7.70 0.00153   
7 solka  Max    317.      20.0        15.8  0.0000929 
8 solka  k        0.1      0.0321      3.12 0.0356    
9 solka  Lag      7.61     1.53        4.98 0.00758

This may not be a problem for the poster but note that this effectively is using `pool=FALSE` whereas the code in the question is using the default which is `pool=TRUE`. — G. Grothendieck, Oct 12 '22 at 12:21
I've just noticed that k = 0.06 for solka. Looks like the lower limit was ignored. — hnguyen, Oct 12 '22 at 13:09

Fernando Miguez · Answer 2 · 2022-10-14T21:03:16.090

2

One suggestion is to use nlsLMList in package nlraa. It will use the algorithm in 'minpack.lm::nlsLM'. Notice that I reorder the names in the lower argument. (You need the latest version from github for this to work).

library(nlme)
library(nlraa)
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)
long <- tidyr::pivot_longer(dat, -1, values_to = "y")
long$name <- factor(long$name)
st0 <- list(Max = 200, k = 0.1, Lag = 0.5)
nlsLMList(y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))) | name,
 long, 
algorithm="port",
lower=c(Max =-Inf, k = 0.1, Lag = -Inf), 
start = st0)

However, a better approach is to fit an NLME

longG <- groupedData(y ~ time | name, data = long)

fitL <- nlsLMList(y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))),
                  longG, 
                  start = st0)

fm <- nlme(fitL, random = pdDiag(Lag + Max + k ~ 1))

plot(augPred(fm, level = 0:1))

coef(fm)

If a fixed-effects only model is preferred then this might work

fgm <- gnls(y ~ (time > Lag) * Max * (1-exp(-k * (time - Lag))),
            data = long,
            params = Lag + Max + k ~ name,
            start = c(Lag = 2.4, 0, 0, Max = 336, 0, 0, k = 0.12, 0, 0))

edited Oct 14 '22 at 21:03

answered Oct 14 '22 at 19:38

Fernando Miguez

21
2

A mixed model may or may not be better (if there are only three levels of the grouping variable it might not be very useful; if there is lots of data per level, and you *don't* want shrinkage estimators, a split fixed-effect model might be better) – Ben Bolker Oct 14 '22 at 20:00
Good point! I edited the response with a fixed-effects model – Fernando Miguez Oct 14 '22 at 21:03
But neither `gnls` nor `nlme` allow bounds on parameters (AFAICS), which was the primary point (I think) of the original question ... – Ben Bolker Oct 14 '22 at 21:05
Yes, you are right! I was thinking that if the concern is for a parameter estimate to be below zero (when this is not expected) then the NLME - with shrinkage - would allow to estimate that parameter better by combining the different groups. I guess an approach often used is to re-express the model so that the estimate is forced to be positive as done in a number of SS functions. – Fernando Miguez Oct 14 '22 at 21:12

Use of "algorithm = 'port' and control lower limit in nlsList

2 Answers2