nlme backend does not work correctly with additive and proportional error #428

jranke · 2020-10-29T15:57:01Z

nlme can now do this if you apply the patch at https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17954

times <- c(0, 1, 3, 7, 14, 28, 56, 90, 120)
k_in <- c(0.1, 0.15, 0.16, 0.18, 0.2)
lrc_in <- mean(log(k_in))
n_sample <- length(times) * 2

additive <- 1
proportional <- 0.05

pred <- as.numeric(sapply(k_in, function(k) rep(100 * exp(- k * times), each = 2)))

set.seed(123456L)
d_syn <- data.frame(
  ID = rep(1:5, each = n_sample),
  TIME = rep(times, 5, each = 2),
  AMT = 0,
  EVID = 0,
  CMT = 1,
  DV = rnorm(length(pred), pred, sqrt(additive^2 + pred^2 * proportional^2))
)

library(nlmixr)
# SFO is for simple first-order
sfo_add_prop <- function() {
  ini({
    lK <- 0.1
    lv0 <- log(100)
    eta.lK ~ 0.1
    eta.lv0 ~ 0.1
    add.err <- 1
    prop.err <- 0.05
  })
  model({
    K <- exp(lK + eta.lK)
    value_0 <- exp(lv0 + eta.lv0)

    d/dt(value) = -K * value
    value(0) = value_0

    value ~ add(add.err) + prop(prop.err)
  })
}
f_sfo_add_prop <- nlmixr(sfo_add_prop, data = d_syn, est = "nlme")
#> 
#> **Iteration 1
#> LME step: Loglik: -393.3028, nlminb iterations: 19
#> reStruct  parameters:
#>      ID1      ID2 
#> 5.470245 6.559779 
#> varStruct  parameters:
#>    const 
#> 23.78338 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  7.296427e-18 
#>  fixed effects: -1.866188  4.614154  
#>  iterations: 4 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>        fixed     reStruct    varStruct 
#> 1.0535851581 0.0004293941 1.0968148962 
#> 
#> **Iteration 2
#> LME step: Loglik: -290.0553, nlminb iterations: 1
#> reStruct  parameters:
#>      ID1      ID2 
#> 5.470245 6.559779 
#> varStruct  parameters:
#>    const 
#> 23.78338 
#>  Beginning PNLS step: ..  completed fit_nlme() step.
#> PNLS step: RSS =  7.296427e-18 
#>  fixed effects: -1.866188  4.614154  
#>  iterations: 1 
#> Convergence crit. (must all become <= tolerance = 1e-05):
#>     fixed  reStruct varStruct 
#>         0         0         0
#> Load into sympy...Using sympy via SnakeCharmR
#> done
#> Freeing Python/SymPy memory...done
#> ################################################################################
#> Optimizing expressions in Predictions/EBE model...done
#> Compiling Predictions/EBE model...done
#> Standardized prediction/ebe models produced.
#> Calculating residuals/tables
#> done.
print(f_sfo_add_prop)
#> ── nlmixr nlme by maximum likelihood (ODE; μ-ref & covs) nlme OBF fit ────────── 
#>          OBJF      AIC      BIC Log-likelihood Condition Number
#> nlme 414.7016 592.1105 607.1094      -290.0553         7.472195
#> 
#> ── Time (sec; $time): ────────────────────────────────────────────────────────── 
#>          nlme    setup table    other
#> elapsed 1.393 5.769284 0.028 0.426716
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ─────────────────────────── 
#>              Est.     SE  %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> lK          -1.87 0.0418  2.24    0.155 (0.143, 0.168)        0      100.% 
#> lv0          4.61 0.0153 0.331         101 (97.9, 104)        0      100.% 
#> add.err  2.13e+10                             2.13e+10                     
#> prop.err 2.85e-10                             2.85e-10                     
#>  
#> 
#>   Covariance Type ($covMethod): nlme
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Minimization message ($message):  
#>     Non-positive definite approximate variance-covariance 
#> 
#> ── Fit Data (object is a modified tibble): ───────────────────────────────────── 
#> # A tibble: 90 x 14
#>   ID     TIME    DV  EVID  PRED   RES IPRED  IRES     IWRES    eta.lK  eta.lv0
#>   <fct> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>    <dbl>
#> 1 1         0 104.      0 101.   3.35 101.   3.35  1.57e-10 -1.29e-22 2.00e-23
#> 2 1         0  98.6     0 101.  -2.31 101.  -2.31 -1.08e-10 -1.29e-22 2.00e-23
#> 3 1         1  88.8     0  86.4  2.40  86.4  2.40  1.12e-10 -1.29e-22 2.00e-23
#> # … with 87 more rows, and 3 more variables: K <dbl>, value_0 <dbl>,
#> #   value <dbl>

