Search code examples
rlistnlsnlme

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


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?


Solution

  • 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