nlme backend does not work correctly with additive and proportional error

[jranke]

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

``` r 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]

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]

These naming conventions come from monolix:

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

[mattfidler]

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

[mattfidler]

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

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

[mattfidler]

The nlme variance is defined here:

https://github.com/nlmixrdevelopment/nlmixr/blob/bbab91e55a1a5ede334be8b549f7439b161db260/R/ui.R#L2826-L2855

[mattfidler]

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]

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]

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]

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]

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]

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]

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]

Great! I will look into using it for nlmixr too

[jranke]

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]

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

[jranke]

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]

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]

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

https://github.com/nlmixrdevelopment/nlmixr/blob/b7c7ccd6b472b1647e72aa876c49672ffa7c76e5/R/ui.R#L2864

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]

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

[jranke]

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

[jranke]

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

[mattfidler]

Awesome!

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

On Wed, Dec 16, 2020, 3:07 AM Johannes Ranke notifications@github.com 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 https://github.com/nlmixrdevelopment/nlmixr/issues/428#issuecomment-745956215, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAD5VWQCWA7BEOPZOUNATNTSVB2DPANCNFSM4TD5SPXA .