library(nlme)
f_nlsList <- nlsList(DV ~ SSasymp(TIME, 0, value_0, lK) | ID, d_syn,
  start = list(value_0 = 100, lK = log(0.1)))
f_nlme <- nlme(f_nlsList, weights = varConstProp(), control = list(sigma = 1))
#> Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
#> Iteration 4, LME step: nlminb() did not converge (code = 1). PORT message: false
#> convergence (8)
coef(f_nlme$modelStruct$varStruct, unconstrained = FALSE)
#>      const       prop 
#> 0.93657508 0.04897923

^{Created on 2020-10-29 by the reprex package (v0.3.0)}

Session info

devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.3 (2020-10-10)
#>  os       Debian GNU/Linux 10 (buster)
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language en                          
#>  collate  de_DE.UTF-8                 
#>  ctype    de_DE.UTF-8                 
#>  tz       Europe/Berlin               
#>  date     2020-10-29                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package       * version    date       lib source                              
#>  assertthat      0.2.1      2019-03-21 [1] CRAN (R 4.0.0)                      
#>  backports       1.1.10     2020-09-15 [1] CRAN (R 4.0.2)                      
#>  brew            1.0-6      2011-04-13 [1] CRAN (R 4.0.0)                      
#>  callr           3.5.1      2020-10-13 [1] CRAN (R 4.0.3)                      
#>  cli             2.1.0      2020-10-12 [1] CRAN (R 4.0.3)                      
#>  colorspace      1.4-1      2019-03-18 [1] CRAN (R 4.0.0)                      
#>  crayon          1.3.4      2017-09-16 [1] CRAN (R 4.0.0)                      
#>  data.table      1.13.2     2020-10-19 [1] CRAN (R 4.0.3)                      
#>  desc            1.2.0      2018-05-01 [1] CRAN (R 4.0.0)                      
#>  devtools        2.3.2      2020-09-18 [1] CRAN (R 4.0.2)                      
#>  digest          0.6.27     2020-10-24 [1] CRAN (R 4.0.3)                      
#>  dparser         0.1.8      2017-11-13 [1] CRAN (R 4.0.0)                      
#>  dplyr           1.0.2      2020-08-18 [1] CRAN (R 4.0.2)                      
#>  ellipsis        0.3.1      2020-05-15 [1] CRAN (R 4.0.1)                      
#>  evaluate        0.14       2019-05-28 [1] CRAN (R 4.0.0)                      
#>  fansi           0.4.1      2020-01-08 [1] CRAN (R 4.0.0)                      
#>  fs              1.5.0      2020-07-31 [1] CRAN (R 4.0.2)                      
#>  generics        0.0.2      2018-11-29 [1] CRAN (R 4.0.0)                      
#>  ggplot2         3.3.2      2020-06-19 [1] CRAN (R 4.0.2)                      
#>  glue            1.4.2      2020-08-27 [1] CRAN (R 4.0.2)                      
#>  gtable          0.3.0      2019-03-25 [1] CRAN (R 4.0.0)                      
#>  highr           0.8        2019-03-20 [1] CRAN (R 4.0.0)                      
#>  htmltools       0.5.0      2020-06-16 [1] CRAN (R 4.0.2)                      
#>  jsonlite        1.7.1      2020-09-07 [1] CRAN (R 4.0.2)                      
#>  knitr           1.30       2020-09-22 [1] CRAN (R 4.0.2)                      
#>  lattice         0.20-41    2020-04-02 [1] CRAN (R 4.0.0)                      
#>  lazyeval        0.2.2      2019-03-15 [1] CRAN (R 4.0.0)                      
#>  lifecycle       0.2.0      2020-03-06 [1] CRAN (R 4.0.0)                      
#>  lotri           0.2.2      2020-05-29 [1] CRAN (R 4.0.2)                      
#>  magrittr        1.5        2014-11-22 [1] CRAN (R 4.0.0)                      
#>  Matrix          1.2-18     2019-11-27 [3] CRAN (R 4.0.3)                      
#>  memoise         1.1.0      2017-04-21 [1] CRAN (R 4.0.0)                      
#>  munsell         0.5.0      2018-06-12 [1] CRAN (R 4.0.0)                      
#>  mvnfast         0.2.5.1    2020-10-14 [1] CRAN (R 4.0.3)                      
#>  n1qn1           6.0.1-9    2020-07-02 [1] CRAN (R 4.0.2)                      
#>  nlme          * 3.1-150.1  2020-10-29 [1] local                               
#>  nlmixr        * 1.1.1-9    2020-10-29 [1] local                               
#>  pillar          1.4.6      2020-07-10 [1] CRAN (R 4.0.2)                      
#>  pkgbuild        1.1.0      2020-07-13 [1] CRAN (R 4.0.2)                      
#>  pkgconfig       2.0.3      2019-09-22 [1] CRAN (R 4.0.0)                      
#>  pkgload         1.1.0      2020-05-29 [1] CRAN (R 4.0.1)                      
#>  PreciseSums     0.4        2020-07-10 [1] CRAN (R 4.0.2)                      
#>  prettyunits     1.1.1      2020-01-24 [1] CRAN (R 4.0.0)                      
#>  processx        3.4.4      2020-09-03 [1] CRAN (R 4.0.2)                      
#>  ps              1.4.0      2020-10-07 [1] CRAN (R 4.0.2)                      
#>  purrr           0.3.4      2020-04-17 [1] CRAN (R 4.0.0)                      
#>  R6              2.5.0      2020-10-28 [1] CRAN (R 4.0.3)                      
#>  Rcpp            1.0.5      2020-07-06 [1] CRAN (R 4.0.2)                      
#>  RcppArmadillo   0.10.1.0.0 2020-10-20 [1] CRAN (R 4.0.3)                      
#>  remotes         2.2.0      2020-07-21 [1] CRAN (R 4.0.2)                      
#>  rex             1.2.0      2020-04-21 [1] CRAN (R 4.0.0)                      
#>  rlang           0.4.8      2020-10-08 [1] CRAN (R 4.0.3)                      
#>  rmarkdown       2.5        2020-10-21 [1] CRAN (R 4.0.3)                      
#>  rprojroot       1.3-2      2018-01-03 [1] CRAN (R 4.0.0)                      
#>  rstudioapi      0.11       2020-02-07 [1] CRAN (R 4.0.0)                      
#>  RxODE         * 0.9.2-1    2020-10-15 [1] CRAN (R 4.0.3)                      
#>  scales          1.1.1      2020-05-11 [1] CRAN (R 4.0.2)                      
#>  sessioninfo     1.1.1      2018-11-05 [1] CRAN (R 4.0.0)                      
#>  SnakeCharmR     1.0.8      2020-05-04 [1] Github (asieira/SnakeCharmR@26002e0)
#>  stringi         1.5.3      2020-09-09 [1] CRAN (R 4.0.2)                      
#>  stringr         1.4.0      2019-02-10 [1] CRAN (R 4.0.0)                      
#>  sys             3.4        2020-07-23 [1] CRAN (R 4.0.2)                      
#>  testthat        2.3.2      2020-03-02 [1] CRAN (R 4.0.0)                      
#>  tibble          3.0.4      2020-10-12 [1] CRAN (R 4.0.3)                      
#>  tidyselect      1.1.0      2020-05-11 [1] CRAN (R 4.0.2)                      
#>  units           0.6-7      2020-06-13 [1] CRAN (R 4.0.2)                      
#>  usethis         1.6.3      2020-09-17 [1] CRAN (R 4.0.2)                      
#>  utf8            1.1.4      2018-05-24 [1] CRAN (R 4.0.0)                      
#>  vctrs           0.3.4      2020-08-29 [1] CRAN (R 4.0.2)                      
#>  withr           2.3.0      2020-09-22 [1] CRAN (R 4.0.2)                      
#>  xfun            0.18       2020-09-29 [1] CRAN (R 4.0.2)                      
#>  yaml            2.2.1      2020-02-01 [1] CRAN (R 4.0.0)                      
#> 
#> [1] /usr/local/lib/R/site-library
#> [2] /usr/lib/R/site-library
#> [3] /usr/lib/R/library

mattfidler · 2020-10-29T18:07:24Z

There seems to be 2 different types of additive + proportional error models defined:

combined1

s2(v) = t1^2 + t2^2*v^2

And

combined2

s2(v) = (t1 + abs(v)^t2)^2

You are saying the patch referred to above allows combined1 style of additive + proportional errors?

mattfidler · 2020-10-29T18:08:49Z

These naming conventions come from monolix:

http://monolix.lixoft.com/data-and-models/errormodel/

mattfidler · 2020-10-29T18:11:19Z

In the next release you can switch between the two error models.

mattfidler · 2020-10-29T18:13:19Z

See https://github.com/nlmixrdevelopment/nlmixr/blob/foceiDur/NEWS.md

It would be great if we could switch in nlme as well.

mattfidler · 2020-10-29T18:19:28Z

The nlme variance is defined here:

nlmixr/R/ui.R

Lines 2826 to 2855 in bbab91e

    
           nlmixrUI.nlme.var <- function(object) { 
        
             ## Get the variance for the nlme object 
        
             add.prop.errs <- object$add.prop.errs 
        
             w.no.add <- which(!add.prop.errs$add) 
        
             w.no.prop <- which(!add.prop.errs$prop) 
        
             const <- grp <- "" 
        
             power <- ", fixed=c(1)" 
        
             powera <- ", fixed=list(power=1)" 
        
             if (length(add.prop.errs$y) > 1) { 
        
               grp <- " | nlmixr.grp" 
        
             } 
        
             if (length(w.no.add) > 0) { 
        
               const <- sprintf(", fixed=list(%s)", paste(paste0(add.prop.errs$y[w.no.add], "=0"), collapse = ", ")) 
        
             } 
        
             if (length(w.no.prop) > 0) { 
        
               power <- sprintf(", fixed=list(%s)", paste(paste0(add.prop.errs$y, "=", ifelse(add.prop.errs$prop, 1, 0)), collapse = ", ")) 
        
               powera <- sprintf(", fixed=list(power=list(%s))", paste(paste0(add.prop.errs$y, "=", ifelse(add.prop.errs$prop, 1, 0)), collapse = ", ")) 
        
             } 
        
             tmp <- sprintf("varConstPower(form=~fitted(.)%s%s)", grp, powera) 
        
             if (all(!add.prop.errs$prop)) { 
        
               tmp <- sprintf("varIdent(form = ~ 1%s)", grp) 
        
               if (tmp == "varIdent(form = ~ 1)") { 
        
                 warning("initial condition for additive error ignored with nlme") 
        
                 return(NULL) 
        
               } 
        
             } else if (all(!add.prop.errs$add)) { 
        
               tmp <- sprintf("varPower(form = ~ fitted(.)%s%s)", grp, power) 
        
             } 
        
             return(eval(parse(text = tmp))) 
        
           }

mattfidler · 2020-10-29T18:27:57Z

Your model above produces:

varConstPower(form=~fitted(.), fixed=list(power=1))

But in your comparison you use:

varConstProp()

I could add this when the error model is combined1; I believe the first is correct for combined2.

mattfidler · 2020-10-29T18:53:08Z

It seems that are two different schools of thoughts with combined1 and combined2. The next release defaults to combined2 which is more commonly used in NONMEM, but perhaps it should be a nlmixr option you can choose at the beginning as well. If you are a fan of the combined1 approach, it would be a pain to always have to specify it instead of combined2

jranke · 2020-10-29T18:57:41Z

There seems to be 2 different types of additive + proportional error models defined:

combined1
s2(v) = t1^2 + t2^2*v^2
And

combined2
s2(v) = (t1 + abs(v)^t2)^2
You are saying the patch refered to above allows combined1 style of additive + proportional errors?
Yes, indeed.

mattfidler · 2020-10-29T19:06:57Z

Ok @jranke

Let me know when the patch is accepted. When it is, let me know the version of nlme that is used so I an add combined1 as an option for nlme.

jranke · 2020-10-29T19:13:17Z

Sorry, I was too quick. The nlme patch provides

s2(v) = t1^2 + t2^2*v^2

which is equivalent to combined2 in Monolix terms. Note that s2(v) is on a variance scale. See also my post to R-devel
https://stat.ethz.ch/pipermail/r-devel/2020-October/080046.html and the references cited there. The patch is currently being reviewed by Martin Maechler - it would be great if it would make it into nlme - I will let you know!

mattfidler · 2020-10-29T21:34:54Z

So, misspoke:

combined1

s2(v) = (t1 + abs(v)^t2)^2

combined2

s2(v) = t1^2 + t2^2*v^2

It would good to have combined2 in nlmixr nlme, especially since it is the default for the other routines.

jranke · 2020-10-31T07:38:10Z

The patch to nlme was merged yesterday, yay! You can easily get it from their subversion repo if you have subversion installed:

svn co https://svn.r-project.org/R-packages/trunk/nlme

mattfidler · 2020-10-31T16:37:50Z

Great! I will look into using it for nlmixr too

jranke · 2020-10-31T16:52:56Z

Probably all you need to do is change varConstPow to varConstProp, fix sigma in the call to nlme ..., control = list(sigma = 1) and initialise const with the additive error and prop with the proportional error.... good luck!

mattfidler · 2020-11-19T15:37:48Z

This has been integrated; Do you know when nlme will be released (only with a new version of R)?

jranke · 2020-11-19T15:44:25Z

This has been integrated; Do you know when nlme will be released (only with a new version of R)?

Nice! nlme has its own release cycle, and has releases every couple of months, so it depends.

jranke · 2020-11-19T16:03:32Z

I think you need to check the formula in the news entry you wrote for this, as nlme::varConstProp() does not include parameter c in the error model

y = f + sqrt(a^2+b^2*f^(2c))*err

only the additive component a and the proportionality constant b.

mattfidler · 2020-11-19T16:21:00Z

Indeed you are right; I have an error issued if this occurs:

nlmixr/R/ui.R

Line 2864 in b7c7ccd

stop("add+prop combined2 does not support power for nlme currently")

I was testing and noticed that b can be negative, I will adjust this to be positive just for consistency between the methods.

Note you can always get the underlying nlme by as.nlme

mattfidler · 2020-11-19T16:24:55Z

I misread your comment, I will adjust the news item.

jranke · 2020-11-19T16:28:13Z

Ah, I see, but your explanation was nevertheless helpful as I hadn't looked at the code!

jranke · 2020-12-16T09:06:48Z

Hi, just in case you did not notice, nlme version 3.1.151 containing varConstProp has been published a couple of days ago.

mattfidler · 2020-12-16T12:49:27Z

Awesome! We will be publishing soon too so it should work perfectly.

…

On Wed, Dec 16, 2020, 3:07 AM Johannes Ranke ***@***.***> wrote: Hi, just in case you did not notice, nlme version 3.1.151 containing varConstProp has been published a couple of days ago. — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#428 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAD5VWQCWA7BEOPZOUNATNTSVB2DPANCNFSM4TD5SPXA> .

mattfidler added a commit that referenced this issue Nov 4, 2020

Start using combined2 for nlme See #428

9ec5931

mattfidler added a commit that referenced this issue Nov 19, 2020

Updated news for #428

c286ef3

mattfidler closed this as completed Nov 19, 2020

mattfidler added a commit that referenced this issue Nov 19, 2020

Add abs(prop) for nlme #428

cddd57f

mattfidler added a commit that referenced this issue Nov 19, 2020

Full name #428

c2c92f9

mattfidler added a commit that referenced this issue Nov 19, 2020

word-smith error #428

afe1d10

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

nlme backend does not work correctly with additive and proportional error #428

nlme backend does not work correctly with additive and proportional error #428

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020 •

edited

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020 •

edited

mattfidler commented Oct 29, 2020

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020

jranke commented Oct 31, 2020

mattfidler commented Oct 31, 2020

jranke commented Oct 31, 2020

mattfidler commented Nov 19, 2020

jranke commented Nov 19, 2020

jranke commented Nov 19, 2020

mattfidler commented Nov 19, 2020

mattfidler commented Nov 19, 2020

jranke commented Nov 19, 2020

jranke commented Dec 16, 2020

mattfidler commented Dec 16, 2020 via email

nlme backend does not work correctly with additive and proportional error #428

nlme backend does not work correctly with additive and proportional error #428

Comments

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020 • edited

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020

mattfidler commented Oct 29, 2020 • edited

mattfidler commented Oct 29, 2020

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020

jranke commented Oct 29, 2020

mattfidler commented Oct 29, 2020

jranke commented Oct 31, 2020

mattfidler commented Oct 31, 2020

jranke commented Oct 31, 2020

mattfidler commented Nov 19, 2020

jranke commented Nov 19, 2020

jranke commented Nov 19, 2020

mattfidler commented Nov 19, 2020

mattfidler commented Nov 19, 2020

jranke commented Nov 19, 2020

jranke commented Dec 16, 2020

mattfidler commented Dec 16, 2020 via email

mattfidler commented Oct 29, 2020 •

edited

mattfidler commented Oct 29, 2020 •

edited